Znajdź prawdziwe korzenie wielomianu

Napisz samodzielny program, który po otrzymaniu wielomianu i powiązania znajdzie wszystkie rzeczywiste pierwiastki tego wielomianu na błąd absolutny nieprzekraczający granicy.

Ograniczenia

Wiem, że Mathematica i prawdopodobnie niektóre inne języki mają rozwiązanie z jednym symbolem, co jest nudne, więc powinieneś trzymać się prymitywnych operacji (dodawanie, odejmowanie, mnożenie, dzielenie).

Istnieje pewna elastyczność formatów wejściowych i wyjściowych. Możesz przyjmować dane wejściowe za pomocą argumentów stdin lub wiersza poleceń w dowolnym rozsądnym formacie. Możesz zezwolić na liczbę zmiennoprzecinkową lub wymagać użycia pewnej reprezentacji liczb wymiernych. Możesz wziąć granicę lub odwrotność granicy, a jeśli używasz zmiennoprzecinkowego, możesz założyć, że granica nie będzie mniejsza niż 2 ulp. Wielomian powinien być wyrażony jako lista współczynników monomialnych, ale może być dużym lub małym endianem.

Musisz być w stanie uzasadnić, dlaczego twój program zawsze będzie działał (zagadnienia numeryczne modulo), chociaż nie jest konieczne dostarczanie pełnych dowodów w linii.

Program musi obsługiwać wielomiany z powtarzanymi pierwiastkami.

Przykład

x^2 - 2 = 0 (error bound 0.01)

Dane wejściowe mogą być np

-2 0 1 0.01
100 1 0 -2
1/100 ; x^2-2

Dane wyjściowe mogą być np

-1.41 1.42

ale nie

-1.40 1.40

ponieważ ma to błędy bezwzględne około 0,014 ...

Przypadki testowe

Prosty:

x^2 - 2 = 0 (error bound 0.01)

x^4 + 0.81 x^2 - 0.47 x + 0.06 (error bound 10^-6)

Wiele korzeni:

x^4 - 8 x^3 + 18 x^2 - 27 (error bound 10^-6)

Wielomian Wilkinsona:

x^20 - 210 x^19 + 20615 x^18 - 1256850 x^17 + 53327946 x^16 -1672280820 x^15 +
    40171771630 x^14 - 756111184500 x^13 + 11310276995381 x^12 - 135585182899530 x^11 +
    1307535010540395 x^10 - 10142299865511450 x^9 + 63030812099294896 x^8 -
    311333643161390640 x^7 + 1206647803780373360 x^6 -3599979517947607200 x^5 +
    8037811822645051776 x^4 - 12870931245150988800 x^3 + 13803759753640704000 x^2 -
    8752948036761600000 x + 2432902008176640000  (error bound 2^-32)

Uwaga: To pytanie było w piaskownicy przez około 3 miesiące. Jeśli uważasz, że wymagało to ulepszenia przed opublikowaniem, odwiedź Sandbox i skomentuj pozostałe proponowane pytania, zanim zostaną opublikowane w Main.

Mathematica, 223

r[p_,t_]:=Module[{l},l=Exponent[p[x],x];Re@Select[NestWhile[Table[#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}],{i,l}]&,Table[(0.9+0.1*I)^i,{i,l}],2*Max[Table[Abs[#1[[i]]-#2[[i]]],{i,l}]]>t&,2],Abs[Im[#]]<t^.5&]]

To rozwiązanie implementuje metodę Durand – Kerner do rozwiązywania wielomianów. Zauważ, że nie jest to kompletne rozwiązanie (jak zostanie pokazane poniżej), ponieważ nie mogę jeszcze obsłużyć wielomianu Wilkinsona z określoną precyzją. Najpierw wyjaśnienie tego, co robię:

#[[i]]-p[#[[i]]]/Product[If[i!=j,#[[i]]-#[[j]],1],{j,l}]&: W ten sposób funkcja oblicza dla każdego indeksu inastępne przybliżenie Duranda-Kernera. Następnie ta linia jest hermetyzowana w tabeli i stosowana za pomocą NestWhile do punktów wejściowych generowanych przez Table[(0.9+0.1*I)^i,{i,l}]. Warunkiem na NestWhile jest to, że maksymalna zmiana (we wszystkich warunkach) z jednej iteracji do następnej jest większa niż określona precyzja. Kiedy wszystkie warunki zmieniły się mniej niż to, NestWhile kończy i Re@Selectusuwa zera, które nie spadają na rzeczywistą linię.

Przykładowe dane wyjściowe:

> p[x_] := x^2 - 2
> r[p, .01]
{1.41421, -1.41421}

> p[x_] := x^4 - 8 x^3 + 18 x^2 - 27
> r[p, 10^-6]
{2.99999, 3., 3.00001, -1.}

> p[x_] := x^20 - 210 x^19 + ... + 2432902008176640000 (Wilkinson's)
> Sort[r[p, 2^-32]]
{1., 2., 3., 4., 5., 6., 7.00001, 7.99994, 9.00018, 10.0002, 11.0007, \
11.9809, 13.0043, 14.0227, 14.9878, 16.0158, 16.9959, 17.9992, \
19.0001, 20.}

Jak zapewne widać, gdy stopień rośnie, metoda ta zaczyna odbijać się od prawidłowych wartości, nigdy tak naprawdę nie bazując całkowicie. Jeśli ustawię warunek zatrzymania mojego kodu na coś bardziej rygorystycznego niż „od jednej iteracji do następnej domysły zmienione o nie więcej niż epsilon”, algorytm nigdy się nie zatrzyma. Chyba powinienem po prostu użyć Duranda-Kernera jako danych wejściowych do metody Newtona?