Come vengono calcolati gli esponenti?
Domanda
Sto cercando di determinare il runtime asintotico di uno dei miei algoritmi, che utilizza esponenti, ma non sono sicuro di come gli esponenti vengano calcolati a livello di programmazione.
Sto specificatamente cercando l'algoritmo pow () usato per i numeri a virgola mobile e precisione doppia.
Soluzione
Ho avuto la possibilità di esaminare l'implementazione di fdlibm. I commenti descrivono l'algoritmo utilizzato:
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 53-24 = 29 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'|<=0.5.
* 3. Return x**y = 2**n*exp(y'*log2)
seguito da un elenco di tutti i casi speciali gestiti (0, 1, inf, nan).
Le sezioni più intense del codice, dopo tutta la gestione dei casi speciali, riguardano i calcoli log2
e 2 **
. E non ci sono anelli in nessuno di questi. Quindi, nonostante la complessità delle primitive a virgola mobile, sembra un algoritmo asintoticamente a tempo costante.
Gli esperti in virgola mobile (di cui non sono uno) sono invitati a commentare. : -)
Altri suggerimenti
A meno che non abbiano scoperto un modo migliore per farlo, credo che i valori approssimativi di trig, funzioni logaritmiche ed esponenziali (ad esempio per crescita esponenziale e decadimento) siano generalmente calcolati usando le regole aritmetiche e Taylor Series espansioni per produrre un risultato approssimativo accurato con la precisione richiesta. (Vedi qualsiasi libro di calcolo per i dettagli sulle espansioni di funzioni delle serie di potenza, serie Taylor e serie Maclaurin.) Si noti che è passato un po 'di tempo da quando ho fatto tutto questo, quindi non potrei dirti, ad esempio, esattamente come calcolare il numero di termini della serie che è necessario includere garantiscono un errore sufficientemente piccolo da essere trascurabile in un calcolo a doppia precisione.
Ad esempio, l'espansione della serie Taylor / Maclaurin per e ^ x è questa:
+inf [ x^k ] x^2 x^3 x^4 x^5
e^x = SUM [ --- ] = 1 + x + --- + ----- + ------- + --------- + ....
k=0 [ k! ] 2*1 3*2*1 4*3*2*1 5*4*3*2*1
Se prendi tutti i termini (k da 0 a infinito), questa espansione è esatta e completa (nessun errore).
Tuttavia, se non prendi tutti i termini all'infinito, ma ti fermi dopo aver detto 5 termini o 50 termini o qualsiasi altra cosa, produci un risultato approssimativo che differisce dal reale e ^ x valore della funzione per un resto che è abbastanza facile da calcolare.
La buona notizia per gli esponenziali è che converge bene e i termini della sua espansione polinomiale sono abbastanza facili da codificare in modo iterativo, quindi potresti (ripetere, POTERE - ricorda , è passato del tempo), non è nemmeno necessario pre-calcolare quanti termini è necessario per garantire che il tuo errore sia inferiore alla precisione perché puoi testare la dimensione del contributo ad ogni iterazione e fermarti quando diventa abbastanza vicino a zero. In pratica, non so se questa strategia sia praticabile o meno - dovrei provarla. Ci sono dettagli importanti che ho da tempo dimenticato. Roba come: precisione della macchina, errore della macchina ed errore di arrotondamento, ecc.
Inoltre, tieni presente che se non stai usando e ^ x, ma stai facendo crescita / decadimento con un'altra base come 2 ^ x o 10 ^ x, la funzione polinomiale approssimativa cambia.
L'approccio abituale, per elevare a a b, per un esponente intero, va in questo modo:
result = 1
while b > 0
if b is odd
result *= a
b -= 1
b /= 2
a = a * a
È generalmente logaritmico nella dimensione dell'esponente. L'algoritmo si basa sull'invariante "a ^ b * risultato = a0 ^ b0", dove a0 e b0 sono i valori iniziali di a e b.
Per esponenti negativi o non interi, sono necessari logaritmi e approssimazioni e analisi numerica. Il tempo di esecuzione dipenderà dall'algoritmo utilizzato e dalla precisione per cui è sintonizzata la libreria.
Modifica: poiché sembra esserci un certo interesse, ecco una versione senza moltiplicazioni extra.
result = 1
while b > 0
while b is even
a = a * a
b = b / 2
result = result * a
b = b - 1
Se scrivessi una funzione pow destinata a Intel, restituirei exp2 (log2 (x) * y). Il microcodice Intel per log2 è sicuramente più veloce di qualsiasi altra cosa che sarei in grado di codificare, anche se potessi ricordare il calcolo del primo anno e l'analisi numerica della scuola elementare.
Puoi usare exp (n * ln (x)) per calcolare x n . Sia x che n possono essere numeri in virgola mobile a precisione doppia. Il logaritmo naturale e la funzione esponenziale possono essere calcolati usando la serie di Taylor. Qui puoi trovare le formule: http://en.wikipedia.org/wiki/Taylor_series