Calcolo in virgola mobile Powers (PHP/BCMath)

Question

Sto scrivendo un wrapper per l'estensione bcmath, e bug #10116 riguardo bcpow() è particolarmente fastidioso - proietta l'$right_operand ($exp) ad un (PHP nativo, non lunghezza arbitraria) intero, quindi quando si tenta di calcolare la radice quadrata (o qualsiasi altra radice superiore a 1) di un numero si finisce sempre con 1 invece del risultato corretto.

ho iniziato la ricerca di algoritmi che mi permettesse di calcolare la radice ennesima di un numero e mi found this answer che sembra piuttosto solido, in realtà ho expanded the formula usando WolframAlpha e sono stato in grado di migliorare la sua velocità di circa il 5%, mantenendo la precisione dei risultati.

Ecco un'implementazione PHP puro mimando mia implementazione BCMath e suoi limiti:

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; 
} 

Quanto sopra seems to work greattranne quando 1/fmod($exp, 1) non cedere un numero intero. Ad esempio, se $exp è 0.123456, il suo inverso sarà 8.10005 e l'esito della pow() e _pow() sarà un po 'diverso (demo):

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

Come posso ottenere lo stesso livello di accuratezza usando i calcoli esponenziali "manuali"?

Answer 1

L'algoritmo impiegato per trovare il n ^esimo radice di una (positiva) il numero a è l'algoritmo di Newton per trovare lo zero di

f(x) = x^n - a. 

che coinvolge solo i poteri con i numeri naturali come esponenti, quindi è semplice da implementare.

Calcolo potenza con un esponente 0 < y < 1 dove y non è nella forma 1/n con un numero intero n è più complicato. Facendo l'analogo, risolvendo

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

sarebbe di nuovo coinvolgere calcolare una potenza con esponente non integrale, il vero problema che stiamo cercando di risolvere.

Se y = n/d è razionale con piccolo denominatore d, il problema è facilmente risolvibile calcolando

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

ma per la maggior parte razionale 0 < y < 1, numeratore e denominatore sono piuttosto grandi, e l'intermedio x^n sarebbe enorme, in modo che il il calcolo userebbe molta memoria e impiegherebbe un tempo (relativamente) lungo. (Per l'esempio esponente del 0.123456 = 1929/15625, non è troppo male, ma 0.1234567 sarebbe piuttosto faticoso.)

Un modo per calcolare la potenza per il generale razionale 0 < y < 1 è quello di scrivere

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

con numeri interi a < b < c < ... < q e moltiplicare/dividere l'individuo x^(1/k). (Ogni razionale 0 < y < 1 ha tali rappresentazioni, e più brevi tali rappresentazioni in genere non comportano molti termini, ad esempio,

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

utilizzando solo le aggiunte nella decomposizione porta a rappresentazioni più lunghi con denominatori più grandi, ad esempio

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

in modo da richiedere più lavoro.)

Alcuni miglioramenti sono possibili mescolando gli approcci, trovare prima un'approssimazione razionale vicina a y con denominatore piccolo v l'espansione della frazione continua di - per l'esponente di esempio 1929/15625 = [0;8,9,1,192] e l'utilizzo dei primi quattro quozienti parziali restituisce l'approssimazione 10/81 = 0.123456790123... [nota che 10/81 = 1/8 - 1/648, le somme parziali della scomposizione più breve in frazioni pure sono convergenti] - e quindi scomporre il resto in puro frazioni.

Tuttavia, in generale, tale approccio porta a calcolare n ^th radici per grandi n, che è anche lenta e molta memoria se la precisione desiderata del risultato finale è alta.

Tutto sommato, è probabilmente più semplice e veloce da implementare e explog e utilizzare

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