Funzione ERF in C++

[KRYHAWOK]

Ciao gente.

Qualcuno mi sa dire dove posso recuperare una libreria C++ che contenga la funzione ERF, ossia l'integrale della gaussiana?

Grazie!

[KRYHAWOK]

Dopo un po' di mosse penso di esserci arrivato:
il problema sta nel fatto che la funzione gaussiana non ha primitiva, quindi il calcolo dell'area sottesa è solo numerico: erf(x) non ha espressione analitica.

La gaussiana standardizzata è un'esponenziale, quindi trasformabile in una serie uniformemente convergente su tutto R. Ho poi integrato per serie per approssimare la funzione a mio piacimento. Per avere una precisione fino a 10^(-4) è necessario sommare non meno di 15 termini della serie ottenuta, ma devo ammettere che non è poi tanto male.

Ecco quindi l'espressione approssimata dell'area sottesa dalla gaussiana standardizzata in funzione di x (con f[0]=0). A partire da questa basta un bel ciclo FOR e il calcolo è fatto.
A chi può servire:

http://utenti.lycos.it/plmnko/Normale.jpg

Se qualcuno ha un'idea migliore mi faccia sapere!

[KRYHAWOK]

E' sorto un altro problema; pensavo che una volta trovata una formula analitica per il calcolo il gioco fosse finito.

Ho provato a tradurre in C++ la funzione che ho riportato sopra nel modo più ovvio, ossia:

*****************************************
long double normale(float x) {
long double sum=0;
for(i=0; i<=15; i++)
sum=sum+(pow(2,-i)*pow(x,2*i+1)*pow(-1,i))/((2*i+1)*fatt(i));
return sum/sqrt(1/(2*3.1415926)); }
*****************************************

..che traduce la formula. Ma i valori ottenuti non sono soddisfacenti: la stessa equazione, in Derive, dà un risultato mooolto migliore.

Ho provato allora ad eliminare il ciclo FOR, espcitando la sommatoria e scrivendo esplicitamente i coefficienti. In questo modo i valori sono nettamente più simili a quelli di Derive, ma non oltre x=0.7. Considerando che la normale è praticamente costante dopo x=3.7.... non ci siamo ancora.

Come posso tradurre in C++ la formula dell'immagine in modo decente?????

Grazie dell'aiuto!

[/\/\@®¢Ø]

Se non ricordo male ci sono degli algoritmi appositi per l'integrazione numerica, per renderla sia piu' veloce che piu' precisa.
Ora non chiedermi quali :D, pero' se ripesco il libro di calcolo numerico magari posso esserti d'aiuto.

Hai provato a compilare senza le ottimizzazioni sui floating point ?

[KRYHAWOK]

Cosa intendi "senza le ottimizzazioni sui floating point"? Scusa ma sto perdendo lucidità....

[/\/\@®¢Ø]

Intendevo dire che tra le ottimizzazioni del compilatore, spesso ce ne sono alcune che aumentano la velocita' sacrificando un po' la precisione sui floating point. Potresti vedere se nel tuo caso sono stato attivate ed eventualmente toglierle
( che compilatore usi ? )

[KRYHAWOK]

All'inizio pensavo che fosse colpa di approssimazioni del compilatore... poi mi son messo a spezzare la formula, monitorando man mano come cambiava il valore ottenuto dal ciclo FOR. I termini che si aggiungono alla serie dovrebbero essere più piccoli man mano che si procede, se voglio che la serie converga.... ma invece ho visto che fluttuavano abbastanza casualmente....

In sostanza, l'errore non era dovuto nè al compilatore, nè all'algoritmo, ma semplicemente alla funzione fatt(x) che avevo importato da una vecchia libreria che avevo fatto..... in quella libreria avevo impostato la funzione fattoriale così:

int fatt(int x) { ............}

Ovviamente andava bene finchè i fattoriali da calcolare erano piccoli... ma nell'algoritmo compare anche 15!
bastava modificare la funzione in modo che restituisse valori non interi<65535 (già 9!>65535............)

Se dio vuole adesso la funzione restituisce i valori corretti.... basta stare attenti a ricordarsi che da x=3.9 circa in poi la funzione è identicamente = 0.5, mente l'algoritmo di sopra non è in grado di reggere oltre un certo intervallo perchè la successione non converge più a ERF(x).

Ecco il grafico della serie calcolata fino al quindicesimo termine in rosso, che si sovrappone perfettamente al grafico della ERF(x) vera e propria (in nero, visibile agli estremi).

http://utenti.lycos.it/plmnko/erf.jpg

Cmq grazie dell'aiuto!