/****************************************************************** !* Purpose: This program computes the confluent * !* hypergeometric function M(a,b,x) using * !* subroutine CHGM * !* Input : a --- Parameter * !* b --- Parameter ( b <> 0,-1,-2,... ) * !* x --- Argument * !* Output: HG --- M(a,b,x) * !* Example: * !* a b x M(a,b,x) * !* ----------------------------------------- * !* 1.5 2.0 20.0 .1208527185D+09 * !* 4.5 2.0 20.0 .1103561117D+12 * !* -1.5 2.0 20.0 .1004836854D+05 * !* -4.5 2.0 20.0 -.3936045244D+03 * !* 1.5 2.0 50.0 .8231906643D+21 * !* 4.5 2.0 50.0 .9310512715D+25 * !* -1.5 2.0 50.0 .2998660728D+16 * !* -4.5 2.0 50.0 -.1807604439D+13 * !* -------------------------------------------------------------- * !* 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 void GAMMA(double X, double *GA) { /* ================================================== ! Purpose: Compute gamma function â(x) ! Input : x --- Argument of â(x) ! ( x is not equal to 0,-1,-2,úúú) ! Output: GA --- â(x) ! ================================================== */ double G[27]; double GR, PI, R, Z; int K, M, M1; PI=3.141592653589793; if (X == int(X)) if (X > 0.0) { *GA=1.0; M1=X-1; for (K=2; K<=M1; K++) *GA *= K; } else *GA=1.0e+300; else { if (fabs(X) > 1.0) { Z=fabs(X); M=int(Z); R=1.0; for (K=1; K<=M; K++) R=R*(Z-K); Z=Z-M; } else Z=X; G[1] = 1.0; G[2] = 0.5772156649015329; G[3] = -0.6558780715202538; G[4] = -0.420026350340952e-1; G[5] = 0.1665386113822915; G[6] = -0.421977345555443e-1; G[7] = -0.96219715278770e-2; G[8] = 0.72189432466630e-2; G[9] = -0.11651675918591e-2; G[10] = -0.2152416741149e-3; G[11] = 0.1280502823882e-3; G[12] = -0.201348547807e-4; G[13] = -0.12504934821e-5; G[14] = 0.11330272320e-5; G[15] = -0.2056338417e-6; G[16] = 0.61160950e-8; G[17] = 0.50020075e-8; G[18] = -0.11812746e-8; G[19] = 0.1043427e-9; G[20] = 0.77823e-11; G[21] = -0.36968e-11; G[22] = 0.51e-12; G[23] = -0.206e-13; G[24] = -0.54e-14; G[25] = 14e-14; G[26] = 1e-15; GR=G[26]; for (K=25; K>0; K--) GR=GR*Z+G[K]; *GA=1.0/(GR*Z); if (fabs(X) > 1.0) { *GA=(*GA)*R; if (X < 0.0) *GA=-PI/(X*(*GA)*sin(PI*X)); } } } void CHGM(double A, double B, double X, double *HG) { /* =================================================== ! Purpose: Compute confluent hypergeometric function ! M(a,b,x) ! Input : a --- Parameter ! b --- Parameter ( b <> 0,-1,-2,... ) ! x --- Argument ! Output: HG --- M(a,b,x) ! Routine called: GAMMA for computing â(x) ! =================================================== */ double A0,A1,HG1,HG2,PI,R,R1,R2,RG,SUM1,SUM2,TA,TB,TBA,X0,XG,Y0,Y1; int I, J, K, LA, M, N, NL; PI=3.141592653589793; A0=A; A1=A; X0=X; *HG=0.0; if (B == 0.0 || B == -fabs(int(B))) *HG=1.0e+300; else if(A == 0.0 || X == 0.0) *HG=1.0; else if(A == -1.0) *HG=1.0-X/B; else if(A == B) *HG=exp(X); else if(A-B == 1.0) *HG=(1.0+X/B)*exp(X); else if(A == 1.0 && B == 2.0) *HG=(exp(X)-1.0)/X; else if(A == int(A) && A < 0.0) { M=int(-A); R=1.0; *HG=1.0; for (K=1; K<=M; K++) { R=R*(A+K-1.0)/K/(B+K-1.0)*X; *HG=(*HG)+R; } } if (*HG != 0.0) return; if (X < 0.0) { A=B-A; A0=A; X=fabs(X); } if (A < 2.0) NL=0; if (A >= 2.0) { NL=1; LA=int(A); A=A-LA-1.0; } for (N=0; N<=NL; N++) { if (A0 >= 2.0) A=A+1.0; if (X <= 30.0+fabs(B) || A < 0.0) { *HG=1.0; RG=1.0; for (J=1; J<501; J++) { RG=RG*(A+J-1.0)/(J*(B+J-1.0))*X; *HG=(*HG)+RG; if (fabs(RG/(*HG)) < 1.0e-15) goto e25; } } else { GAMMA(A,&TA); GAMMA(B,&TB); XG=B-A; GAMMA(XG,&TBA); SUM1=1.0; SUM2=1.0; R1=1.0; R2=1.0; for (I=1; I<9; I++) { R1=-R1*(A+I-1.0)*(A-B+I)/(X*I); R2=-R2*(B-A+I-1.0)*(A-I)/(X*I); SUM1=SUM1+R1; SUM2=SUM2+R2; } HG1=TB/TBA*pow(X,-A)*cos(PI*A)*SUM1; HG2=TB/TA*exp(X)*pow(X,A-B)*SUM2; *HG=HG1+HG2; } e25: if (N == 0) Y0=*HG; if (N == 1) Y1=*HG; } if (A0 >= 2.0) for (I=1; I