Siewanie iteracji Newtona dla korzenia kostki efektywnie

Question

Jak mogę efektywnie znaleźć pierwiastek kostki z liczby? Myślę, że metoda Newtona-Raphsona może być używana, ale nie wiem, jak odgadnąć początkowe rozwiązanie programowo, aby zminimalizować liczbę iteracji.

Answer 1

Jest to złudnie skomplikowane pytanie. Here to miła ankieta dotycząca niektórych możliwych podejść.

Answer 2

W związku z "zgnilizną łącza", która wyprzedziła Akceptowaną odpowiedź, udzielę bardziej samodzielnej odpowiedzi, koncentrując się na temacie szybkiego uzyskania wstępnego przypuszczenia odpowiedniego dla itelacji superliniowej.

"Badanie" wykonane przez metamerystę (Wayback link) dostarczyło pewnych porównań czasowych dla różnych kombinacji wartości początkowych/iteracji (uwzględniono obie metody Newtona i Halleya). Jego referencje dotyczą prac W. Kahan, "Computing a Real Cube Root" oraz K. Turkowski, "Computing the Cube Root".

Metamorysta aktualizuje technikę bit-fiddlingową w stylu DEC-VAX W. Kahana za pomocą tego fragmentu, który "zakłada 32-bitowe liczby całkowite" i opiera się na formacie IEEE 754 dla debelów "w celu wygenerowania wstępnych oszacowań z 5 bitami precyzji" :

inline double cbrt_5d(double d) 
{ 
    const unsigned int B1 = 715094163; 
    double t = 0.0; 
    unsigned int* pt = (unsigned int*) &t; 
    unsigned int* px = (unsigned int*) &d; 
    pt[1]=px[1]/3+B1; 
    return t; 
} 

Kod K. Turkowski zapewnia nieznacznie większą precyzję („około 6 bitów”) przez konwencjonalne potęgi dwójki skalowaniem float fr, a następnie kwadratowej zbliżania do jego pierwiastek na przedziale [0,125 , 1.0):

/* Compute seed with a quadratic qpproximation */ 
fr = (-0.46946116F * fr + 1.072302F) * fr + 0.3812513F;/* 0.5<=fr<1 */ 

, a następnie przywrócenie wykładnika potęgi dwóch (dostosowanej do jednej trzeciej). Eksploatacja i odtworzenie wykładnika/mantysy używają math library calls do frexp i ldexp.

Porównanie z innymi pierwiastek „nasienie” przybliżeń

Aby docenić tych przybliżeń pierwiastek musimy porównać je z innych możliwych form. Najpierw kryteria oceny: uważamy przybliżenie przedziału [1/8,1] i używamy najlepszego (minimalizującego maksimum) błędu względnego.

Oznacza to, że jeśli f(x) jest proponowana przybliżeniem x^{1/3}, znaleźć swój względny błąd:

 error_rel = max | f(x)/x^(1/3) - 1 | on [1/8,1] 

Najprostszym przybliżeniem byłoby oczywiście użyć pojedynczego stałą w przedziale, a najlepiej względny błąd w takim przypadku osiąga się przez wybranie f_0(x) = sqrt(2)/2, średniej geometrycznej wartości w punktach końcowych. Daje to 1,27 bitów względnej dokładności, szybki, ale brudny punkt początkowy dla iteracji Newtona.

lepsze przybliżenie byłby najlepszy wielomian pierwszego stopnia:

f_1(x) = 0.6042181313*x + 0.4531635984 

Daje 4.12 bitów względną dokładność, duża poprawa, ale brakuje 5-6 bitów względnej dokładności obiecanych przez odpowiednich metod Kahana i Turkowskiego. Ale jest w parku i używa tylko jednego mnożenia (i jednego dodatku).

Wreszcie, jeśli pozwolimy sobie na podział zamiast na mnożenie?Okazuje się, że z jednej dywizji i dwóch „dodatkami”, możemy mieć najlepszy liniowy-cząstkowej funkcji:

f_M(x) = 1.4774329094 - 0.8414323527/(x+0.7387320679) 

co daje 7.265 bitów względnej dokładności.

Na pierwszy rzut oka wydaje się to atrakcyjnym podejściem, ale stara zasada polegała na traktowaniu kosztu dywizji FP jak trzy multiplikacje FP (i na ogół ignorowanie dodatków i odejmowań). Jednak przy obecnych projektach FPU jest to nierealne. Podczas gdy względny koszt mnożenia do dodawania/odejmowania spadł, w większości przypadków do współczynnika równego dwa lub nawet równości, koszt podziału nie spadł, ale często wzrósł do 7-10 razy więcej niż koszt mnożenia. Dlatego musimy być skąpi wobec naszych operacji podziału.

Answer 3

static double cubeRoot(double num) { 
    double x = num; 

    if(num >= 0) { 
     for(int i = 0; i < 10 ; i++) { 
      x = ((2 * x * x * x) + num)/(3 * x * x); 
     } 
    } 
    return x; 
} 

Answer 4

Wydaje się, że w kwestii optymalizacji zostało już omówione, ale chciałbym dodać do poprawy funkcji Cuberoot() Zamieszczone tu dla innych ludzi potykając się na tej stronie szukających korzeni szybki algorytm kostki .

Istniejący algorytm działa dobrze, ale poza zakresem od 0 do 100 daje nieprawidłowe wyniki.

Oto poprawiona wersja, która działa z liczbami między -/+ 1 biliardem (1E15). Jeśli potrzebujesz pracować z większymi liczbami, po prostu użyj więcej iteracji.

static double cubeRoot(double num){ 
    boolean neg = (num < 0); 
    double x = Math.abs(num); 
    for(int i = 0, iterations = 60; i < iterations; i++){ 
     x = ((2 * x * x * x) + num)/(3 * x * x); 
    } 
    if(neg){ return 0 - x; } 
    return x; 
} 

Odnośnie optymalizacji Zgaduję oryginalny plakat był pytaniem, jak przewidują minimalną liczbę iteracji dla dokładnego wyniku, biorąc pod uwagę dowolny rozmiar wejściowy. Wydaje się jednak, że w przypadku większości ogólnych korzyści z optymalizacji nie są warte dodatkowej złożoności. Nawet przy powyższej funkcji 100 iteracji zajmuje mniej niż 0,2 ms przeciętnego sprzętu komputerowego. Gdyby szybkość była niezwykle ważna, rozważałabym użycie tabel obliczeń obliczanych wstępnie. Ale to pochodzi od programisty komputerowego, a nie od inżyniera systemów wbudowanych.