Funzione ipergeometrica di Gauss in C

Vedi allegato.

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).

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)

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.

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).

Vero, ma non so come rimediare. hmin viene passato come argomento alla odeint, ma cambiarlo non risolve.

Ho visto che la GNU scientific library ha la funzione ipergeometrica (vedi http://www.gnu.org/software/gsl/manual/ ... tions.html); ti conviene provarla.

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:

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


Compilando con

Codice: Seleziona tutto

gcc -c -lm ipergeom.c -o ipergeom.out

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

Codice: Seleziona tutto

pkg-config --libs gsl_sf_hyperg.h

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

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


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.

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


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.

Da dove ti salta fuori hyperg_2F1_series? I nomi di tutte le funzioni della libreria gsl cominciano con gsl_

Svista.

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.

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).

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

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.

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

Temo sia impossibile. odeint sta per Ordinary DIfferential Equation INTegrator, mentre le equazioni differenziali costitutive delle funzioni di Appell sono alle derivate parziali.
Prima di fasciarti la testa, però, verifica qual è il raggio di convergenza della serie, può anche darsi che ti basti.

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.

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.

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.

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.

Mi sa che questa e' l'idea giusta . Il dominio di convergenza della serie non era sufficiente, quello dell'integrale pare di si' invece.
Al momento l'ho calcolato semplicemente con il metodo di Simpson, devo solo perfezionarlo per avere una convergenza piu' rapida. Ti aggiornero'.

Al momento l'ho calcolato semplicemente con il metodo di Simpson, devo solo perfezionarlo per avere una convergenza piu' rapida. Ti aggiornero'.

Attenzione però che l'integrando può avere delle singolarità a seconda dei valori dei parametri.

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à)

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


appellf1_integrand.c

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


beta.c

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


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

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.

... non si capisce se abbiano reso pubblico il loro codice.

Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il link
Ciao

... non si capisce se abbiano reso pubblico il loro codice.

Sì, qualcuno lo ha recentemente convertito dal Fortran in R (pacchetto statistico), con la collaborazione degli autori. Ecco il link
Ciao

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.

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.

... il tuo link si riferisce alla traduzione (o meglio wrapping) in R del codice di Colavecchia...

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.

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                                                         


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.

Quella che hai riportato però mi sembra la 2F1, non la F1 di Appell.

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


Ciao

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')


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?

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