Znajdź prawdziwe korzenie wielomianu

24

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.

code-golf math polynomials Peter Taylor
źródło
„Musisz być w stanie uzasadnić, dlaczego twój program zawsze będzie działał” Na wszelki wypadek, gdy ktoś potrzebuje pomocy
dr belisarius
@belisarius, ??
Peter Taylor
3
miał być żartem :(
Dr. Belisarius
Wiem, że to stare wyzwanie, więc nie czuję się zobowiązana do odpowiedzi, jeśli nie masz ochoty na ponowne otwarcie. (a) Czy możemy napisać funkcję, czy tylko pełny program? (b) W przypadku, gdy możemy napisać funkcję, czy możemy założyć, że dane wejściowe wykorzystują jakiś wygodny typ danych, np. Python fractions.Fraction(typ racjonalny)? (c) Czy musimy radzić sobie z wielomianami stopnia <1? (d) Czy możemy założyć, że wiodącym współczynnikiem jest 1?
Ell
(e) W odniesieniu do wielomianów z powtarzającymi się pierwiastkami, warto wprowadzić rozróżnienie między pierwiastkami nieparzystych i parzystych wielokrotności (przypadki testowe mają tylko pierwiastki nieparzystych wielokrotności.) Chociaż pierwiastki nieparzystej wielokrotności nie są zbyt trudne do pokonania, ja ' Nie jestem pewien, jak namacalne jest prawidłowe posługiwanie się liczbami pierwiastkowymi nawet wielokrotności, zwłaszcza, że ​​określa się margines błędu dla wartości pierwiastków, a nie dla ich istnienia. (...)
Ell

Odpowiedzi:

8

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ę: kod w formacie matematycznym

#[[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?

Kaya
źródło
Durand-Kerner ma również potencjalne problemy z wieloma źródłami. (Metoda Newtona też może niewiele pomóc - wielomian Wilkinsona jest specjalnie wybierany jako źle uwarunkowany).
Peter Taylor
Masz całkowitą rację: zrezygnowałem z tego działania po powiększeniu przybliżenia Wilkinsona do x = 17, to absolutny bałagan. Martwię się, że będę musiał wybrać symboliczne rozwiązanie oparte na Groebner, aby uzyskać znacznie większą dokładność.
Kaya