Gamma function: implementation?
mccaugh at s.cs.uiuc.edu
mccaugh at s.cs.uiuc.edu
Thu Sep 21 13:40:00 AEST 1989
Brian:
I am posting this in case mailing fails (also in case anyone wants to correct):
#include <math.h>
/* Gamma - xx = argument, */
/* ier = return-code (0: no error; 1: negative; 2: overflow */
double Gamma (xx,ier)
double xx;
int *ier;
{
double x,y, Gx,Gy, err;
#define MAX_FLOAT 1.0E75
if(xx > 57.0) { *ier = 2; /* overflow */
return(MAX_FLOAT);
}
/* else xx <= 57.0 */
x = xx;
err = 1.0E-6; /* margin of error */
*ier = 0;
Gx = 1.0;
if(x > 2.0) { do { x = x - 1.0;
Gx = Gx*x;
}
while(x > 2.0);
/* x <= 2.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x <= 2.0 */
if(x == 1.0) return(Gx);
if(x > 1.0) { y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x < 1.0 */
if(x > err) { do { Gx = Gx/x;
x = x + 1.0;
}
while(x <= 1.0);
/* x > 1.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
/* else x <= err */
y = floor(x) - x;
if(fabs(y) <= err) { *ier = 1;
return(Gx);
}
/* else fabs(y) > err */
if(1-y <= err) { *ier = 1;
return(Gx);
}
/* else 1-y > err */
while(x <= 1.0) { Gx = Gx/x;
x = x + 1.0;
}
/* x > 1.0 */
y = x - 1.0;
Gy = 1.0 - 0.57710166*(y-1) + 0.98585399*(y-1)^2 -
0.87642182*(y-1)^3 + 0.83282120*(y-1)^4 -
0.56847290*(y-1)^5 + 0.25482049*(y-1)^6 -
0.05149930*(y-1)^7;
Gx = Gx*Gy;
return(Gx);
}
My home-brewed C takes exponentiation ('^') while other C's may not; other-
wise, this should work. Hope this helps, and good luck with your porting!
Scott McCaughrin
University of Illinois
Urbana, Illinois.
More information about the Comp.lang.c
mailing list