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

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<NA; J++) {
        GR1 += 0.5*log((X+J)*(X+J)+Y*Y);
        GI1 += atan(Y/(X+J));
      }
      *GR = *GR - GR1;
      *GI = *GI - GI1;
    }
    if (X1 < 0.0) {
      Z1=sqrt(X*X+Y*Y);
      TH1=atan(Y/X);
      SR=-sin(PI*X)*COSH(PI*Y);
      SI=-cos(PI*X)*SINH(PI*Y);
      Z2=sqrt(SR*SR+SI*SI);
      TH2=atan(SI/SR);
      if (SR < 0.0)  TH2 += PI;
      *GR=log(PI/(Z1*Z2))-*GR;
      *GI=-TH1-TH2-*GI;
      X=X1;
      Y=Y1;
    }
    if (KF == 1) {
      G0=exp(*GR);
      *GR=G0*cos(*GI);
      *GI=G0*sin(*GI);
    }
  }


  void main()  {

    printf("\n  Please enter KF, x and y: ");
    scanf("%d %lf %lf", &KF, &X, &Y);
    printf("\n");
    if (KF == 1)
      printf("        x         y           Re[G(z)]           Im[G(z)]\n");
    else
      printf("        x         y          Re[LnG(z)]          Im[LnG(z)]\n");
    printf("    ---------------------------------------------------------\n");

    CGAMA(X,Y,KF,&GR,&GI);

    printf(" %10.2f%10.2f%19.10f%19.10f\n\n", X, Y, GR, GI);

  }

// end of file mcgama.cpp