/******************************************************************
!*      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 <stdio.h>
#include <math.h>


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<LA; I++) {
              *HG=((2.0*A-B+X)*Y1+(B-A)*Y0)/A;
              Y0=Y1;
              Y1=*HG;
              A=A+1.0;
            }
        if (X0 < 0.0) *HG=(*HG)*exp(X0);
        A=A1;
        X=X0;
}
 

void main() {
        
        double A, B, X, HG;
	
        printf("\n");
        printf("   a       b      x          M(a,b,x)     \n");
        printf(" -----------------------------------------\n");

        A=1.5; B=2.0; X=20.0;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=4.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=-1.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=-4.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=1.5; X=50.0;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=4.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=-1.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n", A,B,X,HG);

        A=-4.5;
		CHGM(A,B,X,&HG);
        printf(" %5.1f  %5.1f  %5.1f    %e\n\n", A,B,X,HG);

}

// end of file mchgm.cpp