Funzione ipergeometrica di Gauss in C
Re: Funzione ipergeometrica di Gauss in C
Vedi allegato.
- Allegati
-
- hypergeometric.tgz
- (8.48 KiB) Scaricato 40 volte
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Di nuovo grazie, ottimo lavoro.
Ho fatto un po' di prove confrontando con i risultati di wolphram alpha: i risultati sono coincidenti, variando al più l'ultima cifra decimale.
Due osservazioni:
1, di carattere matematico: mi aspetterei che se a, b, c, z sono tutti reali, la parte immaginaria della funzione deve essere nulla, ma non è così nemmeno per wolphram (per esempio, plot per a=1, b=2, c=4, z tra 0 e 10)
2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Ho fatto un po' di prove confrontando con i risultati di wolphram alpha: i risultati sono coincidenti, variando al più l'ultima cifra decimale.
Due osservazioni:
1, di carattere matematico: mi aspetterei che se a, b, c, z sono tutti reali, la parte immaginaria della funzione deve essere nulla, ma non è così nemmeno per wolphram (per esempio, plot per a=1, b=2, c=4, z tra 0 e 10)
2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Re: Funzione ipergeometrica di Gauss in C
La funzione ipergeometrica è definita come la soluzione dell'equazione differenziale ipergeometrica di Eulero, e non è stupefacente che un'equazione differenziale del secondo ordine a coefficienti reali abbia una soluzione esprimibile in forma complessa: pensa ad esempio all'oscillatore armonico. La serie ipergeometrica risolve l'equazione ipergeometrica ma ha raggio di convergenza 1 (non 0.5 come erroneamente avevo scritto). Per z reale compreso tra 0 e 1 il valore di 2F1 dovrebbe sempre essere reale.giuseppe morello ha scritto:Di nuovo grazie, ottimo lavoro.
Ho fatto un po' di prove confrontando con i risultati di wolphram alpha: i risultati sono coincidenti, variando al più l'ultima cifra decimale.
Due osservazioni:
1, di carattere matematico: mi aspetterei che se a, b, c, z sono tutti reali, la parte immaginaria della funzione deve essere nulla, ma non è così nemmeno per wolphram (per esempio, plot per a=1, b=2, c=4, z tra 0 e 10)
Vero, ma non so come rimediare. hmin viene passato come argomento alla odeint, ma cambiarlo non risolve.
2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Ho visto che la GNU scientific library ha la funzione ipergeometrica (vedi http://www.gnu.org/software/gsl/manual/ ... tions.html); ti conviene provarla.
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Scusa l'imbranataggine, ma non sono ancora molto pratico. Ho provato a usare La GNU Scientific Library, anche se non so se alla fine c'è tutto quello che mi serve, comunque l'ho installata sul mio portatile e ho provato a fare girare il semplice codice ipergeom.c:
Compilando con apparentemente si produce l'eseguibile, ma provando ad eseguirlo mi appare "Permission denied".
Preciso che nel portatile non ho Ubuntu, ma Mac OS X, però ho provato anche su un sistema Ubuntu, con medesimo risultato.
Ho letto che forse dovrei specificare il path degli header usati, ho provato a cercarli con; nel Mac non riconosce il comando, mentre su Ubuntu non trova il pacchetto.
Probabilmente mi sto perdendo in un bicchiere d'acqua, ma non conosco bene la logica di questi sistemi.
Codice: Seleziona tutto
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<gsl/gsl_sf_hyperg.h>
double a, b, c, x, iper;
int main()
{
printf("Inserisci i valori a, b, c, x:\n");
scanf("%lf %lf %lf %lf", &a, &b, &c, &x);
iper=hyperg_2F1_series(a,b,c,x);
printf("%18.15f\n", iper);
return 0;
}
Codice: Seleziona tutto
gcc -c -lm ipergeom.c -o ipergeom.out
Preciso che nel portatile non ho Ubuntu, ma Mac OS X, però ho provato anche su un sistema Ubuntu, con medesimo risultato.
Ho letto che forse dovrei specificare il path degli header usati, ho provato a cercarli con
Codice: Seleziona tutto
pkg-config --libs gsl_sf_hyperg.h
Probabilmente mi sto perdendo in un bicchiere d'acqua, ma non conosco bene la logica di questi sistemi.
Re: Funzione ipergeometrica di Gauss in C
Se si produce un eseguibile significa che gli header sono stati trovati, e anche le librerie.
Mi sembra strano, però, che abbia compilato con il comando che hai riportato. Io ho dovuto dare:
Da dove ti salta fuori hyperg_2F1_series? I nomi di tutte le funzioni della libreria gsl cominciano con gsl_
Io ho usato gsl_sf_hyperg_2F1 (che però converge solo per |x|<1).
... e ahimé, a ben guardare il manuale _tutte_ le funzioni ipergeometriche della gsl hanno questa limitazione.
Comunque casi come quello che hai citato nel tuo messaggio precedente, con a o b interi negativi, li puoi gestire facilmente come casi speciale, perché i simboli di Pochhammer per interi negativi diventano 0 dopo un certo indice e quindi la serie ipergeometrica risulta avere un numero finito di termini.
Mi sembra strano, però, che abbia compilato con il comando che hai riportato. Io ho dovuto dare:
Codice: Seleziona tutto
gcc -lm -lgsl -lgslcblas -o eseguibile sorgente.c
Io ho usato gsl_sf_hyperg_2F1 (che però converge solo per |x|<1).
... e ahimé, a ben guardare il manuale _tutte_ le funzioni ipergeometriche della gsl hanno questa limitazione.
Comunque casi come quello che hai citato nel tuo messaggio precedente, con a o b interi negativi, li puoi gestire facilmente come casi speciale, perché i simboli di Pochhammer per interi negativi diventano 0 dopo un certo indice e quindi la serie ipergeometrica risulta avere un numero finito di termini.
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Non so cosa facesse, ovviamente il con il tuo comando invece funziona. Devo assolutamente studiare un po' o seguire qualche corso, perché a cercare come si compila sul web spesso mi confondo.bite ha scritto:Se si produce un eseguibile significa che gli header sono stati trovati, e anche le librerie.
Mi sembra strano, però, che abbia compilato con il comando che hai riportato. Io ho dovuto dare:
Codice: Seleziona tutto
gcc -lm -lgsl -lgslcblas -o eseguibile sorgente.c
Svista.bite ha scritto: Da dove ti salta fuori hyperg_2F1_series? I nomi di tutte le funzioni della libreria gsl cominciano con gsl_
Gli esempi che facevo erano del tutto casuali, non lavoro solo con interi. Potrebbe anche darsi che il codice che mi hai inviato mi basti (sulla base delle poche prove fatte). Per verificarlo devo prima terminare il programma, per il quale mi mancano ancora una funzione abbastanza facile da computare e in più l'ipergeometrica F1 di Appell, per la quale ancora non ho trovato moltissimo (non mi risulta né su GSL né su Numerical Recipes).bite ha scritto: Io ho usato gsl_sf_hyperg_2F1 (che però converge solo per |x|<1).
... e ahimé, a ben guardare il manuale _tutte_ le funzioni ipergeometriche della gsl hanno questa limitazione.
Comunque casi come quello che hai citato nel tuo messaggio precedente, con a o b interi negativi, li puoi gestire facilmente come casi speciale, perché i simboli di Pochhammer per interi negativi diventano 0 dopo un certo indice e quindi la serie ipergeometrica risulta avere un numero finito di termini.
Re: Funzione ipergeometrica di Gauss in C
http://mathworld.wolfram.com/AppellHype ... ction.html
arxiv.org/pdf/1105.3565
arxiv.org/pdf/hep-ph/0110083
Buon divertimento
Un giorno poi ci spiegherai che cosa stai combinando
arxiv.org/pdf/1105.3565
arxiv.org/pdf/hep-ph/0110083
Buon divertimento
Un giorno poi ci spiegherai che cosa stai combinando
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Lo sviluppo in serie e' abbastanza facile, l'ho fatto sulla falsariga di quello di 2F1. Ora devo sperare di riuscire ad adattare qualcosa tipo odeint
Quel che sto facendo non e' niente di eccezionale: semplicemente sto cercando di tradurre in C e adattare meglio ai miei scopi un codice che simula curve di luce di transiti, solo che l'originale non ricorre a queste funzioni.
Quel che sto facendo non e' niente di eccezionale: semplicemente sto cercando di tradurre in C e adattare meglio ai miei scopi un codice che simula curve di luce di transiti, solo che l'originale non ricorre a queste funzioni.
Re: Funzione ipergeometrica di Gauss in C
Temo sia impossibile. odeint sta per Ordinary DIfferential Equation INTegrator, mentre le equazioni differenziali costitutive delle funzioni di Appell sono alle derivate parziali.giuseppe morello ha scritto:Lo sviluppo in serie e' abbastanza facile, l'ho fatto sulla falsariga di quello di 2F1. Ora devo sperare di riuscire ad adattare qualcosa tipo odeint
Prima di fasciarti la testa, però, verifica qual è il raggio di convergenza della serie, può anche darsi che ti basti.
Il secondo paper che ti ho linkato parla di "higher order radiative corrections to scattering amplitudes"; per questa applicazione sembra che sia sufficiente uno sviluppo della F1 in un parametro piccolo attorno ad un valore noto. Se pensi che possa andar bene anche per te, dai uno sguardo ai paragrafi 1 e 4.2.Quel che sto facendo non e' niente di eccezionale: semplicemente sto cercando di tradurre in C e adattare meglio ai miei scopi un codice che simula curve di luce di transiti, solo che l'originale non ricorre a queste funzioni.
In alternativa, visto che la F1 ha anche una rappresentazione in termini di un integrale semplice (nel senso di non doppio ) può darsi che ti convenga seguire questa strada.
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Mi sa che questa e' l'idea giusta . Il dominio di convergenza della serie non era sufficiente, quello dell'integrale pare di si' invece.bite ha scritto: In alternativa, visto che la F1 ha anche una rappresentazione in termini di un integrale semplice (nel senso di non doppio ) può darsi che ti convenga seguire questa strada.
Al momento l'ho calcolato semplicemente con il metodo di Simpson, devo solo perfezionarlo per avere una convergenza piu' rapida. Ti aggiornero'.
Re: Funzione ipergeometrica di Gauss in C
Attenzione però che l'integrando può avere delle singolarità a seconda dei valori dei parametri.giuseppe morello ha scritto: Al momento l'ho calcolato semplicemente con il metodo di Simpson, devo solo perfezionarlo per avere una convergenza piu' rapida. Ti aggiornero'.
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Avevo notato le singolarità, speravo non creassero particolari problemi, invece li creano eccome, la convergenza è lentissima.
Posto il mio codice per il calcolo dell'integrale:
appell_integ.c (l'else commentato indica la prima strada, che è quella più naturale, con passo costante; poi invece ho optato per un passo variabile in modo da accorciarlo in prossimità delle possibili singolarità)
appellf1_integrand.c
beta.c
Valori dei parametri che mi possono interessare:
0.5, 1, 0.5, 2.75, -399, 0.12
0.5, -2, 1, 1, 0.12, -48
Posto il mio codice per il calcolo dell'integrale:
appell_integ.c (l'else commentato indica la prima strada, che è quella più naturale, con passo costante; poi invece ho optato per un passo variabile in modo da accorciarlo in prossimità delle possibili singolarità)
Codice: Seleziona tutto
#include<math.h>
//#include<complex.h>
#include<stdio.h>
#include"beta.h"
#include"appellf1_integrand.h"
int main()
{
double integ[30], a, b, c, d, x, y, t, f, h, lim;
int i, n;
//eps=pow(10.,-15.);
//printf("eps=%18.15lf\n", eps);
printf("Inserisci i valori di a, b, c, d, x, y:\n");
scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &x, &y);
if(a<=0.)
{
printf("Can't accept a<=0 in appell_integ\n");
}
else if(d-a<=0.)
{
//nerror("Can't accept d<=a in appell_integ\n");
printf("Can't accept d<=a in appell_integ\n");
}
else if(x>1. || y>1.)
{
//nerror("1 or 2 poles along real axis in appell_integ\n");
printf("Can't accept x>1 or y>1 in appell_integ\n");
}
/*else
{
integ=0.;
//if(a==1.) integ+=1.;
for(i=1; i<1000000; i++)
{
t=0.000001*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
if(i%2==1) integ+=4.*f;
else integ+=2.*f;
}
if(a-1.<0.)
{
f=appellf1_integrand(0.000001,a,b,c,d,x,y);
}
else
{
f=appellf1_integrand(0.,a,b,c,d,x,y);
}
integ+=f;
if(x!=1. && y!=1. && d-a-1.>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x!=1. && y!=1. && d-a-1.<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x==1. && y!=1. && d-a-1.-b>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x==1. && y!=1. && d-a-1.-b<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x!=1. && y==1. && d-a-1.-c>=0.)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x!=1. && y==1. && d-a-1.-c<0.)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
if(x==1. && y==1. && d-a-1.-b-c>=0)
{
f=appellf1_integrand(1.,a,b,c,d,x,y);
}
if(x==1. && y==1. && d-a-1.-b-c<0)
{
f=appellf1_integrand(0.999999,a,b,c,d,x,y);
}
integ+=f;
integ*=0.000001/3;
integ/=beta(a,d-a);
printf("appell=%18.15lf\n", integ);
}*/
else
{
integ[0]=0.;
for(n=10; n<=27; n++)
{
integ[n]=0.;
h=pow(10.,-n);
for(i=1; i<=10; i++)
{
t=h*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
//printf("f(%lf)=%lf\n", t, f);
if(i==1 || i==10) integ[n]+=f;
else if(i%2==1) integ[n]+=2.*f;
else if(i%2==0) integ[n]+=4.*f;
t=1.-(10.-i)*h;
f=appellf1_integrand(t,a,b,c,d,x,y);
//printf("f(%lf)=%lf\n", t, f);
if(i==1 || i==10) integ[n]+=f;
else if(i%2==1) integ[n]+=2.*f;
else if(i%2==0) integ[n]+=4.*f;
}
integ[n]*=h/3.;
integ[0]+=integ[n];
}
h=pow(10.,-9);
lim=pow(10.,9);
integ[9]=0.;
for(i=1; i<lim; i++)
{
t=h*i;
f=appellf1_integrand(t,a,b,c,d,x,y);
if(i==1 || i==lim-1) integ[9]+=f;
else if(i%2==1) integ[9]+=2.*f;
else if(i%2==0) integ[9]+=4.*f;
}
integ[9]*=h/3.;
integ[0]+=integ[9];
integ[0]/=beta(a,d-a);
printf("appell=%18.15lf\n", integ[0]);
}
return 0;
}
Codice: Seleziona tutto
#include<math.h>
double appellf1_integrand(double t, double a, double b, double c, double d, double x, double y)
{
double f;
if(t==1. && x==1. && y!=1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-b)*pow(1.-y*t,-c);
if(t==1. && x!=1. && y==1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-c)*pow(1.-x*t,-b);
if(t==1. && x==1. && y==1.) f=pow(t,a-1.)*pow(1.-t,d-a-1.-b-c);
else f=pow(t,a-1.)*pow(1.-t,d-a-1.)*pow(1.-x*t,-b)*pow(1.-y*t,-c);
return f;
}
Codice: Seleziona tutto
#include<math.h>
double beta(double x, double y)
{
double B;
B=exp(log(gamma(x))+log(gamma(y))-log(gamma(x+y)));
return B;
}
0.5, 1, 0.5, 2.75, -399, 0.12
0.5, -2, 1, 1, 0.12, -48
Re: Funzione ipergeometrica di Gauss in C
Non puoi gestire così le singolarità, dovresti integrare in campo complesso lungo un percorso che le evita, ad esempio composto da tratti rettilinei e da semicirconferenze che ci girano intorno, e poi vedere se esiste un limite dell'integrale per il raggio delle semicirconferenze che tende a zero.
Però qui trovi del codice fortran già fatto e in questo documento gli autori spiegano i metodi usati.
Questi sono altri autori che hanno fatto la stessa cosa, ma con un approssimante di Padé (razionale) per gestire le singolarità; non si capisce se abbiano reso pubblico il loro codice.
Però qui trovi del codice fortran già fatto e in questo documento gli autori spiegano i metodi usati.
Questi sono altri autori che hanno fatto la stessa cosa, ma con un approssimante di Padé (razionale) per gestire le singolarità; non si capisce se abbiano reso pubblico il loro codice.
Re: Funzione ipergeometrica di Gauss in C
Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il linkbite ha scritto:... non si capisce se abbiano reso pubblico il loro codice.
Ciao
Re: Funzione ipergeometrica di Gauss in C
No, il tuo link si riferisce alla traduzione (o meglio wrapping) in R del codice di Colavecchia, che è il primo dei due link che ho postato io.BlueEyes ha scritto:Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il linkbite ha scritto:... non si capisce se abbiano reso pubblico il loro codice.
Ciao
Nel mio secondo link gli autori sono altri (Cuyt, Driver, Tan e Verdonk). Vi si descrive un algoritmo che costoro avranno (si spera ) provato scrivendo del codice, di cui però io non ho trovato traccia.
Re: Funzione ipergeometrica di Gauss in C
Sì, è vero. C'è però dentro (il link) un file pdf, da dove ho potuto ricavare il seguente codice Fortran, che ritengo possa essere utile allo studente che ha aperto il topic.bite ha scritto:... il tuo link si riferisce alla traduzione (o meglio wrapping) in R del codice di Colavecchia...
Codice: Seleziona tutto
c***********************************************************************
c
c program : complex hypergeometric function
c
c notation : F(a,b;c;z)
c
c usage: The program has been written to be used as a function
c where the user supplies the hypergeometric parameters
c (a,b, and c), the independent variable (z).
c
c warning: This code is under construction. Warning messages
c will appear in regions where the code has not yet
c been completed.
c
c written by:
c
c Robert C. Forrey
c Institute for Theoretical Atomic and Molecular Physics
c Harvard-Smithsonian Center for Astrophysics
c 60 Garden Street Cambridge, MA 02138
c rforrey@cfa.harvard.edu
c
c***********************************************************************
c
c function name - chyp(a,b,c,z)
c
c computation
c performed - calculates the hypergeometric function.
c
c arguments z - the independent variable of the hypergeometric
c function.
c
c a,b,c - parameters of the hypergeometric function.
c
c precision - complex*16
c
c language - fortran
c
c***********************************************************************
complex*16 function chyp(a,b,c,z)
implicit none
real*8 zero,one,pi,tmp1,tmp2,x,y,r,phi
parameter (zero=0.d0,one=1.d0)
integer flag,n
complex*16 i,a,b,c,z,w,f1,f2,cgamm,csum,
# ctmp,ctmp1,ctmp2,ctmp3,ctmp4,ogamm
pi=dacos(-1.d0)
i=dcmplx(0.d0,1.d0)
c Set error messages
tmp1=dreal(a-b)
tmp2=dimag(a-b)
n=nint(tmp1)
tmp1=tmp1+tmp2-n
if(tmp1.eq.zero)then
write(*,*)'transformation error in chyp'
return
endif
tmp1=dreal(c-a-b)
tmp2=dimag(c-a-b)
n=nint(tmp1)
tmp1=tmp1+tmp2-n
if(tmp1.eq.zero)then
write(*,*)'transformation error in chyp'
return
endif
c Handle the special cases when z=0 or z=1
tmp1=cdabs(z)
if (tmp1.eq.zero) then
chyp=dcmplx(one)
return
endif
ctmp=dcmplx(one)
if (z.eq.ctmp) then
c chyp=cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
chyp=cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
return
endif
c Transform to a new variable w whose magnitude is the
c lowest possible value that lies between 0 and 1.
tmp1=cdabs(z)
flag=3
w=1/(1-z)
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=1
endif
w=z/(z-1)
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=2
endif
w=1-z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=4
endif
w=1-1/z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=5
endif
w=1/z
tmp2=cdabs(w)
if(tmp2 .lt. tmp1) then
tmp1=tmp2
flag=6
endif
if(tmp1.ge.one)then
write(*,*)'error in transformation'
return
endif
write(20,*)'flag=',flag
c Compute the hypergeometric function in z via the transformation
c theory and the series representation in w.
if (flag .eq. 1)then
w=1/(1-z)
f1=csum(a,c-b,a-b+1,w)
f2=csum(b,c-a,b-a+1,w)
x=dreal(1-z)
y=dimag(1-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 1'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)/cgamm(b)/cgamm(c-a)
c # +f2*ctmp2*cgamm(c)*cgamm(a-b)/cgamm(a)/cgamm(c-b)
chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)*ogamm(b)*ogamm(c-a)
# +f2*ctmp2*cgamm(c)*cgamm(a-b)*ogamm(a)*ogamm(c-b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-w**a
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-w**b
elseif (flag .eq. 2)then
w=z/(z-1)
chyp=(1-w)**a*csum(a,c-b,c,w)
elseif (flag .eq. 3)then
w=z
chyp=csum(a,b,c,w)
elseif (flag .eq. 4)then
w=1-z
f1=csum(a,b,a+b-c+1,w)
f2=csum(c-a,c-b,c-a-b+1,w)
x=dreal(w)
y=dimag(w)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 4'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp2=ctmp**(c-a-b)
c chyp=f1*cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
c # +f2*ctmp2*cgamm(c)*cgamm(a+b-c)/cgamm(a)/cgamm(b)
chyp=f1*cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
# +f2*ctmp2*cgamm(c)*cgamm(a+b-c)*ogamm(a)*ogamm(b)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-w**(c-a-b)
elseif (flag .eq. 5)then
w=1-1/z
f1=csum(a,a-c+1,a+b-c+1,w)
f2=csum(c-a,1-a,c-a-b+1,w)
x=dreal(z)
y=dimag(z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.lt.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 5'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(a-c)
x=dreal(1-z)
y=dimag(1-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.le.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 5'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp2=ctmp2*ctmp**(c-a-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(c-a-b)/cgamm(c-a)/cgamm(c-b)
c # +f2*ctmp2*cgamm(c)*cgamm(a+b-c)/cgamm(a)/cgamm(b)
chyp=f1*ctmp1*cgamm(c)*cgamm(c-a-b)*ogamm(c-a)*ogamm(c-b)
# +f2*ctmp2*cgamm(c)*cgamm(a+b-c)*ogamm(a)*ogamm(b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-z**(-a)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-z**(a-c)*(1-z)**(c-a-b)
elseif (flag .eq. 6)then
w=1/z
f1=csum(a,a-c+1,a-b+1,w)
f2=csum(b-c+1,b,b-a+1,w)
x=dreal(-z)
y=dimag(-z)
if(x.gt.zero .and. y.ge.zero)then
phi=datan(dabs(y/x))
elseif(x.lt.zero .and. y.gt.zero)then
phi=pi-datan(dabs(y/x))
elseif(x.lt.zero .and. y.le.zero)then
phi=-pi+datan(dabs(y/x))
elseif(x.gt.zero .and. y.lt.zero)then
phi=-datan(dabs(y/x))
elseif(x.eq.zero .and. y.gt.zero)then
phi=pi/2
elseif(x.eq.zero .and. y.lt.zero)then
phi=-pi/2
else
write(*,*)'error in transformation 6'
return
endif
r=dsqrt(x**2+y**2)
ctmp=r*cdexp(i*phi)
ctmp1=ctmp**(-a)
ctmp2=ctmp**(-b)
c chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)/cgamm(b)/cgamm(c-a)
c # +f2*ctmp2*cgamm(c)*cgamm(a-b)/cgamm(a)/cgamm(c-b)
chyp=f1*ctmp1*cgamm(c)*cgamm(b-a)*ogamm(b)*ogamm(c-a)
# +f2*ctmp2*cgamm(c)*cgamm(a-b)*ogamm(a)*ogamm(c-b)
c write(*,'(2(2 e13.6,2x))')ctmp1,ctmp1-(-z)**(-a)
c write(*,'(2(2 e13.6,2x))')ctmp2,ctmp2-(-z)**(-b)
endif
c write(20,'(i2,2x,2(e13.6,2x))')flag,cdabs(z),cdabs(w)
return
end
c*******************************************************************
c
c234567
complex*16 function csum(a,b,c,z)
implicit none
integer n
real*8 test,eps
complex*16 a,b,c,z,num,den,tmp,sum
eps=1.d-14
tmp=dcmplx(1.d0)
sum=dcmplx(0.d0)
do 10 n=1,300
sum=sum+tmp
test=cdabs(tmp/sum)
if(test.gt.eps)then
num=(a+n-1)*(b+n-1)
den=(c+n-1)*n
tmp=tmp*num*z/den
else
go to 20
endif
10 continue
write(*,*)'warning in csum: sum not completely converged,
#|z|=',cdabs(z)
20 continue
csum=sum
return
end
c*******************************************************************
c
c function name - ogamm
c
c computation
c performed - calculates one over the gamma function
c
c usage - ogamm(z)
c
c argument z - any complex number
c
c precision - double
c
c language - fortran
c
c*******************************************************************
complex*16 function ogamm(z)
implicit none
real*8 one
parameter (one=1.d0)
complex*16 z,ztmp,coeff,cgamm
ztmp=z
coeff=dcmplx(one)
1 if (dreal(ztmp).lt.one) then
coeff=coeff*ztmp
ztmp=ztmp+one
go to 1
endif
ogamm=coeff/cgamm(ztmp)
return
end
c*******************************************************************
c
c function name - cgamm
c
c computation
c performed - calculates the complex gamma function
c
c usage - cgamm(z)
c
c argument z - any complex number
c
c precision - double
c
c language - fortran
c
c*******************************************************************
complex*16 function cgamm(z)
implicit none
integer k
real*8 c(0:15),pi,one
parameter (one=1.d0)
complex*16 z,ztmp,sum,tmp,const
pi=acos(-1.d0)
call coeff(c,15)
ztmp=z
const=dcmplx(one)
1 if (dreal(ztmp).lt.one) then
const=const*ztmp
ztmp=ztmp+one
go to 1
endif
sum=cmplx(0.d0)
tmp=cmplx(1.d0)
do 10 k=0,15
sum=sum+c(k)*tmp
tmp=tmp*(ztmp-k-1)/(ztmp+k)
10 continue
cgamm=sum*dsqrt(2*pi)*(ztmp+9.d0/2)**(ztmp-.5d0)
# *exp(-ztmp-9.d0/2)/const
return
end
c*******************************************************************
c
c subroutine name - coeff
c
c computation
c performed - tabulates the gamma function coefficients
c (see luke vol.2 p.304)
c
c usage - call coeff(c,n)
c
c arguments c - the array (output) which contains the
c coefficients.
c
c n - the dimension (input) of the array 'c'.
c
c precision - double
c
c language - fortran 77
c
c*******************************************************************
subroutine coeff(c,n)
real*8 c(0:n)
c(0) = 41.624436916439068d0
c(1) =-51.224241022374774d0
c(2) = 11.338755813488977d0
c(3) = -0.747732687772388d0
c(4) = 0.008782877493061d0
c(5) = -0.000001899030264d0
c(6) = 0.000000001946335d0
c(7) = -0.000000000199345d0
c(8) = 0.000000000008433d0
c(9) = 0.000000000001486d0
c(10)= -0.000000000000806d0
c(11)= 0.000000000000293d0
c(12)= -0.000000000000102d0
c(13)= 0.000000000000037d0
c(14)= -0.000000000000014d0
c(15)= 0.000000000000006d0
return
end
Re: Funzione ipergeometrica di Gauss in C
Quella che hai riportato però mi sembra la 2F1, non la F1 di Appell.BlueEyes ha scritto:C'è però dentro (il link) un file pdf, da dove ho potuto ricavare il seguente codice Fortran, che ritengo possa essere utile allo studente che ha aperto il topic.
Re: Funzione ipergeometrica di Gauss in C
giuseppe morello ha scritto: 2, sull'algoritmo: per a=-1, b=5, c=3, z=4 e anche in alcuni altri casi, il codice non converge, ritenendo che hnext risulti troppo piccolo. Per gli stessi valori wolphram non mostra problemi. In effetti, l'ho fatto stampare e risulta nullo almeno fino alla 15^a cifra, mentre hmin=0.0001 (anche se non ho ancora trovato dove viene assegnato questo valore).
Usando questa libreria di Python, si ottiene la convergenza della funzione ipergeometrica, che coincide con quella di Wolfram (-5.66666). Ecco l'output
Codice: Seleziona tutto
$ python
Python 2.7.1+ (r271:86832, Apr 11 2011, 18:05:24)
[GCC 4.5.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import mpmath
>>> float(mpmath.hyp2f1(float(-1), float(5), float(3), float(4)))
-5.666666666666667
>>>
EDIT. Con z avente la parte immaginaria non nulla si ottiene:
Codice: Seleziona tutto
>>> import mpmath
>>> mpmath.hyp2f1(-1, 5, 3, 4+3j)
mpc(real='-5.666666666666667', imag='-5.0')
-
- Prode Principiante
- Messaggi: 68
- Iscrizione: giovedì 5 aprile 2012, 9:56
Re: Funzione ipergeometrica di Gauss in C
Ok, grazie ancora per i contributi. Sto cercando di studiarmeli per vedere quale riesco a riscrivere.
Eventualmente, anche se al momento è una soluzione che non mi piace moltissimo, è possibile linkare librerie scritte in altri linguaggi, per esempio python o fortran, da un sorgente scritto in C?
Eventualmente, anche se al momento è una soluzione che non mi piace moltissimo, è possibile linkare librerie scritte in altri linguaggi, per esempio python o fortran, da un sorgente scritto in C?
Re: Funzione ipergeometrica di Gauss in C
Fortran sicuramente sì e probabilmente anche python, anche se mi sembra un po' un carcabaggio che un linguaggio compilato ne chiami uno interpretato.
Ma vedo che in questo momento c'è un esperto di python che sta visualizzando questa pagina e spero che intervenga
Ma vedo che in questo momento c'è un esperto di python che sta visualizzando questa pagina e spero che intervenga
Chi c’è in linea
Visualizzano questa sezione: 0 utenti iscritti e 11 ospiti