Реализация Durand-kerner не работает

Question

Что случилось с этой реализацией алгоритма Дюран-Кернера (here)?

def durand_kerner(poly, start=complex(.4, .9), epsilon=10**-16):#float('-inf')): 
    roots = [] 
    for e in xrange(poly.degree): 
     roots.append(start ** e) 
    while True: 
     new = [] 
     for i, r in enumerate(roots): 
      new_r = r - (poly(r))/(reduce(operator.mul, [(r - r_1) for j, r_1 in enumerate(roots) if i != j])) 
      new.append(new_r) 
     if all(n == roots[i] or abs(n - roots[i]) < epsilon for i, n in enumerate(new)): 
      return new 
     roots = new 

Когда я пытаюсь это, я должен остановить его KeyboardInterrupt, потому что это не мешает!
poly - это полиномиальный экземпляр библиотеки pypol.

Спасибо заранее, рубик

EDIT: Используя NumPy полином он занимает 9 итераций:

In [1]: import numpy as np 

In [2]: roots.d1(np.poly1d([1, -3, 3, -5])) 
3 
[(1.3607734793516519+2.0222302921553128j), (-1.3982133295376746-0.69356635962504309j), (3.0374398501860234-1.3286639325302696j)] 
[(0.98096328371966801+1.3474626910848715j), (-0.3352519326012724-0.64406860772816388j), (2.3542886488816044-0.70339408335670761j)] 
[(0.31718054925650596+0.93649454851955749j), (0.49001572078718736-0.9661410790307261j), (2.1928037299563066+0.029646530511168612j)] 
[(0.20901563897345796+1.5727420147652911j), (0.041206038662691125-1.5275192097633465j), (2.7497783223638508-0.045222805001944255j)] 
[(0.21297050700971876+1.3948274731404162j), (0.18467846583682396-1.3845653821841168j), (2.6023510271534573-0.010262090956299326j)] 
[(0.20653075193800668+1.374878742771485j), (0.20600107336130213-1.3746529207714699j), (2.5874681747006911-0.00022582200001499547j)] 
[(0.20629950692533283+1.3747296033941407j), (0.20629947661265013-1.374729584400741j), (2.5874010164620169-1.899339978055233e-08j)] 
[(0.20629947401589896+1.3747296369986031j), (0.20629947401590082-1.3747296369986042j), (2.5874010519682002+9.1830687539942581e-16j)] 
[(0.20629947401590029+1.3747296369986026j), (0.20629947401590026-1.3747296369986026j), (2.5874010519681994+1.1832913578315177e-30j)] 
Out[2]: 
[(0.20629947401590029+1.3747296369986026j), 
(0.20629947401590029-1.3747296369986026j), 
(2.5874010519681994+0j)] 

Используя pypol полином он никогда не заканчивает (это, вероятно, ошибка в pypol):

In [3]: roots.d2(poly1d([1, -3, 3, -5])) 
^C--------------------------------------------------------------------------- 
KeyboardInterrupt 

, но я не могу найти ошибку !!

EDIT2: Сравнивая метод __call__ с Poly Мартина:

>>> p = Poly(-5, 3, -3, 1) 
>>> from pypol import poly1d 
>>> p2 = poly1d([1, -3, 3, -5]) 

>>> for i in xrange(-100000, 100000): 
    assert p(i) == p2(i) 


>>> 
>>> for i in xrange(-10000, 10000): 
    assert p(complex(1, i)) == p2(complex(1, i)) 


>>> for i in xrange(-10000, 10000): 
    assert p(complex(i, i)) == p2(complex(i, i)) 


>>> 

EDIT3: pypol отлично работает, если корни не являются комплексными числами:

In [1]: p = pypol.funcs.from_roots([4, -2, 443, -11212]) 

In [2]: durand_kerner(p) 
Out[2]: [(4+0j), (443+0j), (-2+0j), (-11212+0j)] 

Так что Безразлично» t работает только тогда, когда корни являются комплексными числами!

EDIT4: Я написал несколько иную реализацию для Numpy полиномов и увидел, что после одной итерации корни (полинома Википедии) различны:

In [4]: d1(numpyp.poly1d([1, -3, 3, -5])) 
Out[4]: 
[(0.98096328371966801+1.3474626910848715j), 
(-0.3352519326012724-0.64406860772816388j), 
(2.3542886488816044-0.70339408335670761j)] 

In [5]: d2(pypol.poly1d([1, -3, 3, -5])) 
Out[5]: 
[(0.9809632837196679+1.3474626910848717j), 
(-0.33525193260127306-0.64406860772816377j), 
(2.3542886488816048-0.70339408335670772j)] ## here 

EDIT5: Эй! Если я меняю линию: if all(n == roots[i] ...) на if all(str(n) == str(roots[i]) ...), она заканчивается и возвращает правильные корни !!!

In [9]: p = pypol.poly1d([1, -3, 3, -5]) 

In [10]: roots.durand_kerner(p) 
Out[10]: 
[(0.20629947401590029+1.3747296369986026j), 
(0.20629947401590013-1.3747296369986026j), 
(2.5874010519681994+0j)] 

Но вопрос в том, почему он работает с другим комплексом комплексных чисел?

UPDATE
Теперь он работает, и я сделал несколько тестов:

In [1]: p = pypol.poly1d([1, -3, 3, -1]) 

In [2]: p 
Out[2]: + x^3 - 3x^2 + 3x - 1 

In [3]: pypol.roots.cubic(p) 
Out[3]: (1.0, 1.0, 1.0) 

In [4]: durand_kerner(p) 
Out[4]: 
((1+0j), 
(1.0000002484566535-2.708692281244913e-17j), 
(0.99999975147728026+2.9792265510301965e-17j)) 

In [5]: q = x ** 3 - 1 

In [6]: q 
Out[6]: + x^3 - 1 

In [7]: pypol.roots.cubic(q) 
Out[7]: (1.0, (-0.5+0.8660254037844386j), (-0.5-0.8660254037844386j)) 

In [8]: durand_kerner(q) 
Out[8]: ((1+0j), (-0.5-0.8660254037844386j), (-0.5+0.8660254037844386j)) 

Answer 1

Ваш алгоритм выглядит хорошо, и это работает для меня пример в википедии

import operator 
class Poly: 
    def __init__(self, *koeff): 
     self.koeff = koeff 
     self.degree = len(koeff)-1 

    def __call__(self, val): 
     res = 0 
     x = 1 
     for k in self.koeff: 
      res += x*k 
      x *= val 
     return res 

def durand_kerner(poly, start=complex(.4, .9), epsilon=10**-16):#float('-inf')): 
    roots = [] 
    for e in xrange(poly.degree): 
     roots.append(start ** e) 
    while True: 
     new = [] 
     for i, r in enumerate(roots): 
      new_r = r - (poly(r))/(reduce(operator.mul, [(r - r_1) 
            for j, r_1 in enumerate(roots) if i != j])) 
      new.append(new_r) 
     if all((n == roots[i] or abs(n - roots[i]) < epsilon) for i, n in enumerate(new)): 
      return new 
     roots = new 

print durand_kerner(Poly(-5,3,-3,1)) 

[(0.20629947401590026+1.3747296369986026j), 
(0.20629947401590026-1.3747296369986026j), 
(2.5874010519681994+8.6361685550944446e-78j)] 

Answer 2

О вашем «EDIT 5»: это происходит потому, что str() не форматирует номера до полной точности.

>>> print str((2.5874010519681994+8.636168555094445e-78j)) 
(2.58740105197+8.63616855509e-78j) 
>>> print repr((2.5874010519681994+8.636168555094445e-78j)) 
(2.5874010519681994+8.636168555094445e-78j) 
>>> 

Так что не делайте этого.

В любом случае, проверка на равенство в вашем коде:

if all(n == roots[i] or abs(n - roots[i]) < epsilon for i, n in enumerate(new)): 

является излишним; если n == roots[i], то abs(n - roots[i]) будет равна нулю, так что вы могли бы сделать только

if all(abs(n - roots[i]) < epsilon for i, n in enumerate(new)): 

и положить немного усилий в разработке, что по умолчанию для epsilon должно быть; как я указал в комментарии, решение X**3 == 1 сходится, но ваш эпсилон по умолчанию слишком мал, чтобы остановить его навсегда. 1.12e-16 выглядит как лучшая ставка для эпсилона по умолчанию.

Для развлечения попробуйте неподходящий Poly ([- 1, 3, -3, 1]) ... три корня все равны (1 + 0j) ... он принимает более 600 итераций, а ошибки в последние 10 или около того итераций идут удивительно, пока он не придет к разумному решению из дальнего левого поля.