/************************************************************* * Calculate Function Gamma with a complex argument * * ---------------------------------------------------------- * * EXPLANATION: * * Purpose: This program computes the gamma function G(z) * * or Ln[G(z)] for a complex argument using * * subroutine CGAMA * * Input : x --- Real part of z * * y --- Imaginary part of z * * KF --- Function code * * KF=0 for Ln[G(z)] * * KF=1 for G(z) * * Output: GR --- Real part of Ln[G(z)] or G(z) * * GI --- Imaginary part of Ln[G(z)] or G(z) * * Examples: * * x y Re[G(z)] Im[G(z)] * * -------------------------------------------------------- * * 2.50 5.00 .2267360319E-01 -.1172284404E-01 * * 5.00 10.00 .1327696517E-01 .3639011746E-02 * * 2.50 -5.00 .2267360319E-01 .1172284404E-01 * * 5.00 -10.00 .1327696517E-01 -.3639011746E-02 * * * * x y Re[LnG(z)] Im[LnG(z)] * * -------------------------------------------------------- * * 2.50 5.00 -.3668103262E+01 .5806009801E+01 * * 5.00 10.00 -.4285507444E+01 .1911707090E+02 * * 2.50 -5.00 -.3668103262E+01 -.5806009801E+01 * * 5.00 -10.00 -.4285507444E+01 -.1911707090E+02 * * ---------------------------------------------------------- * * SAMPLE RUNS: * * * * Please enter KF, x and y: 1 2.5 5.0 * * * * x y Re[G(z)] Im[G(z)] * * -------------------------------------------------------- * * 2.50 5.00 0.0226736032 -0.0117228440 * * * * Please enter KF, zx and zy: 0 2.5 5.0 * * * * x y Re[LnG(z)] Im[LnG(z)] * * -------------------------------------------------------- * * 2.50 5.00 -3.6681032624 5.8060098006 * * * * ---------------------------------------------------------- * * REFERENCE: * * "Fortran Routines for Computation of Special Functions, * * jin.ece.uiuc.edu/routines/routines.html". * * * * C++ Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *************************************************************/ #include #include double X,Y,GR,GI; int KF; //auxiliary hyperbolic functions double SINH(double x) { double expx; expx = exp(x); return (0.5*(expx-1.0/expx)); } double COSH(double x) { double expx; expx = exp(x); return (0.5*(expx+1.0/expx)); } void CGAMA(double X, double Y, int KF, double *GR, double *GI) { /*=========================================================== Purpose: Compute the gamma function G(z) or Ln[G(z)] for a complex argument Input : x --- Real part of z y --- Imaginary part of z KF --- Function code KF=0 for Ln[G(z)] KF=1 for G(z) Output: GR --- Real part of Ln[G(z)] or G(z) GI --- Imaginary part of Ln[G(z)] or G(z) ===========================================================*/ double A[11]; //index 0 not used double G0,GR1,GI1,SR,SI,T,TH,TH1,TH2,X0,X1,Y1,Z1,Z2; double PI = 4.0*atan(1.0); int J,K,NA; //Initialize table A A[1]= 8.333333333333333e-02; A[2]= -2.777777777777778e-03; A[3]= 7.936507936507937e-04; A[4]= -5.952380952380952e-04; A[5]= 8.417508417508418e-04; A[6]= -1.917526917526918e-03; A[7]= 6.410256410256410e-03; A[8]= -2.955065359477124e-02; A[9]= 1.796443723688307e-01; A[10]=-1.39243221690590; if (Y == 0.0 && X == floor(X) && X <= 0.0) { *GR=1e+300; //arbitrary big number *GI=0.0; return; } else if (X < 0.0) { X1=X; Y1=Y; X=-X; Y=-Y; } X0=X; if (X <= 7.0) { NA=(int) floor(7.0-X); X0=X + 1.0*NA; } Z1=sqrt(X0*X0+Y*Y); TH=atan(Y/X0); *GR=(X0-0.5)*log(Z1)-TH*Y-X0+0.5*log(2.0*PI); *GI=TH*(X0-0.5)+Y*log(Z1)-Y; for (K=1; K<=10; K++) { T=pow(Z1,(1-2*K)); *GR=*GR+A[K]*T*cos((2.0*K-1.0)*TH); *GI=*GI-A[K]*T*sin((2.0*K-1.0)*TH); } if (X <= 7.0) { GR1=0.0; GI1=0.0; for (J=0; J