2011-09-18 7 views

risposta

11

Questa è una domanda apparentemente complessa. Here è una bella panoramica di alcuni possibili approcci.

+2

Potrei sembrare maleducato, ma penso che nel caso in cui il link non funzioni, dovresti aggiungere * alcuni * contenuti alla tua risposta. (Copia e incolla se vuoi) – Astrobleme

6

In considerazione del "link rot" che ha superato la risposta accettata, darò una risposta più autonoma concentrandosi sull'argomento di ottenere rapidamente un'ipotesi iniziale adatta all'iterazione superlineare.

Il "rilevamento" di metamerist (Wayback link) ha fornito alcuni confronti temporali per varie combinazioni di valore iniziale/iterazione (sono inclusi entrambi i metodi di Newton e Halley). I suoi riferimenti sono alle opere di W. Kahan, "Calcolo di una vera cubatura di root" e da K. Turkowski, "Calcolo del cubo di root".

metamarist aggiorna la tecnica bit-vaning dell'era DEC-VAX di W. Kahan con questo snippet, che "assume numeri interi a 32 bit" e si basa sul formato IEEE 754 per i doppi "per generare stime iniziali con 5 bit di precisione" :

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

Il codice di K. Turkowski fornisce leggermente maggiore precisione ("circa 6 bit") da un ridimensionamento convenzionale poteri-di-due su float fr, seguita da un'approssimazione quadratica alla sua radice cubica sull'intervallo [0.125 , 1.0):

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

e un successivo restauro dell'esponente di due (adattato a un terzo). L'estrazione esponente/mantissa e il restauro fanno uso di a frexp e ldexp.

Confronto con altri radice cubica "seme" approssimazioni

Per apprezzare quelle approssimazioni radice cubica di cui abbiamo bisogno per confrontarli con altre forme possibili. Prima i criteri per giudicare: consideriamo l'approssimazione nell'intervallo [1/8,1], e utilizziamo l'errore relativo migliore (minimizzando il massimo).

Cioè, se f(x) è una proposta di approssimazione alla x^{1/3}, troviamo il suo errore relativo:

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

L'approssimazione più semplice sarebbe ovviamente quella di utilizzare una sola costante nell'intervallo, e il miglior errore relativo in tal caso si ottiene selezionando f_0(x) = sqrt(2)/2, la media geometrica dei valori ai punti finali. Questo dà 1,27 bit di precisione relativa, un punto di partenza rapido ma sporco per un'iterazione di Newton.

Un'approssimazione migliore sarebbe la migliore polinomio di primo grado:

f_1(x) = 0.6042181313*x + 0.4531635984 

Questo dà 4.12 bit di precisione relativa, un grande miglioramento ma breve dei 5-6 bit di precisione relativa promesse dai rispettivi metodi di Kahan e Turkowski. Ma è nel ballpark e usa solo una moltiplicazione (e una aggiunta).

Infine, cosa succede se ci concediamo una divisione anziché una moltiplicazione?Risulta che con una divisione e due "aggiunte" possiamo sfruttare al meglio la funzione lineare frazionaria:

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

che dà 7.265 bit di precisione relativa.

A prima vista, questo sembra un approccio attraente, ma una vecchia regola empirica consisteva nel trattare il costo di una divisione FP come tre moltiplicazioni FP (e soprattutto ignorare le aggiunte e le sottrazioni). Tuttavia con gli attuali progetti FPU questo non è realistico. Mentre il costo relativo delle moltiplicazioni da aggiungere/sottrarre è sceso, nella maggior parte dei casi a un fattore pari a due o addirittura a parità, il costo della divisione non è diminuito, ma spesso è salito a 7-10 volte il costo della moltiplicazione. Quindi dobbiamo essere avari con le nostre operazioni di divisione.

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

Sembra che la questione di ottimizzazione è già stato affrontato, ma vorrei aggiungere un miglioramento alla funzione cubeRoot() postato qui, per le altre persone inciampo in questa pagina in cerca di un algoritmo di radice cubica rapida .

L'algoritmo esistente funziona bene, ma al di fuori dell'intervallo 0-100 fornisce risultati errati.

Ecco una versione modificata che funziona con numeri compresi tra -/+ 1 quadrilione (1E15). Se hai bisogno di lavorare con numeri più grandi, usa solo più iterazioni.

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

Per quanto riguarda l'ottimizzazione, sto cercando di indovinare il manifesto originale è stato chiesto come prevedere il numero minimo di iterazioni per un risultato preciso, data una dimensione di ingresso arbitrario. Ma sembra che per i casi più generali il guadagno ottenuto dall'ottimizzazione non valga la complessità aggiunta. Anche con la funzione sopra riportata, 100 iterazioni richiedono meno di 0,2 ms sull'hardware di consumo medio. Se la velocità fosse della massima importanza, prenderei in considerazione l'utilizzo di tabelle di ricerca pre-calcolate. Ma questo proviene da uno sviluppatore desktop, non da un ingegnere di sistemi embedded.

Problemi correlati