Funzione ipergeometrica di Gauss in C

Linguaggi di programmazione: php, perl, python, C, bash e tutti gli altri.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Salve a tutti,
ho bisogno di un codice in C che sia capace di calcolare i valori della funzione ipergeometrica di Gauss 2F1. Ho provato a copiare quello che ho trovato in un Numerical Recipes in C (Second Edition), ma mi dà degli errori, molti dei quali sembrano legati al fatto che non ho la libreria "complex.h". Sapete dove posso scaricarla?
In alternativa esistono dei codici più semplici?
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

Lo header (non libreria!) complex.h è di sistema e dovresti semplicemente scrivere

Codice: Seleziona tutto

#include <complex.h>
per vederlo.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Grazie del suggerimento. Nel libro suggeriva di scrivere

Codice: Seleziona tutto

#include"complex.h"
non avevo mai visto l'uso delle virgolette e non so quale sia la differenza.
Anche aggiustando come hai detto, continuano a restare molti (quasi tutti) degli errori precedenti; riporto un pezzo di codice per mostrarne uno:

Codice: Seleziona tutto

#include<math.h>
#include<complex.h>
#include"nrutil.h"
#define EPS 1.0e-6

fcomplex aa, bb, cc, z0, dz;
Compilando mi risulta errore al rigo 6:
error: expected '=', ',', ';', 'asm' or '__attribute__' before 'aa'
Perché questo errore? Pare non riconosca fcomplex.
Ultima modifica di giuseppe morello il mercoledì 18 aprile 2012, 16:05, modificato 1 volta in totale.
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

Se il libro suggerisce di scrivere "complex.h" è possibile che non intenda lo header di sistema (che in effetti se non sbaglio c'è solo dal C99 mentre il Numerical Recipes è piuttosto vecchiottto) ma che si riferisca ad un suo header che trovi da qualche parte nel libro o nel cd.

Per tagliare la testa al toro basta vedere se viene utilizzata la keyword _Complex. In questo caso si riferisce al complex.h di sistema.

In questo momento non ho sottomano la mia copia di Numerical Recipes e non posso verificare. Se hai voglia di postare qui il codice ci do un'occhiata.
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

Una breve sgooglata ha scavato fuori http://www.nr.com/public-domain.html

dove trovi complex.h e anche complex.c, che dovrai compilare e linkare assieme al resto.

Se cerchi bene sono convinto che ci siano anche sul libro, con tutte le indicazioni del caso.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Grazie di nuovo. Poi questi complex.h l'ho trovato e l'ho copiato in un file nella cartella di lavoro. Mi risultano pero' esattamente gli stessi errori anche in complex.h.
Comunque in fondo ho messo il sorgente che ho tratto dal Numerical Recipes, se avessi la pazienza di darci un'occhiata (e' un po' elaborato). A me basterebbe lavorare con i numeri reali, questo usa pure i complessi; sto cercando di trovare un codice piu' semplice che eviti i complessi.

Codice: Seleziona tutto

#include<math.h>
#include"complex.h"
#include"nrutil.h"
#define EPS 1.0e-6 //Accuracy parameter

fcomplex aa, bb, cc, z0, dz; //Communicates with hypdrv

int kmax, kount; //Used by odeint
float *xp, **yp, dxsav;

fcomplex hypgeo(fcomplex a, fcomplex b, fcomplex c, fcomplex z)
/*Complex hypergeometric function 2F1 for complexa, b, c, and z, by direct integration of the hypergeometric equation in the complex plane. The branch cut is taken to lie along the real axis, Re z > 1.*/
{
		void bsstep(float y[], float dydx[], int nv, float *xx, float htry, float eps, float yscal[], float *hdid, float *hnext, void (*derivs)(float, float[], float[]);
		void hypdrv(float s, float yy[], float dyyds[]);
		void hypser(fcomplex a, fcomplex b, fcomplex c, fcomplex z, fcomplex *series, fcomplex *deriv);
		void odeint(float ystart[], int nvar, float x1, float x2, float eps, float h1, float hmin, int *nok, int *nbad, void(*derivs)(float, float[], float[]), void(*rkqs)(float[], float[], int, float *, float, float, float[], float *, float *, void(*)(float, float[], float[])));
		int nbad, nok;
		fcomplex ans, y[3];
		float *yy;
		
		kmax=0;
		if(z.r*z.r+z.i*z.i <= 0.25)
		{
				hypser(a,b,c,z,&ans,&y[2]);
				return ans;
		}
		else if(z.r<0.0) z0=Complex(-0.5,0.0);
		else if(z.r<=1.0) z0=Complex(0.5,0.0);
		else z0=Complex(0.0,z.i>0.0?0.5:-0.5);
		aa=a;
		bb=b;
		cc=c;
		dz=Csub(z,z0);
		hypser(aa,bb,cc,z0,&y[1],&y[2]);
		yy=vector(1,4);
		yy[1]=y[1].r;
		yy[2]=y[1].i;
		yy[3]=y[2].r;
		yy[4]=y[2].i;
		odeint(yy,4,0.0,1.0,EPS,0.1,0.0001,&nok,&nbad,hypdrv,bsstep);
/*The arguments to odeint are the vector of indipendent variables, its lenght, the starting and ending values of the dependent variable, the accuracy parameter, an initial guess for stepsize, a minimum stepsize, the (returned) number of good and bad steps taken, and the names of the derivative routine and the (here Bulirsch-Stoer) stepping routine.*/
		y[1]=Complex(yy[1],yy[2]);
		free_vector(yy,1,4);
		return y[1];
}


#define ONE Complex(1.0,0.0)

void hypser(fcomplex a, fcomplex b, fcomplex c, fcomplex z, fcomplex *series, fcomplex *deriv)
/*Returns the hypergeometric series 2F1 and its derivative, iterating to machine accuracy. For |z|<=1/2 the convergence is quite rapid.*/
{
		void nrerror(char error_text[]);
		int n;
		fcomplex aa, bb, cc, fac, temp;

		deriv->r=0.0;
		deriv->i=0.0;
		fac=Complex(1.0,0.0);
		temp=fac;
		aa=a;
		bb=b;
		cc=c;
		for(n=1; n<=1000; n++)
		{
				fac=Cmul(fac,Cmul(aa,Cdiv(bb,cc)));
				deriv->r+=fac.r;
				deriv->i+=fac.i;
				fac=Cmul(fac,RCmul(1.0/n,z));
				*series=Cadd(temp,fac);
				if(series->r == temp.r && series->i == temp.i) return;
				temp=*series;
				aa=Cadd(aa,ONE);
				bb=Cadd(bb,ONE);
				cc=Cadd(cc,ONE);
		}
		nrerror("convergence failure in hypser");
}


//extern fcomplex aa,bb,cc,z0,dz;

void hypdrv(float s, float yy[], float dyyds[])
/*Computes derivatives for the hypergeometric equation, see text equation (5.14.4).*/
{
		fcomplex z, y[3], dyds[3];

		y[1]=Complex(yy[1],yy[2]);
		y[2]=Complex(yy[3],yy[4]);
		z=Cadd(z0,RCmul(s,dz));
		dyds[1]=Cmul(y[2],dz);
		dyds[2]=Cmul(Csub(Cmul(Cmul(aa,bb),y[1]),Cmul(Csub(cc,Cmul(Cadd(Cadd(aa,bb),ONE),z)),y[2])),Cdiv(dz,Cmul(z,Csub(ONE,z))));
		dyyds[1]=dyds[1].r;
		dyyds[2]=dyds[1].i;
		dyyds[3]=dyds[2].r;
		dyyds[4]=dyds[2].i;
}

/*Tratto da Numerical Recipes in C, second edition, Press, Teukolsky, Vetterling, Flannery; pagg.272-273*/
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

Il codice di NR che hai postato ha tre funzioni: hypgeo calcola 2F1 integrando l'equazione differenziale ipergeometrica di Eulero, hypser calcola 2F1 sommando fino ad un massimo di 1000 termini della serie, o meno se si converge prima entro l'epsilon macchina; infine hypdrv calcola le derivate, non mi è ben chiaro come (sarà spiegato nel testo).

Penso che hypser dovrebbe bastarti; te l'ho riscritta  facendo uso dello header complex.h di sistema:

Codice: Seleziona tutto

/* Tratto da Numerical Recipes in C, second edition,
	 Press, Teukolsky, Vetterling, Flannery; pagg.272-273
	 
	 Riscritto facendo uso di complex.h del C99
*/

#include <stdio.h>
#include <complex.h>

void hypser
(double complex a, double complex b, double complex c, double complex z,
double complex * series, double complex * deriv)

/* Returns the hypergeometric series 2F1 and its derivative,
   iterating to machine accuracy. For |z|<=1/2 the convergence is quite rapid.*/
   
{
		int n;
		double complex aa, bb, cc, fac, temp;
		double complex temp1, one;
		one = 1.0 + 0.0 * I;
		
		*deriv = 0.0 + 0.0 * I;
		fac = 1.0 + 0.0 * I;
		temp = fac;
		aa = a;
		bb = b;
		cc = c;
		for (n = 1; n <= 1000; n++)
		{
				fac *= (aa * (bb / cc));
				*deriv += fac;
				temp1 = 1.0 / n + 0.0 * I;
				fac *= temp1 * z;
				
				*series = temp + fac;
				if (*series == temp) return;
				temp = *series;
				aa += one;
				bb += one;
				cc += one;
		}
		printf ("convergence failure in hypser\n");
}


int main ()
{
	complex a, b, c, z, result, deriv;	
	a = 1.0 + 0.0 * I;
	b = 2.0 + 0.0 * I;
	c = 3.0 + 0.0 * I;
	z = 0.5 + 0.0 * I;
	hypser (a, b, c, z, &result, &deriv);
	printf ("%f + %fI\n", creal (result), cimag (result));
}

Verificala chiamandola con diversi valori e confrontando il risultato con quello che ti dà wolfram alpha (www.wolframalpha.com)
Ad esempio, per verificare il risultato della chiamata che ho messo nel main ho scritto

Codice: Seleziona tutto

Hypergeometric2F1[1,2,3,0.5+0I]
in wolfram alpha.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Grazie 1000, lo proverò al più presto e ti farò sapere.
Avevo visto che erano 3 funzioni, ma pensavo che venissero usate tutte e 3 a seconda dei casi.
Comunque se funziona quello che hai scritto è perfetto. (Io ne ho scritto uno più semplice per il caso reale, ma converge troppo troppo lentamente).
Avatar utente
BlueEyes
Entusiasta Emergente
Entusiasta Emergente
Messaggi: 1330
Iscrizione: giovedì 15 marzo 2012, 14:08

Re: Funzione ipergeometrica di Gauss in C

Messaggio da BlueEyes »

Ottimo codice, bite! Ciao
Immagine
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

L'ho solo riscritto  ;D

Un commento degli autori dice
For |z| 1/2, la convergenza non è lenta... proprio non c'è.

Ecco perché ci sono tre funzioni: per valori fuori dal raggio di convergenza della serie viene utilizzato l'altro metodo, la risoluzione dell'equazione differenziale... che però immagino sia piuttosto lento.

Quindi se servono valori |z| > 1/2 diventa indispensabile utilizzare anche la hypgeo e le funzioni che questa chiama, odeint e bsstep. Al massimo si può fare a meno della hypdrv se non serve la derivata.

Sulla base della mia traduzione non dovrebbe essere difficile riscrivere le funzioni che servono. Comunque se c'è bisogno di una mano sono qui.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

bite ha scritto: L'ho solo riscritto  ;D
Pero' mi risulta molto piu' leggibile. Quello originale ancora non mi e' chiaro.
Solo, da ignorante in materia, qual e' il vantaggio di definire *series e *deriv come puntatori, anziche' come tutte le altre variabili?

L'ho verificato stamattina, pare funzionare bene per |z|10, quindi mi sa che devo puntare su hypgeo.
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

giuseppe morello ha scritto: qual e' il vantaggio di definire *series e *deriv come puntatori, anziche' come tutte le altre variabili?
Il C, a differenza del C++, non consente di passare argomenti per riferimento, quindi passare un puntatore è l'unica possibilità se la funzione chiamata deve scrivere un risultato nell'argomento.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

bite ha scritto: Il C, a differenza del C++, non consente di passare argomenti per riferimento, quindi passare un puntatore è l'unica possibilità se la funzione chiamata deve scrivere un risultato nell'argomento.
Ecco, a volte confondo i due linguaggi.

Ho provato a tradurre hypgeo, basandomi sulla tua traduzione (mi manca qualcosa).

Codice: Seleziona tutto

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

#define EPS=1.0e-6 //Accuracy parameter

double complex aa, bb, cc, z0, dz;

int kmax, kount;

double *xp, **yp, dxsav;

double complex hypgeo (double complez a, double complex b, double complex c, double complex z)
{
		void bsstep(double y[], double dydx[], int nv, double *xx, double htry, double eps, double yscal[], double *hdid, double *hnext, void(*derivs)(double, double[], double[]);
		void hypdrv(double s, double yy[], double dyyds[]);
		void hypser(double complex a, double complex b, double complex c, double complex z, double complex *series, double complex *deriv);
		void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1, double hmin, int *nok, int *nbad, void(*derivs)(double, double[], double[]),
void(*rkqs)(double[], double[], int, double *, double, double, double[], double *, double *, void(*)(double, double[], double[])));
		int nbad, nok;
		double complex ans, y[3];
		double *yy;

		kmax=0;
		if(creal(z)*creal(z)+cimag(z)*cimag(z) <= 0.25)
		{
				hypser(a, b, c, z, &ans, &y[2]);
				return ans;
		}
		else if(creal(z)<0.0)
		{
				z0=-0.5+0.0*I;
		}
		else if(creal(z)<=1.0)
		{
				z0=0.5+0.0*I;
		}
		else
		{
				z0=...;
		}
		aa=a;
		bb=b;
		cc=c;
		dz=...;
		hypser(aa, bb, cc, z0, &y[1], &y[2]);
		yy=vector(1,4);
		yy[1]=creal(y[1]);
		yy[2]=cimag(y[1]);
		yy[3]=creal(y[2]);
		yy[4]=cimag(y[2]);
		odeint(yy, 4, 0.0, 1.0, EPS, 0.1, 0.0001, &nok, &nbad, hypdrv, bsstep);
		y[1]=...;
		free_vector(yy, 1, 4);
		return y[1];
}
Immagino di dovere trovare anche bsstep, odeint, e aggiungerli, assieme anche a hypser e hypdrv, affinche' possa funzionare.
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

hypdrv solo se ti serve la derivata, altrimenti puoi scartarlo.

La dichiarazione delle funzioni esterne dentro il corpo della funzione è un arcaismo. Per pulizia potresti scrivere i sorgenti in file .c separati: bsstep.c, odeint.c, hypser.c, hypgeo.c, ognuno con un corrispondente header: bsstep.h, odeint.h, hypser.h, hypgeo.h, contenente solo la dichiarazione; poi ad esempio hypgeo.c includerà hypgeo.h, hypser.h, odeint.h, bsstep.h e così via.

Numerical Recipes, essendo scritto da fortranisti, non si è mai adeguato agli array del C che partono dall'indice 0 e ha un suo inutilissimo allocatore/deallocatore di vettori per poterli usare con indici da 1 in su. Qesto poteva avere un senso una volta, quando la ram era rara e cara; adesso non vale più la pena, meglio semplificare gestendo gli array semplicemente con malloc/free e allocandoli con un posto in più, che resterà inutilizzato. Questo ti dovrebbe permettere di fare a meno dello header nrutil.h contribuendo alla pulizia del codice.

Invece di

Codice: Seleziona tutto

if(creal(z)*creal(z)+cimag(z)*cimag(z) <= 0.25)
puoi scrivere più semplicemente

Codice: Seleziona tutto

if(cabs(z) <= 0.5)
Quando hai completato la riscrittura posta che ci do un'occhiata.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Bene, a meno di errori (che sicuramente ci sono) ho finito. Sono 4 file:

hypgeo

Codice: Seleziona tutto

#include<math.h>
#include<stdio.h>
#include<complex.h>
//#include"nrutil.h"

#define EPS=1.0e-6 //Accuracy parameter

double complex aa, bb, cc, z0, dz;

int kmax, kount; //Usati da odeint

double *xp, **yp, dxsav;

double complex hypgeo (double complex a, double complex b, double complex c, double complex z)
{
		void bsstep(double y[], double dydx[], int nv, double *xx, double htry, double eps, double yscal[], double *hdid, double *hnext, void(*derivs)(double, double[], double[]);
void hypser(double complex a, double complex b, double complex c, double complex z, double complex *series, double complex *deriv);
		void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1, double hmin, int *nok, int *nbad, void(*derivs)(double, double[], double[]),
void(*rkqs)(double[], double[], int, double *, double, double, double[], double *, double *, void(*)(double, double[], double[])));
		int nbad, nok;
		double complex ans, y[3];
		double *yy;

		kmax=0;
		if(cabs(z) <= 0.5)
		{
				hypser(a, b, c, z, &ans, &y[2]);
				return ans;
		}
		else if(creal(z)<0.0)
		{
				z0=-0.5+0.0*I;
		}
		else if(creal(z)<=1.0)
		{
				z0=0.5+0.0*I;
		}
		else
		{
				z0=0.0+(cimag(z)>0.0 ? 0.5 : -0.5);
		}
		aa=a;
		bb=b;
		cc=c;
		dz=z-z0;
		hypser(aa, bb, cc, z0, &y[1], &y[2]);
		yy=vector(1,4);
		yy[1]=creal(y[1]);
		yy[2]=cimag(y[1]);
		yy[3]=creal(y[2]);
		yy[4]=cimag(y[2]);
		odeint(yy, 4, 0.0, 1.0, EPS, 0.1, 0.0001, &nok, &nbad, hypdrv, bsstep);
		y[1]=yy[1]+yy[2]*I;
		free_vector(yy, 1, 4);
		return y[1];
}
hypser

Codice: Seleziona tutto

/* Tratto da Numerical Recipes in C, second edition,
	 Press, Teukolsky, Vetterling, Flannery; pagg.272-273
	 
	 Riscritto facendo uso di complex.h del C99
*/

#include <stdio.h>
#include <complex.h>

void hypser
(double complex a, double complex b, double complex c, double complex z,
double complex * series, double complex * deriv)

/* Returns the hypergeometric series 2F1 and its derivative,
   iterating to machine accuracy. For |z|<=1/2 the convergence is quite rapid.*/
   
{
		int n;
		double complex aa, bb, cc, fac, temp;
		double complex temp1, one;
		one = 1.0 + 0.0 * I;
		
		*deriv = 0.0 + 0.0 * I;
		fac = 1.0 + 0.0 * I;
		temp = fac;
		aa = a;
		bb = b;
		cc = c;
		for (n = 1; n <= 1000; n++)
		{
				fac *= (aa * (bb / cc));
				*deriv += fac;
				temp1 = 1.0 / n + 0.0 * I;
				fac *= temp1 * z;
				
				*series = temp + fac;
				if (*series == temp) return;
				temp = *series;
				aa += one;
				bb += one;
				cc += one;
		}
		printf ("convergence failure in hypser\n");
}


int main ()
{
	complex a, b, c, z, result, deriv;
	printf("Inserisci a, b, c e z:\n");
	scanf("%lf %lf %lf %lf", &a, &b, &c, &z);	
	hypser (a, b, c, z, &result, &deriv);
	printf ("%18.15lf + %18.15lfI\n", creal (result), cimag (result));
}
odeint

Codice: Seleziona tutto

#include<math.h>
#define NRANSI
//#include "nrutil.h"
#define MAXSTP 10000
#define TINY 1.0e-30

extern int kmax, kount;
extern double *xp, **yp, dxsav;

void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1, double hmin, int *nok, int *nbad, void(*derivs)(double, double[], double[]), void(*rkqs)(double[], double[], int, double *, double, double, double[], double *, double *, void(*)(double, double[], double[])))
{
		int nstp, i;
		double xsav, x, hnext, hdid, h;
		double *yscal, *y, *dydx;

		yscal=vector(1,nvar);
		y=vector(1,nvar);
		dydx=vector(1,nvar);
		x=x1;
		h=SIGN(h1, x2-x1);
		*nok=(*nbad)=kount=0;
		for(i=1; i<=nvar; i++)
		{
				y[i]=ystart[i];
		}
		if(kmax>0)
		{
				xsav=x-dxsav*2.0;
		}
		for(nstp=1; nstp<=MAXSTP; nstp++)
		{
				(*derivs)(x,y,dydx);
				for(i=1; i<=nvar; i++)
				{
						yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
				}
				if(kmax>0 && kount<kmax-1 && fabs(x-xsav)>fabs(dxsav))
				{
						xp[++kount]=x;
						for(i=1; i<=nvar; i++)
						{
								yp[i][kount]=y[i];
						}
						xsav=x;
				}
				if((x+h-x2)*(x+h-x1)>0.0)
				{
						h=x2-x;
				}
				(*rkqs)(y, dydx, nvar, &x, h, eps, yscal, &hdid, &hnext, derivs);
				if(hdid==h)
				{
						++(*nok);
				}
				else
				{
						++(*nbad);
				}
				if((x-x2)*(x2-x1)>=0.0)
				{
						for(i=1; i<=nvar; i++)
						{
								ystart[i]=y[i];
						}
						if(kmax)
						{
								xp[++kount]=x;
								for(i=1; i<=nvar; i++)
								{
										yp[i][kount]=y[i];
								}
						}
						free_vector(dydx,1,nvar);
						free_vector(y,1,nvar);
						free_vector(y,1,nvar);
						return;
				}
				if(fabs(hnext)<=hmin) nrerror("Step size too small in odeint");
				h=hnext;
		}
		nrerror("Too many steps in routine odeint");
}

#undef MAXSTP
#undef TINY
#undef NRANSI
bsstep

Codice: Seleziona tutto

#include <math.h>

#define IMAX 11
#define NUSE 7
#define SHRINK 0.95
#define GROW 1.2

float **d=0,*x=0;	/* defining declaration */

void bsstep(y,dydx,nv,xx,htry,eps,yscal,hdid,hnext,derivs)
float y[],dydx[],*xx,htry,eps,yscal[],*hdid,*hnext;
void (*derivs)();
int nv;
{
	int i,j;
	double xsav,xest,h,errmax,temp;
	double *ysav,*dysav,*yseq,*yerr,*vector(),**matrix();
	static int nseq[IMAX+1]={0,2,4,6,8,12,16,24,32,48,64,96};
	void mmid(),rzextr(),nrerror(),free_matrix(),free_vector();

	ysav=vector(1,nv);
	dysav=vector(1,nv);
	yseq=vector(1,nv);
	yerr=vector(1,nv);
	x=vector(1,IMAX);
	d=matrix(1,nv,1,NUSE);
	h=htry;
	xsav=(*xx);
	for (i=1;i<=nv;i++) {
		ysav[i]=y[i];
		dysav[i]=dydx[i];
	}
	for (;;) {
		for (i=1;i<=IMAX;i++) {
			mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq,derivs);
			xest=(temp=h/nseq[i],temp*temp);
			rzextr(i,xest,yseq,y,yerr,nv,NUSE);
			errmax=0.0;
			for (j=1;j<=nv;j++)
				if (errmax < fabs(yerr[j]/yscal[j]))
					errmax=fabs(yerr[j]/yscal[j]);
			errmax /= eps;
			if (errmax < 1.0) {
				*xx += h;
				*hdid=h;
				*hnext = i==NUSE? h*SHRINK : i==NUSE-1?
					h*GROW : (h*nseq[NUSE-1])/nseq[i];
				free_matrix(d,1,nv,1,NUSE);
				free_vector(x,1,IMAX);
				free_vector(yerr,1,nv);
				free_vector(yseq,1,nv);
				free_vector(dysav,1,nv);
				free_vector(ysav,1,nv);
				return;
			}
		}
		h *= 0.25;
		for (i=1;i<=(IMAX-NUSE)/2;i++) h /= 2.0;
		if ((*xx+h) == (*xx)) nrerror("Step size underflow in BSSTEP");
	}
}

#undef IMAX
#undef NUSE
#undef SHRINK
#undef GROW
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

La derivata ti serve?
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Per ora no, infatti non ho tradotto hypdrv. L'ho lasciata stare in hypster, perche' non penso dia fastidio, e se dovesse servirmi poi c'e' gia'.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Mi da' molti errori, per lo piu' di dichiarazione:

Codice: Seleziona tutto

hypgeo.c:6:12: warning: missing whitespace after the macro name
hypgeo.c: In function ‘hypgeo’:
hypgeo.c:20: error: expected declaration specifiers or ‘...’ before ‘nok’
hypgeo.c:21: error: expected declaration specifiers or ‘...’ before ‘y’
hypgeo.c:24: error: expected declaration specifiers or ‘...’ before ‘kmax’
hypgeo.c:25: error: expected declaration specifiers or ‘...’ before ‘if’
hypgeo.c:43: error: expected declaration specifiers or ‘...’ before ‘bb’
hypgeo.c:44: error: expected declaration specifiers or ‘...’ before ‘cc’
hypgeo.c:45: error: expected declaration specifiers or ‘...’ before ‘dz’
hypgeo.c:46: error: expected declaration specifiers or ‘...’ before ‘hypser’
hypgeo.c:47: error: expected declaration specifiers or ‘...’ before ‘yy’
hypgeo.c:48: error: expected declaration specifiers or ‘...’ before ‘yy’
hypgeo.c:49: error: expected declaration specifiers or ‘...’ before ‘yy’
hypgeo.c:50: error: expected declaration specifiers or ‘...’ before ‘yy’
hypgeo.c:51: error: expected declaration specifiers or ‘...’ before ‘yy’
hypgeo.c:52: error: expected declaration specifiers or ‘...’ before ‘odeint’
hypgeo.c:53: error: expected declaration specifiers or ‘...’ before ‘y’
hypgeo.c:54: error: expected declaration specifiers or ‘...’ before ‘free_vector’
hypgeo.c:55: error: expected declaration specifiers or ‘...’ before ‘return’
hypgeo.c:56: error: expected declaration specifiers or ‘...’ before ‘}’ token
hypgeo.c:56: error: expected ‘;’, ‘,’ or ‘)’ before ‘}’ token
hypgeo.c:56: error: expected declaration or statement at end of input
Alcuni di questi ho pensato siano dovuti a problemi con gli header, anche se non tutti. Assumendo che hypgeo.c sia il main file, mi potresti fare cosa dovrebbe contenere un header, per esempio hypser.h? Grazie.
Avatar utente
bite
Imperturbabile Insigne
Imperturbabile Insigne
Messaggi: 3798
Iscrizione: sabato 19 maggio 2007, 22:10

Re: Funzione ipergeometrica di Gauss in C

Messaggio da bite »

Ho visto che la derivata è necessaria a odeint. Inoltre mancano alcune altre funzioni (rkck, rkqs, mmid, rzextr) richiamate direttamente o indirettamente da odeint. Se vuoi nel fine settimana sistemo il tutto e poi ti posto un piccolo tgz.
giuseppe morello
Prode Principiante
Messaggi: 68
Iscrizione: giovedì 5 aprile 2012, 9:56

Re: Funzione ipergeometrica di Gauss in C

Messaggio da giuseppe morello »

Se puoi farlo mi saresti davvero di grande aiuto. Io ovviamente ci sbattero' pure la testa, ma questo codice sta risultando molto piu' difficile di quelli che ho scritto finora e piu' difficile del previsto. Grazie.
Scrivi risposta

Ritorna a “Programmazione”

Chi c’è in linea

Visualizzano questa sezione: DoctorStrange e 10 ospiti