# Attempt at solving roots of polynomial

## Recommended Posts

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

```#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()
\$array = Stringsplit(\$userinput, "x")
Select
Case StringinStr(\$array[2], "^2") = 1
second()
Case StringinStr(\$array[2], "^3") = 1
third()
EndSelect
EndFunc

Func second()
\$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()
\$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))))
\$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))))
\$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))))
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.

```#include <Math.au3>
\$pi = 3.141592653589793

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 = 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 by Ultima2

##### Share on other sites

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 by Wus

Thank you, Wus

## Create an account

Register a new account

×

• Wiki

• Back

• Git