Welcome, guest | Sign In | My Account | Store | Cart
```# Durand-Kerner method for solving polynomial equations
# http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method
# http://en.wikipedia.org/wiki/Horner%27s_method
# http://en.wikipedia.org/wiki/Properties_of_polynomial_roots
# FB - 20130428
import random
import math

def bound(coefficients):
coefficients.reverse()
n = len(coefficients) - 1
b = 0.0
for i in range(n):
b += abs(coefficients[i] / coefficients[i + 1])
coefficients.reverse()
return b

def polynomial(z, coefficients): # Horner method
t = complex(0, 0)
for c in reversed(coefficients):
t = t * z + c
return t

eps = 1e-7 # max error allowed
def DurandKerner(coefficients):
n = len(coefficients) - 1
roots = [complex(0, 0)] * n
bnd = bound(coefficients)
retry = True
while retry:
retry = False
# set initial roots as random points within bounding circle
for k in range(n):
r = bnd * random.random()
theta = 2.0 * math.pi * random.random()
roots[k] = complex(r * math.cos(theta), r * math.sin(theta))

itCtr = 0
rootsNew = roots[:]
flag = True
while flag:
flag = False
for k in range(n):
temp = complex(1.0, 0.0)
for j in range(n):
if j != k:
temp *= roots[k] - roots[j]
rootsNew[k] = roots[k] - polynomial(roots[k], coefficients) / temp
if abs(roots[k] - rootsNew[k]) > eps:
# print abs(roots[k] - rootsNew[k])
flag = True
if math.isnan(rootsNew[k].real) or math.isnan(rootsNew[k].imag):
flag = False
retry = True
print 'retrying...'
break
roots = rootsNew[:]
itCtr += 1

print "iteration count: " + str(itCtr)
return roots

# example
# x**3-3*x**2+3*x-5=0
coefficients = [complex(-5, 0), complex(3, 0), complex(-3, 0), complex(1, 0)]
print "coefficients: " + str(coefficients)
print "roots: " + str(DurandKerner(coefficients))
```

#### Diff to Previous Revision

```--- revision 1 2011-09-05 01:40:40
+++ revision 2 2013-04-28 20:47:23
@@ -1,35 +1,67 @@
# Durand-Kerner method for solving polynomial equations
# http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method
-# FB - 201109046
+# http://en.wikipedia.org/wiki/Horner%27s_method
+# http://en.wikipedia.org/wiki/Properties_of_polynomial_roots
+# FB - 20130428
import random
+import math

-def polynomial(z, coefficients):
+def bound(coefficients):
+    coefficients.reverse()
+    n = len(coefficients) - 1
+    b = 0.0
+    for i in range(n):
+        b += abs(coefficients[i] / coefficients[i + 1])
+    coefficients.reverse()
+    return b
+
+def polynomial(z, coefficients): # Horner method
t = complex(0, 0)
-    for k in range(len(coefficients)):
-        t += coefficients[k] * z ** k
+    for c in reversed(coefficients):
+        t = t * z + c
return t

eps = 1e-7 # max error allowed
def DurandKerner(coefficients):
n = len(coefficients) - 1
roots = [complex(0, 0)] * n
-    for k in range(n):
-        roots[k] = complex(random.random(), random.random())
-    rootsNew = roots[:]
-    flag = True
-    while flag:
-        flag = False
+    bnd = bound(coefficients)
+    retry = True
+    while retry:
+        retry = False
+        # set initial roots as random points within bounding circle
for k in range(n):
-            temp = complex(1, 0)
-            for j in range(n):
-                if j != k:
-                    temp *= roots[k] - roots[j]
-            rootsNew[k] = roots[k] - polynomial(roots[k], coefficients) / temp
-            if abs(roots[k] - rootsNew[k]) > eps:
-                flag = True
-        roots = rootsNew[:]
+            r = bnd * random.random()
+            theta = 2.0 * math.pi * random.random()
+            roots[k] = complex(r * math.cos(theta), r * math.sin(theta))
+
+        itCtr = 0
+        rootsNew = roots[:]
+        flag = True
+        while flag:
+            flag = False
+            for k in range(n):
+                temp = complex(1.0, 0.0)
+                for j in range(n):
+                    if j != k:
+                        temp *= roots[k] - roots[j]
+                rootsNew[k] = roots[k] - polynomial(roots[k], coefficients) / temp
+                if abs(roots[k] - rootsNew[k]) > eps:
+                    # print abs(roots[k] - rootsNew[k])
+                    flag = True
+                if math.isnan(rootsNew[k].real) or math.isnan(rootsNew[k].imag):
+                    flag = False
+                    retry = True
+                    print 'retrying...'
+                    break
+            roots = rootsNew[:]
+            itCtr += 1
+
+    print "iteration count: " + str(itCtr)
return roots

-# main
-coefficients = [complex(-5, 0), complex(3, 0),complex(-3, 0),complex(1, 0)]
-print DurandKerner(coefficients)
+# example
+# x**3-3*x**2+3*x-5=0
+coefficients = [complex(-5, 0), complex(3, 0), complex(-3, 0), complex(1, 0)]
+print "coefficients: " + str(coefficients)
+print "roots: " + str(DurandKerner(coefficients))
```