Obliczenia zmiennoprzecinkowe Powers (PHP/BCMath)

Question

Piszę otoki dla rozszerzenia bcmath i bug #10116 dotyczące bcpow() jest szczególnie irytujące - to rzuca $right_operand ($exp) do (native PHP, a nie dowolnej długości) integer, więc gdy spróbujesz obliczyć pierwiastek kwadratowy (lub dowolny inny root wyższy niż 1) liczby, zawsze otrzymasz numer 1 zamiast poprawnego wyniku.

zacząłem poszukiwania algorytmów, które pozwoliłyby mi obliczyć pierwiastkowanie szeregu i ja found this answer który wygląda całkiem solidne, tak naprawdę expanded the formula użyciu WolframAlpha i byłem w stanie poprawić jego prędkość o około 5%, przy jednoczesnym zachowaniu dokładności wyników.

Oto czysta realizacja PHP naśladując moje wykonanie bcmath i jego ograniczeń:

function _pow($n, $exp) 
{ 
    $result = pow($n, intval($exp)); // bcmath casts $exp to (int) 

    if (fmod($exp, 1) > 0) // does $exp have a fracional part higher than 0? 
    { 
     $exp = 1/fmod($exp, 1); // convert the modulo into a root (2.5 -> 1/0.5 = 2) 

     $x = 1; 
     $y = (($n * _pow($x, 1 - $exp))/$exp) - ($x/$exp) + $x; 

     do 
     { 
      $x = $y; 
      $y = (($n * _pow($x, 1 - $exp))/$exp) - ($x/$exp) + $x; 
     } while ($x > $y); 

     return $result * $x; // 4^2.5 = 4^2 * 4^0.5 = 16 * 2 = 32 
    } 

    return $result; 
} 

Powyższy seems to work greatwyjątkiem gdy 1/fmod($exp, 1) nie daje liczbę całkowitą. Na przykład, jeśli $exp jest 0.123456 jej odwrotne będą 8.10005 i wyniki pow() i _pow() będzie nieco różne (demo):

• pow(2, 0.123456) = 1.0893412745953
• _pow(2, 0.123456) = 1.0905077326653
• _pow(2, 1/8) = _pow(2, 0.125) = 1.0905077326653

Jak osiągnąć ten sam poziom dokładności przy użyciu "ręcznych" obliczeń wykładniczych?

Answer 1

Zastosowany algorytm, aby znaleźć n ^th korzeń (pozytywny) Numer a jest algorytm Newton znajdowania zero

f(x) = x^n - a. 

, która obejmuje tylko uprawnienia z liczb naturalnych jak propagatorów, a tym samym jest proste do wdrożenia.

Obliczanie mocy wykładnikiem 0 < y < 1 gdzie y nie jest w postaci 1/n z całkowitą n jest bardziej skomplikowane. Robi analog, rozwiązywanie

x^(1/y) - a == 0 

znowu wiązać obliczania mocy z niezintegrowanych wykładnik samego problemu, my próbujesz rozwiązać.

Jeśli y = n/d jest racjonalne z małym mianownika d, problem można łatwo rozwiązać poprzez obliczenie

x^(n/d) = (x^n)^(1/d), 

ale dla najbardziej racjonalnym 0 < y < 1, licznik i mianownik są dość duże, a pośredni x^n byłby ogromny, więc obliczenia zużywają dużo pamięci i zajmują (stosunkowo) długi czas. (dla przykładu wykładnik 0.123456 = 1929/15625, to nie jest tak źle, ale 0.1234567 byłoby raczej opodatkowanie).

Jednym ze sposobów, aby obliczyć moc do ogólnego racjonalnego 0 < y < 1 jest napisać

y = 1/a ± 1/b ± 1/c ± ... ± 1/q 

z całkowitymi a < b < c < ... < q i do pomnożenia/podzielenia indywidualnego x^(1/k). (Każdy racjonalny 0 < y < 1 ma takie reprezentacje, a najkrótsze takie reprezentacje generalnie nie wymagają wiele terminów, np

1929/15625 = 1/8 - 1/648 - 1/1265625; 

używając tylko dodatki w dekompozycja prowadzi do dłuższych przedstawień z większymi mianowników, np

1929/15625 = 1/9 + 1/82 + 1/6678 + 1/46501020 + 1/2210396922562500, 

, co wymagałoby więcej pracy).

Pewne ulepszenie jest możliwe dzięki połączeniu podejść, najpierw znajdź bliskie racjonalne przybliżenie do y z małym mianownikiem v np. kontynuacja ekspansji frakcji y - dla przykładu wykładnik 1929/15625 = [0;8,9,1,192] i przy użyciu pierwszych czterech częściowych ilorazów uzyskuje się aproksymację 10/81 = 0.123456790123... [zauważ, że 10/81 = 1/8 - 1/648, częściowe sumy najkrótszego rozkładu w czyste frakcje są konwergentami] - a następnie rozkładają pozostałość na czysty frakcje.

Ogólnie rzecz biorąc to podejście prowadzi do obliczania korzeni dla dużych , które również jest wolne i wymaga dużej ilości pamięci, jeśli pożądana dokładność końcowego wyniku jest wysoka.

W sumie, to chyba prostsze i szybsze do wdrożenia exp i log i używać

x^y = exp(y*log(x))