/********************************************************************
*         Calculate the Psi function for a complex argument         *
* ----------------------------------------------------------------- *
* EXPLANATION:                                                      *
*                                                                   *
*      Purpose: This program computes the psi function psi(z)       *
*               for a complex argument using subroutine CPSI        *
*      Input :  x   --- Real part of z                              *
*               y   --- Imaginary part of z                         *
*      Output:  PSR --- Real part of psi(z)                         *
*               PSI --- Imaginary part of psi(z)                    *
*      Examples:                                                    *
*                  x       y      Re[psi(z)]     Im[psi(z)]         *
*                -------------------------------------------        *
*                 3.0     2.0     1.16459152      .67080728         *
*                 3.0    -2.0     1.16459152     -.67080728         *
*                -3.0     2.0     1.39536075     2.62465344         *
*                -3.0    -2.0     1.39536075    -2.62465344         *
*                                                                   *
* NOTE:          Psi(z) = Gamma'(z) / Gamma(z)                      *
* ----------------------------------------------------------------- *
* SAMPLE RUN:                                                       *
*                                                                   *
*  Please enter x and y (z=x+iy): 3 2                               *
*                                                                   *
*     x       y      Re[Psi(z)]      Im[Psi(z)]                     *
*   ----------------------------------------------                  *
*    3.0     2.0     1.16459152      0.67080728                     *
*                                                                   *
* ----------------------------------------------------------------- *
* REFERENCE:"Fortran Routines for Computation of Special Functions, *
*            jin.ece.uiuc.edu/routines/routines.html".              *
*                                                                   *
*                                 C++ Release By J-P Moreau, Paris. *
********************************************************************/
#include <stdio.h>
#include <math.h>

  double aa,E,X,Y,Z,PSR,PSI;
  double A[10], W[10];  //index 0 not used


// *** Begin utility functions for TANH
void WCoeffs(double X) {
  int I;
  aa=0.5;
  Z=X;
  for (I=1; I<10; I++) {
    W[I]=0.0;
    if (Z>aa)  W[I]=1.0;
    Z -= W[I]*aa;
    aa /= 2.0;
  }
}

void Coeffs() {
  E=exp(1.0);
  A[1]=1.648721270700128;
  A[2]=1.284025416687742;
  A[3]=1.133148453066826;
  A[4]=1.064494458917859;
  A[5]=1.031743407499103;
  A[6]=1.015747708586686;
  A[7]=1.007843097206448;
  A[8]=1.003913889338348;
  A[9]=1.001955033591003;
}

// Modified cordic exponential subroutine
double XExp(double X) {
// This subroutine takes an input value and returns Y=EXP(X)
// X may be any positive or negative real value.
  double Y;
  int I,K;
// Get coefficients
  Coeffs();
// Reduce the range of X
  K=(int) floor(X);
  X=X-K;
// Determine the weighting coeffs, W(I)
  WCoeffs(X);
// Calculate products
  Y=1.0;
  for (I=1; I<10; I++)
    if (W[I]>0.0)  Y *= A[I];
// Perform residual multiplication
  Y *= (1.0+Z*(1.0+Z/2.0*(1.0+Z/3.0*(1.0+Z/4.0))));
// Account for factor EXP(K)
  if (K<0)  E=1.0/E;
  if (abs(K)<1) return Y;
  for (I=1; I<=abs(K); I++) Y *= E;
// Restore X
  X += K;
  return Y;
}

/*--------------------------------------------*
*          Hyperbolic sine Function           *
* ------------------------------------------- *
* This Procedure uses the definition of the   *
* hyperbolic sine and the modified cordic     *
* exponential Function XExp(X) to approximate *
* SINH(X) over the entire range of real X.    *
* -------------------------------------------*/
double SinH(double X) {
  double Y;
  int I;
// Is X small enough to cause round off error?
  if (fabs(X)<0.35) goto e10;
// Calculate SinH(X) using exponential definition
// Get Y=exp(X)
  Y=XExp(X);
  Y=(Y-(1.0/Y))/2.0;
  return Y;
e10:; //series approximation (for X small)
  Z=1.0; Y=1.0;
  for (I=1; I<9; I++) {
    Z *= X*X/((2*I)*(2*I+1));
    Y += Z;
  }
  Y *= X;
  return Y;
}

// hyperbolic cosine Function
double CosH(double X) {
  double Y;
  Y=XExp(X);
  Y=(Y+(1.0/Y))/2.0;
  return Y;
}

// *** End utility functions for TANH ***

// hyperbolic tangent Function
// TanH(X]=SinH(X)/CosH(X)
double TanH(double X) {
  double V,Y;
  V=SinH(X);
  Y=CosH(X);
  return V/Y;
}

      double Sqr(double x) {
        return x*x;
      }

      void CPSI(double X, double Y, double *PSR, double *PSI) {
/*    =============================================
        Purpose: Compute the psi function for a
                 complex argument
        Input :  x   --- Real part of z
                 y   --- Imaginary part of z
        Output:  PSR --- Real part of psi(z)
                 PSI --- Imaginary part of psi(z)
      ============================================= */
        double A[9]; //index 0 not used
        double CT2,PI,RI,RR,TH,TM,TN,X0,X1,Y1,Z0,Z2;
        int    K,N;

		// Initialize Table A
        A[1]=-0.8333333333333e-01;     A[2]=0.83333333333333333e-02;
        A[3]=-0.39682539682539683e-02; A[4]=0.41666666666666667e-02;
        A[5]=-0.75757575757575758e-02; A[6]=0.21092796092796093e-01;
        A[7]=-0.83333333333333333e-01; A[8]=0.4432598039215686;
        
        PI=4.0*atan(1.0);

        if (Y == 0.0 && X == int(X) && X <= 0.0) {
           *PSR=1e+200;
           *PSI=0.0;
        }
        else {
		   if (X < 0.0) {
             X1=X;
             Y1=Y;
             X=-X;
             Y=-Y;
           }
           X0=X;
           if (X < 8.0) {
             N=8-int(X);
             X0=X+1.0*N;
           }
           if (X0 == 0.0 && Y != 0.0)  TH=0.5*PI;
           if (X0 != 0.0)  TH=atan(Y/X0);
           Z2=X0*X0+Y*Y;
           Z0=sqrt(Z2);
           *PSR=log(Z0)-0.5*X0/Z2;
           *PSI=TH+0.5*Y/Z2;
           for (K=1; K<9; K++) {
             *PSR=*PSR+A[K]*pow(Z2,-K)*cos(2.0*K*TH);
             *PSI=*PSI-A[K]*pow(Z2,-K)*sin(2.0*K*TH);
           }
           if (X < 8.0) {
              RR=0.0;
              RI=0.0;
              for (K=1; K<=N; K++) {
                RR=RR+(X0-K)/(Sqr(X0-K)+Y*Y);
                RI=RI+Y/(Sqr(X0-K)+Y*Y);
              }
              *PSR=*PSR-RR;
              *PSI=*PSI+RI;
           }
           if (X1 < 0.0) {
              TN=sin(PI*X)/cos(PI*X);
              TM=TanH(PI*Y);
              CT2=TN*TN+TM*TM;
              *PSR=*PSR+X/(X*X+Y*Y)+PI*(TN-TN*TM*TM)/CT2;
              *PSI=*PSI-Y/(X*X+Y*Y)-PI*TM*(1.0+TN*TN)/CT2;
              X=X1;
              Y=Y1;
           }
        }
      } //CPSI()


  void main()  {

    printf("\n Please enter x and y (z=x+iy): "); scanf("%lf %lf", &X, &Y);
    printf("\n");
    printf("   x       y      Re[Psi(z)]      Im[Psi(z)]\n");
    printf(" ----------------------------------------------\n");

    CPSI(X,Y,&PSR,&PSI);

    printf("%5.1f   %5.1f  %13.8f   %13.8f\n\n", X, Y, PSR, PSI);

  }

// end of file mcpsi.cpp