Attempt at solving roots of polynomial

[Ultima2]

Here is my attempt at writing a script that will solve any roots. It uses the Aberth method

http://en.wikipedia.org/wiki/Aberth_method And it can even solve complex roots. It only goes up to the third degree for right now but I will improve it.

Here is the code

expandcollapse popup

#include <GUIConstantsEx.au3>
#include "complex.au3"

$gui = GUIcreate("Root Solver", 900, 500)
$exit = GUICtrlCreateButton("EXIT", 100, 100, 50, 50)
$input = GUICtrlCreateinput("",20, 50, 220, 30)
$ok = GUICtrlCreatebutton("OK", 30, 100, 50, 50)
$list = GUICtrlCreateList("", 280, 50, 380, 400)
$list2 = GUICtrlCreateList("",680, 50, 180, 400)
GUICtrlCreateLabel("Iteration Trail",680, 25)
GUICtrlCreateLabel("Put in standard form, put 0x^2 if there is no x^2",15, 22)

GUIsetstate(@SW_SHOW)

While 1
    $msg = GUIgetmsg()
    Select
    Case $msg = $exit
        Exit
    Case $msg = $GUI_EVENT_CLOSE 
        Exit
    Case $msg = $ok
        GUICtrlSetData($list, "")
        GUICtrlSetData($list2, "")
        polynomial()
    Endselect
WEnd

Func polynomial()
    $userinput = GUICtrlRead($input)
    $array = Stringsplit($userinput, "x")
    Select 
    Case StringinStr($array[2], "^2") = 1 
    second()
    Case StringinStr($array[2], "^3") = 1
    third()
    EndSelect
EndFunc

Func second()
    $userinput = GUICtrlRead($input)
    $array = Stringsplit($userinput, "x")
    For $i = 1 to $array[0]
        $array[$i] = Stringstripws($array[$i], 8)
    next

    If Stringinstr($array[1], "-") = 1 Then
        $array[1] = (-1)*(StringTrimleft($array[1], 1))
        If StringTrimLeft($array[1], 1) = "0" Then
            $array[1] = -1
        EndIf
    EndIf
    
    If StringLen($array[1]) = 0 Then
        $array[1] = 1
    EndIf
    
    If StringInStr($array[2], "-") > 0 Then
        $b = (-1)*(StringTrimLeft($array[2], 3))
        If StringTrimleft($b, 1) = "0" Then
            $b = -1
        EndIf
    ElseIf StringinStr($array[2], "+") > 0 Then
        $b = StringTrimLeft($array[2], 3)
    EndIf
    
    If Stringlen($b) = 0 Then
        $b = 1
    EndIf
    
    If StringinStr($array[3], "+") = 1 Then
        $c = StringtrimLeft($array[3], 1) + 0.0
    ElseIf StringinStr($array[3], "-") = 1 Then
        $c = (-1)*(StringTrimLeft($array[3], 1))
    EndIf
    
    If StringTrimLeft($array[3], 1) = 0 Then
        $c = $array[3]
    EndIf
    
    $a = complex($array[1], 0)
    $b = complex($b, 0)
    $c = complex($c, 0)
    $two = complex(2,0)
    $l = complex(4, 0)
    $i = complex(-1,0)
    $g = complex_pow($b, 2)
    $l = complex_mul(complex_mul($a, $c), $l)
    $g = complex_pow(complex_minus($g, $l), 0.5)
    $b = complex_mul($b, $i)
    $x1 = complex_div(complex_add($b, $g), complex_mul($a, $two))
    $x2 = complex_div(complex_minus($b, $g), complex_mul($a, $two))
    GUICtrlSetData($list, "x = "&$x1)
    GUICtrlSetData($list, "x = "&$x2)
EndFunc

;uses Aberth Method, it is much better than Newton or Durand-Kerner
Func third()
    $userinput = GUICtrlRead($input)
    $array = Stringsplit($userinput, "x")
    For $i = 1 to $array[0]
        $array[$i] = Stringstripws($array[$i], 8)
    next
    
    For $i = 2 to $array[0] - 1
    If StringInStr($array[$i], "-") > 0 Then
        $array[$i] = (-1)*(StringTrimLeft($array[$i], 3))
        If StringTrimleft($array[$i], 1) = "0" Then
            $array[$i] = -1
        EndIf
    ElseIf StringinStr($array[$i], "+") > 0 Then
        $array[$i] = StringTrimLeft($array[$i], 3)
    EndIf
    Next

    For $i = 2 to $array[0] - 1
    If Stringlen($array[$i]) = 0 Then
        $array[$i] = 1
    EndIf
    Next
    
    For $i = 1 to $array[0]
    If StringLen($array[$i]) = 0 Then
        $array[$i] = 1
    EndIf
    Next
    
    For $i = 1 to $array[0] step 3
    If StringInStr($array[$i], "-") > 0 Then
        $array[$i] = (-1)*(StringTrimLeft($array[$i], 1))
        If StringTrimLeft($array[$i], 1) = "0" Then
            $array[$i] = -1
        EndIf
    EndIf
    Next
    
    $z = complex(Random(), Random()); Must get better seeds
    $z1 = complex(Random(), Random())
    $z2 = complex(Random(), Random())
;Cannot change the coefficients to negative
    $a = complex($array[1], 0)
    $b = complex($array[2], 0)
    $c = complex($array[3], 0)
    $d = complex($array[4], 0)
    $one = complex(1, 0)
    $l = complex(-1, 0)
    $i = 1
    While abs(complex_extract_re(complex_add(complex_mul($a, complex_pow($z, 3)),complex_mul($b, complex_pow($z, 2)),complex_mul($c, $z), $d))) > 1e-10
        $v = complex_add(complex_div($one, complex_minus($z, $z1)), complex_div($one, complex_minus($z, $z2)))
        $m = complex_div(complex_add(complex_mul($a, complex_pow($z, 3)),complex_mul($b, complex_pow($z, 2)),complex_mul($c, $z), $d) , complex_add(complex_mul(complex(3, 0),$a, complex_pow($z, 2)), complex_mul(complex(2, 0),$b,$z), $c))
        $u = complex_mul($l, complex_div($m, complex_minus($one, complex_mul($m, $v))))
        $z = complex_add($z, $u);
        $v = complex_add(complex_div($one, complex_minus($z1, $z)), complex_div($one, complex_minus($z1, $z2)))
        $m = complex_div(complex_add(complex_mul($a, complex_pow($z1, 3)),complex_mul($b, complex_pow($z1, 2)),complex_mul($c, $z1), $d) , complex_add(complex_mul(complex(3, 0),$a, complex_pow($z1, 2)), complex_mul(complex(2, 0),$b,$z1), $c))
        $u = complex_mul($l, complex_div($m, complex_minus($one, complex_mul($m, $v))))
        $z1 = complex_add($z1, $u);
        $v = complex_add(complex_div($one, complex_minus($z2, $z1)), complex_div($one, complex_minus($z2, $z)))
        $m = complex_div(complex_add(complex_mul($a, complex_pow($z2, 3)),complex_mul($b, complex_pow($z2, 2)),complex_mul($c, $z2), $d) , complex_add(complex_mul(complex(3, 0),$a, complex_pow($z2, 2)), complex_mul(complex(2, 0),$b,$z2), $c))
        $u = complex_mul($l, complex_div($m, complex_minus($one, complex_mul($m, $v))))
        $z2 = complex_add($z2, $u)
        GUICtrlSetData($list2, "Loop "&$i&" = "&$z&", "&$z1&", "&$z2)
        $i = $i + 1
        If abs(complex_extract_re(complex_add(complex_mul($a, complex_pow($z, 3)),complex_mul($b, complex_pow($z, 2)),complex_mul($c, $z), $d))) < 1e-10 or abs(complex_extract_re(complex_add(complex_mul($a, complex_pow($z1, 3)),complex_mul($b, complex_pow($z1, 2)),complex_mul($c, $z1), $d))) < 1e-10 or abs(complex_extract_re(complex_add(complex_mul($a, complex_pow($z2, 3)),complex_mul($b, complex_pow($z2, 2)),complex_mul($c, $z2), $d))) < 1e-10 Then
            ExitLoop
        EndIf
    Wend
        
    GUICtrlSetData($list, "x1 = "&$z)
    GUICtrlSetData($list, "x2 = "&$z1)
    GUICtrlSetData($list, "x3 = "&$z2)
EndFunc

Here is the code for the complex number include file. It is my own UDF that I created because I couldnt find any complex numbers UDFs.

expandcollapse popup

#include <Math.au3>
$pi = 3.141592653589793
$toRad = $pi/180

Global const $ni = complex(0,0)
Global const $gi = complex(1,0)

Func complex($re, $im)
    $complex = ""&$re&","&$im&"i"
    return $complex
EndFunc

;returns the magnitude and the argument(theta). This will be used a lot in conjunction with other
;functions
Func complex_polar($x)
    $x = StringSplit($x, ","); Splits the complex to a, bi
    $x[2] = StringTrimright($x[2], 1); Strips the i off of the imaginary number
    $x[1] = $x[1] + 0.0; needs to add 0.0 to make it a float for squaring, type conversion
    $x[2] = $x[2] + 0.0
    $mag = sqrt($x[1]^2 + $x[2]^2)
    $theta = (atan(abs($x[2])/abs($x[1])))*(180/$pi)
;this is to select the quadrant for the magnitude and theta. 3 -7i would be in the 4th quadrant
    Select
    Case $x[1] > 0 and $x[2] >= 0
        $theta = $theta
    Case $x[1] < 0 and $x[2] >= 0
        $theta = 180 - $theta
    Case $x[1] < 0 and $x[2] <= 0
        $theta = 180 + $theta
    Case $x[1] > 0 and $x[2] <= 0
        $theta = 360 - $theta
    EndSelect
    If $theta <> $theta Then
        $theta = 0
    EndIf
    $complex = ""&$mag&","&$theta&""
    Return $complex
EndFunc

;when this function is used with polar, you can do a lot of things
Func complex_rect($complex)
    $complex = StringSplit($complex, ",")
    $real = $complex[1]*cos($toRad*($complex[2]))
    $imag = $complex[1]*sin($toRad*($complex[2]))
    $real = Round($real, 12); rounds the number, can change if you want more precision
    $imag = Round($imag, 12)
    If abs($real) < 1e-10 Then; more ifs and thens for more precision handling >_>
        $real = 0
    EndIf
    If abs($imag) < 1e-10 Then
        $imag = 0
    EndIf
    $complex = ""&$real&","&$imag&"i"
    Return $complex
EndFunc

;this is an example of using the function polar and rect together.
;this can raise to a complex number to any power, even to float like 2.3
Func complex_pow($complex, $pow)
    $complex = complex_polar($complex)
    $complex = StringSplit($complex, ",")
    $complex[1] = $complex[1]^$pow
    $complex[2] = $complex[2]*$pow
    $complex = ""&$complex[1]&","&$complex[2]&""
    $complex = complex_rect($complex)
    Return $complex
EndFunc

Func complex_conj($complex1)
    $complex1 = StringSplit($complex1, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $complex1[2] = (-1)*$complex1[2] + 0.0
    $real = $complex1[1]
    $imag = $complex1[2]
    $complex = ""&$real&","&$imag&"i"
    Return $complex
EndFunc

Func complex_add($complex1, $complex2, $complex3 = $ni, $complex4 = $ni, $complex5 = $ni)
    $complex1 = StringSplit($complex1, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $complex2 = StringSplit($complex2, ",")
    $complex2[2] = StringTrimRight($complex2[2], 1)
    $complex3 = StringSplit($complex3, ",")
    $complex3[2] = StringTrimRight($complex3[2], 1)
    $complex4 = StringSplit($complex4, ",")
    $complex4[2] = StringTrimRight($complex4[2], 1)
    $complex5 = StringSplit($complex5, ",")
    $complex5[2] = StringTrimRight($complex5[2], 1)
    $real = $complex1[1] + $complex2[1] + $complex3[1] + $complex4[1] + $complex5[1]
    $imag = $complex1[2] + $complex2[2] + $complex3[2] + $complex4[2] + $complex5[2]
    $complex = ""&$real&","&$imag&"i"
    Return $complex
EndFunc

Func complex_minus($complex1, $complex2, $complex3 = $ni, $complex4 = $ni, $complex5 = $ni)
    $complex1 = StringSplit($complex1, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $complex2 = StringSplit($complex2, ",")
    $complex2[2] = StringTrimRight($complex2[2], 1)
    $complex3 = StringSplit($complex3, ",")
    $complex3[2] = StringTrimRight($complex3[2], 1)
    $complex4 = StringSplit($complex4, ",")
    $complex4[2] = StringTrimRight($complex4[2], 1)
    $complex5 = StringSplit($complex5, ",")
    $complex5[2] = StringTrimRight($complex5[2], 1)
    $real = $complex1[1] - $complex2[1] - $complex3[1] - $complex4[1] - $complex5[1]
    $imag = $complex1[2] - $complex2[2] - $complex3[2] - $complex4[2] - $complex5[2]
    $complex = ""&$real&","&$imag&"i"
    Return $complex
EndFunc

Func complex_mul($complex1, $complex2, $complex3 = $gi, $complex4 = $gi, $complex5 = $gi)
    $complex1 = complex_polar($complex1)
    $complex2 = complex_polar($complex2)
    $complex3 = complex_polar($complex3)
    $complex4 = complex_polar($complex4)
    $complex5 = complex_polar($complex5)
    $complex1 = StringSplit($complex1, ",")
    $complex2 = StringSplit($complex2, ",")
    $complex3 = StringSplit($complex3, ",")
    $complex4 = StringSplit($complex4, ",")
    $complex5 = StringSplit($complex5, ",")
    $mag = $complex1[1]*$complex2[1]*$complex3[1]*$complex4[1]*$complex5[1]
    $theta = $complex1[2] + $complex2[2] + $complex3[2] + $complex4[2] + $complex5[2]
    $complex = complex_rect(""&$mag&","&$theta&"")
    Return $complex
EndFunc

Func complex_div($complex1, $complex2)
    $complex3 = complex_conj($complex2)
    $complex1 = complex_mul($complex1, $complex3)
    $complex2 = complex_mul($complex2, $complex3)
    $complex1 = StringSplit($complex1, ",")
    $complex2 = StringSplit($complex2, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $real = $complex1[1]/$complex2[1]
    $imag = $complex1[2]/$complex2[1]
    $complex = ""&$real&","&$imag&"i"
    Return $complex
EndFunc

Func complex_abs($complex1)
    $complex1 = StringSplit($complex1, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $complex1[1] = $complex1[1] + 0.0
    $complex1[2] = $complex1[2] + 0.0
    $mag = sqrt($complex1[1]^2 + $complex1[2]^2)
    Return $mag
EndFunc

Func complex_extract_re($complex1)
    $complex1 = StringSplit($complex1, ",")
    $complex1[1] = $complex1[1] + 0.0
    Return $complex1[1]
EndFunc

Func complex_extract_im($complex1)
    $complex1 = StringSplit($complex1, ",")
    $complex1[2] = StringTrimRight($complex1[2], 1)
    $complex1[2] = $complex1[2] + 0.0
    Return $complex1[2]
EndFunc

Dont be too harsh, it is my first attempt at writing complex numbers UDF. You need to write x^3 -2 as x^3 +0x^2 +0x -2. Btw can anyone show me the optimizing process or equation for choosing the seeds for the Aberth method? Right now the seeds are chosen randomly. Thanks

Edited December 21, 2008 by Ultima2

[Wus]

So as a math major who's taken numerical analysis I enjoyed seeing this-- not much non-trivial math appears on these forums at all. That said, I feel like you'd have an easier time of this if you either (1) Used a language that inherently understood imaginary numbers (i.e. matlab) or (2) Used an object oriented language where you could create Polynomial and Complex classes. It's not that I think that what your doing is at all bad, just that your using AutoIT for something that it really wasn't designed for. You wouldn't use a tack hammer to drive in a 16d nail, so I wouldn't suggest using AutoIT for doing complicated numerical work.

That said I will PM you the original AMS article this approximation was published in. And I note here that this scheme is cubically convergent... which is pretty damn good. Makes me wonder why we didn't cover it in class.

Edited December 21, 2008 by Wus

[Ultima2]

Thank you, Wus