Jeśli szukasz dobrego ograniczenia błędu zaokrąglania, niekoniecznie potrzebujesz biblioteki o precyzji aribtrary. Zamiast tego możesz użyć uruchomionej analizy błędów.
Nie udało mi się znaleźć dobrego odnośnika online, ale wszystko to opisano w rozdziale 3.3 książki Nicka Highama „Dokładność i stabilność algorytmów numerycznych”. Pomysł jest dość prosty:
- Ponownie wprowadź współczynnik do kodu, aby w każdym wierszu było jedno przypisanie jednej operacji arytmetycznej.
- Dla każdej zmiennej np.
xUtwórz zmienną, x_errktóra jest inicjalizowana na zero, gdy xzostanie przypisana stała.
- Dla każdej operacji, np.
z = x * yZaktualizuj zmienną z_errprzy użyciu standardowego modelu arytmetyki zmiennoprzecinkowej oraz wynikających zz niej błędów błędów x_erri y_err.
- Zwracana wartość twojej funkcji powinna wtedy również mieć przypisaną odpowiednią
_errwartość. Jest to zależne od danych ograniczenie całkowitego błędu zaokrąglenia.
Trudna część to krok 3. W przypadku najprostszych operacji arytmetycznych można zastosować następujące reguły:
z = x + y -> z_err = u*abs(z) + x_err + y_errz = x - y -> z_err = u*abs(z) + x_err + y_errz = x * y -> z_err = u*abs(z) + x_err*abs(y) + y_err*abs(x)z = x / y -> z_err = u*abs(z) + (x_err*abs(y) + y_err*abs(x))/y^2z = sqrt(x) -> z_err = u*abs(z) + x_err/(2*abs(z))
gdzie u = eps/2jest zaokrąglenie jednostki. Tak, zasady +i -są takie same. Reguły dla każdej innej operacji op(x)można łatwo wyodrębnić za pomocą rozszerzenia serii Taylora zastosowanego wyniku op(x + x_err). Lub możesz spróbować google. Lub korzystając z książki Nicka Highama.
Jako przykład rozważmy następujący kod Matlab / Octave, który ocenia wielomiany we współczynnikach aw punkcie xprzy użyciu schematu Hornera:
function s = horner ( a , x )
s = a(end);
for k=length(a)-1:-1:1
s = a(k) + x*s;
end
W pierwszym etapie podzieliliśmy dwie operacje na s = a(k) + x*s:
function s = horner ( a , x )
s = a(end);
for k=length(a)-1:-1:1
z = x*s;
s = a(k) + z;
end
Następnie wprowadzamy _errzmienne. Zauważ, że dane wejściowe ai xzakłada się, że są dokładne, ale równie dobrze możemy również wymagać od użytkownika przekazania odpowiednich wartości dla a_erri x_err:
function [ s , s_err ] = horner ( a , x )
s = a(end);
s_err = 0;
for k=length(a)-1:-1:1
z = x*s;
z_err = ...;
s = a(k) + z;
s_err = ...;
end
Na koniec stosujemy zasady opisane powyżej, aby uzyskać warunki błędów:
function [ s , s_err ] = horner ( a , x )
u = eps/2;
s = a(end);
s_err = 0;
for k=length(a)-1:-1:1
z = x*s;
z_err = u*abs(z) + s_err*abs(x);
s = a(k) + z;
s_err = u*abs(s) + z_err;
end
Zauważ, że ponieważ nie mamy a_errlub x_errnp. Przyjmuje się, że są zerowe, odpowiednie wyrażenia są po prostu ignorowane w wyrażeniach błędów.
Zrobione! Mamy teraz schemat Hornera, który zwraca wynik szacunku błędu zależnego od danych (uwaga: jest to górna granica błędu) obok wyniku.
Na marginesie, skoro używasz C ++, możesz rozważyć utworzenie własnej klasy dla wartości zmiennoprzecinkowych, która przenosi ten _errtermin i przeciąża wszystkie operacje arytmetyczne, aby zaktualizować te wartości, jak opisano powyżej. W przypadku dużych kodów może to być łatwiejsza, choć mniej obliczeniowa, trasa. Powiedziawszy to, możesz być w stanie znaleźć taką klasę online. Szybkie wyszukiwanie w Google dało mi ten link .
± ux ( 1 ± u )