molto lento std :: pow() per le basi molto vicino a 1

Question

Ho un codice numerico che risolve un'equazione f(x) = 0, in cui devo sollevare x a una potenza p. Lo risolvo usando un sacco di cose, ma alla fine ho il metodo di Newton. La soluzione sembra essere uguale a x = 1, e quindi è la causa dei miei problemi. Quando la soluzione iterata si avvicina a 1, ad esempio x = 1 + 1e-13, il tempo necessario per calcolare lo std::pow(x, p) cresce enormemente, facilmente di un fattore pari a 100, rendendo il mio codice inutilizzabile.

La macchina che esegue questa operazione è AMD64 (Opteron 6172) su CentOS e il comando è semplice . Un comportamento simile si presenta su tutte le mie macchine, tutte x64. Come documentato here, questo non è solo il mio problema (ad esempio, qualcun altro è incazzato troppo), appare solo su x64 e solo per vicino a 1.0. La cosa simile accade a exp.

Risolvere questo problema è vitale per me. Qualcuno sa se c'è qualche modo per aggirare questa lentezza?

EDIT: John ha sottolineato che ciò è dovuto ai denormali. La domanda è quindi, come risolvere questo problema? Il codice è C++, compilato con g++ per l'uso all'interno di GNU Octave. Sembra che, sebbene io abbia impostato CXXFLAGS per includere -mtune=native e -ffast-math, ciò non sta aiutando e il codice viene eseguito altrettanto lentamente.

SOLUZIONE PSEUDO PER ORA: a tutti coloro che si preoccupano di questo problema, le soluzioni suggerite di seguito non hanno funzionato per me personalmente. Ho davvero bisogno della solita velocità di std::pow(), ma senza la lentezza intorno a x = 1. La soluzione per me è personalmente ricorrendo a un trucco:

inline double mpow(double x, double p) __attribute__ ((const)); 

inline double mpow(double x, double p) 
{ 
    double y(x - 1.0); 
    return (std::abs(y) > 1e-4) ? (std::pow(x, p)) : (1.0 + p * y * (1.0 + (p - 1.0) * y * (0.5 + (1.0/6.0) * (p - 2.0) * y))); 
} 

il limite potrebbe essere cambiato, ma per -40 < p < 40 l'errore è più piccolo di circa 1e-11, che è abbastanza buono. L'overhead è minimo da ciò che ho trovato, quindi questo risolve il problema per me.

Answer 1

Potrebbe essere anche il tuo algoritmo. Forse passare a qualcosa come BFGS invece del metodo di Newton aiuterà.

Non si dice nulla sui criteri di convergenza. Forse anche quelli hanno bisogno di essere aggiustati.

Answer 2

L'ovvia soluzione è quella di notare che nei valori reali, a ** b == exp(log(a) * b), si utilizza invece tale modulo. Dovrai controllare che non influenzi negativamente la precisione dei risultati. Modifica: come discusso, anche questo soffre del rallentamento a un grado quasi altrettanto elevato.

Il problema non sono denormali, almeno non direttamente; il tentativo di calcolare exp(-2.4980018054066093e-15) soffre lo stesso rallentamento e -2.4980018054066093e-15 non è certamente un diniego.

Se non vi interessa circa l'accuratezza dei risultati, quindi il ridimensionamento o l'exponend o l'esponente dovrebbe farti fuori della zona lento:

sqrt(pow(a, b * 2)) 
pow(a * 2, b)/pow(2, b) 
... 

Questo bug è noto per manutentori glibc: http://sourceware.org/bugzilla/show_bug.cgi?id=13932 - se stai cercando una soluzione alternativa a una soluzione alternativa, ti consigliamo di commissionare a un esperto di matematica in virgola mobile un'esperienza open source.

Answer 3

64-bit Linux?

Utilizzare il codice pow() da FreeBSD.

La libreria di Linux C (glibc) ha prestazioni orribili nel caso peggiore per determinati input.

Vedi: http://entropymine.com/imageworsener/slowpow/