Funzione ipergeometrica di Gauss in C

No, non sono un esperto di python, ma un newbie, che ha trovato qui il modo di risolvere quella funzione ipergeometrica. L'algoritmo riesce a farla convergere anche per valori di z=1*(10)^100, come indicato nello screenshot allegato, di cui ho verificato i risultati.
Come ho scritto prima, la mia discussione è mostra solo che altri linguaggi (diversi dal C) riescono ad arrivare all'obiettivo.
Ciao

No, non sono un esperto di python, ma un newbie

Non mi riferivo a te (non ti conosco e quindi non posso sapere se sei esperto di python) ma a un'altra persona che se l'è battuta...

Allora, mi sembra di capire che il codice di Colavecchia in tutti i casi in cui max(|x|,|y|)>=1 si riconduce alla zona di convergenza della serie tramite una trasformazione (anche se all'interno della zona stessa non sempre usa la serie ...).

Vedo che l'unico program è driver.f90 quindi per come è scritto serve solo a produrre tabelle, ma questo non è un problema.
Stavo provando a riscrivere in C i vari pezzi, a partire dal get_transformation2.f90, ma forse prima di fare tutta questa fatica sarebbe forse opportuno controllare che funzioni per i valori che interessano a me. Per farlo dovrei costruire tutti gli header e aggiungerli al driver prima di compilare?

Fortran sicuramente sì e probabilmente anche python, anche se mi sembra un po' un carcabaggio che un linguaggio compilato ne chiami uno interpretato.

Anche a me non piace proprio come soluzione, ma è quasi 1 mese che giro attorno a sto programma senza riuscire a concludere. Inoltre, mi pare di avere notato recentemente che anche le funzioni utilizzate dall'originale, ovvero integrali ellittici e altro, non convergano bene. Tento ancora qualcosa però prima di arrendermi.

Il codice fortran di Colavecchia dove l'hai preso? Ho provato a seguire fino in fondo il link che ti avevo postato e mi ha chiesto dei soldi. Sono subito fuggito terrorizzato

Però qui trovi del codice fortran già fatto e in questo documento gli autori spiegano i metodi usati.

Ora non ricordo esattamente, ma mi sembra che cliccando su "qui" e poi "download" ti chiede di accettare i termini della licenza e possibilmente di dare la tua e-mail, ma non chiede soldi.

A me ha chiesto 80 euro
È vero che non ho voluto dare l'email, ma specificava che era opzionale.

In caso vedo se posso dartene copia.

Intanto una brutta sorpresa per me: Wolphram Alpha non riconosce la funzione F2 di Appell. Siccome questa serve per il calcolo della F1 secondo il metodo di Colavecchia, mi sarebbe utile verificare se quel che ho scritto funziona.
Penso non sia difficile con un eseguibile in python oppure una lista di comandi python da terminale (mai usato). Grazie.

P.s.: Ho trovato questo http://www.lamprechts.de/gerd/php/Rechn ... nktion.php , non so quanto affidabile; e' in accordo con i miei risultati, ma solo per le prime 3 cifre circa. Sospetto sia comunque piu' preciso lui di me.

EDIT: Sono riuscito a essere piu' preciso, ora vado abbastanza d'accordo con i risultati del link, quindi lo ritengo abbastanza affidabile. A volte pero' fallisce, quindi se si trova di meglio ben venga.

In caso vedo se posso dartene copia.

Prova a mandarmela, ci do un'occhiata.

Allora, mi ritrovo con un problema di cui non riesco a trovare la causa.

Ho aggiunto alcune nuove funzioni a quelle precedenti, in modo che dovrei potere sempre calcolare la funzione F1 di Appell, anche se non sempre come suggerito nell'articolo.
In pratica il problema e' che quando lancio il programma ci sono parti di codice che si attivano ma non dovrebbero.

In particolare, cosa vorrei che accadesse:
test.c, contiene il main, do' in input i valori dei 6 parametri della F1, chiama quindi appellf1.c.
appellf1.c, a seconda dei parametri decide come fare il calcolo. Nel mio caso semplicemente decide di calcolare la serie, chiamando in causa appf1_hypser.c.
appf1_hypser.c calcola la serie a doppia entrata come da definizione e restituisce il risultato al livello precedente, che a sua volta lo restituisce al file iniziale.
test.c stampa il risultato

Cosa invece accade (ho inserito vari printf per controllare bene la situazione):
test.c ok
appellf1.c ok
appf1_hypser.c ok, ma opera 3 volte cambiando i valori dei parametri di volta in volta.
Quindi senza motivo si attivano anche horng2.c e appellf2.c, senza essere chiamate da nessuno, perche' non mi ritrovo in questo caso.
Finalmente test.c stampa il suo risultato, che pero' non e' quello giusto.

Quale potrebbe essere il problema?

Ho pensato che possa entrarci il makefile. In sintesi riscontro 2 problemi:
funzioni che si attivano senza essere chiamate;
malfunzionamento del passaparola tra più funzioni in serie.
C'è forse qualche opzione da aggiungere nel makefile per gestire questi problemi?
Sorge un conflitto se una funzione restituisce una variabile ad una superiore in una variabile con lo stesso nome? (entrambe variabili locali). In realtà penso di avere sempre fatto così, ma magari mi sbaglio.

Non è che appf1_hypser è ricorsiva per caso?

Non sembra un problema da makefile.

I nomi delle variabili locali in funzioni diverse sono in scope diversi e non vanno in conflitto.

Non direi, c'è un doppio ciclo for, ma il return è fuori da tutto.

Codice: Seleziona tutto

#include<math.h>
#include<complex.h>
#include<stdio.h>
#include"appf1_hypser.h"


double complex appf1_hypser (double complex a, double complex b, double complex c, double complex d, double complex x, double complex y)
{
   int n, m;
   double complex aa, bb, cc, dd, aaa, bbb, ddd, fac, facc, F1;
   
    printf("a=%lf b=%lf c=%lf d=%lf x=%lf y=%lf\n", creal(a), creal(b), creal(c), creal(d), creal(x), creal(y));
   
   fac = 1.0 + 0.0*I;
   F1 = 0.0 + 0.0*I;
   aa = a-1.;
   bb = b-1.;
   cc = c-1.;
    dd = d-1.;
   for (n = 0; n <= 1000; n++)
   {
        if(n>0)
        {
            fac*=(aa/dd)*cc;
            fac*=y/n;
        }
        facc=fac;
        aaa=aa+1.;
        ddd=dd+1.;
        bbb=bb+1.;
        F1+=fac;
        for(m=1; m<=1000; m++)
        {
            facc*=(aaa/ddd)*bbb;
            facc*=x/m;
            F1+=facc;
            aaa+=1.;
            bbb+=1.;
            ddd+=1.;
            //printf("serie=%18.15f+%18.15f I\n", creal(F1), cimag(F1));
        }
      aa += 1.;
      cc += 1.;
      dd += 1.;
   }
    printf("AppellF1serie=%18.15f+%18.15f I\n", creal(F1), cimag(F1));
    return F1;
}


Senza vedere il codice è difficile indovinare che succede.

Finalmente ho risolto questo problema, anche se non ho capito benissimo perché accadeva.
In pratica il file appellf1.c aveva la seguente struttura:

Codice: Seleziona tutto

#include<math.h>
#include<complex.h>
#include<stdio.h>
#include"appellf1.h"
#include"appf1_hypser.h"
#include"horng2.h"

double complex appellf1 (double complex a, double complex b, double complex c, double complex d, double complex x, double complex y)
{
    double complex F1, aus;
    double u, w, dist, distmax;
    int flag;

    distmax = 1.;
    if(cabs(x) >= cabs(y))  dist = cabs(x);
    else  dist = cabs(y);
   
    if(dist < distmax)
    {
        flag = 0;
    }
    else
    {
          ...........una serie di condizioni che a seconda dei valori di x e y scelgono una flag
    }

    if(flag==0)
    {
        F1 = appf1_hypser(a,b,c,d,x,y);
        //printf("FLAG=%d\n", flag);
        //printf("a=%lf b=%lf c=%lf d=%lf x=%lf y=%lf\n", creal(a), creal(b), creal(c), creal(d), creal(x), creal(y));
    }
    ...........a seconda del valore di flag vengono eseguite delle operazioni su x e y e vengono chiamate diverse funzioni

    return F1;
}


Mettendo il return F1 in fondo a ogni if(flag==...) anzichè in fondo al file, fuori da tutto, funziona correttamente.

Sto cercando qualche formula relativa preferibilmente alla F1 di Appell, ma eventualmente anche la F2 di Appell o la G2 di Horn, che mi dica come comportarmi in qualche caso particolare.
In pratica ho notato che mi trovo spesso a che fare con:
F1(0.5, b, 0.5, 0.5, x, y) fuori dalla zona di convergenza della serie.
In questo caso il codice che uso richiede di calcolare la funzione:
G2(a, 1, 1, −0.5, x, y) (non intendo che x e y siano uguali a quelli di prima), la quale a sua volta mi chiede di calcolare:
F2(0.5, b, 1, 0, 1.5, x, y).
Quel d=0 crea problemi perché è associato a un simbolo di Pochhammer che si annulla al denominatore, facendo divergere tutto.
Tuttavia, ho verificato con Wolphram Alpha che la F1 di partenza converge. Devo quindi trovare una formula particolare che mi permetta di trattare casi tipo questo a parte.
Le formule che ho trovato online sono però troppo restrittive, in quanto sostanzialmente mi vincolano più parametri del dovuto. Non esiste magari una fonte specializzata sulle funzioni di Appell o su quelle di Horn? (Lo stesso handbook of mathematical functions di Abramowitz and Stegun, mi sembra non tratti queste funzioni)