多項式の零点(2)

零点が存在する円の半径Rに対して、複素平面上で原点を中心とした-RからRの正方形を考えて、それを縦横10分割した格子点から最急降下法で零点を求める。
これを、f(x) = 16-52x+48x2+68x3-30x4-39x5+7x6+4x7-10x8-2x9+6x10+x11-x12 に対して適用したところ、次のような12点が求められた。


(2.45095897081+5.69021343773e-039j)
(-1.74686952669+2.64116131885e-035j)
(1.76049219702+7.70659973871e-033j)
(1.18961203019-4.36843759502e-018j)
(0.351601004487-0.259979952645j)
(0.351601004487+0.259979952645j)
(-1.3263717523-0.491161094885j)
(-1.3263717523+0.491161094885j)
(-0.940229729687-1.11637690587j)
(-0.940229729687+1.11637690587j)
(0.587903641832-1.35795913946j)
(0.587903641832+1.35795913946j)



from cmath import *
from poly import Poly

def Newton(f, x):
eps = 1e-9;
g = f.differential();
while abs(f.value(x)) > eps:
x -= f.value(x) / g.value(x)
return x

def convert(f):
p = [ -abs(x) for x in f ]
p[len(p)-1] *= -1;
return Poly(p)

def calcRadius(h):
h = convert(f)
x = 1e-3
while h.value(x) < 0:
x += 1e-3
return Newton(h, x)

def appendZero(v, f, z):
for i in range(0, len(v)):
x = v[i]
dist = abs(x.real - z.real) + abs(x.imag - z.imag)
if dist < 1e-8:
if abs(f.value(z)) < abs(f.value(x)):
v[i] = z
return
v.append(z)

f = Poly([ 16, -52, 48, 68, -30, -39, 7, 4, -10, -2, 6, 1, -1 ]);

R = calcRadius(f) # radius of solutions
zero_list = [ ]

y = -R
while y < R + 1e-3:
x = -R
while x < R + 1e-3:
z = x + 1j * y
z = Newton(f, z)
appendZero(zero_list, f, z)
x += R / 5
y += R / 5

zero_list.sort(key = lambda x: abs(x.imag))
for z in zero_list:
print z