00001 #include "cdflib.h"
00002 double psi(double *xx)
00003 /*
00004 ---------------------------------------------------------------------
00005  
00006                  EVALUATION OF THE DIGAMMA FUNCTION
00007  
00008                            -----------
00009  
00010      PSI(XX) IS ASSIGNED THE VALUE 0 WHEN THE DIGAMMA FUNCTION CANNOT
00011      BE COMPUTED.
00012  
00013      THE MAIN COMPUTATION INVOLVES EVALUATION OF RATIONAL CHEBYSHEV
00014      APPROXIMATIONS PUBLISHED IN MATH. COMP. 27, 123-127(1973) BY
00015      CODY, STRECOK AND THACHER.
00016  
00017 ---------------------------------------------------------------------
00018      PSI WAS WRITTEN AT ARGONNE NATIONAL LABORATORY FOR THE FUNPACK
00019      PACKAGE OF SPECIAL FUNCTION SUBROUTINES. PSI WAS MODIFIED BY
00020      A.H. MORRIS (NSWC).
00021 ---------------------------------------------------------------------
00022 */
00023 {
00024 static double dx0 = 1.461632144968362341262659542325721325e0;
00025 static double piov4 = .785398163397448e0;
00026 static double p1[7] = {
00027     .895385022981970e-02,.477762828042627e+01,.142441585084029e+03,
00028     .118645200713425e+04,.363351846806499e+04,.413810161269013e+04,
00029     .130560269827897e+04
00030 };
00031 static double p2[4] = {
00032     -.212940445131011e+01,-.701677227766759e+01,-.448616543918019e+01,
00033     -.648157123766197e+00
00034 };
00035 static double q1[6] = {
00036     .448452573429826e+02,.520752771467162e+03,.221000799247830e+04,
00037     .364127349079381e+04,.190831076596300e+04,.691091682714533e-05
00038 };
00039 static double q2[4] = {
00040     .322703493791143e+02,.892920700481861e+02,.546117738103215e+02,
00041     .777788548522962e+01
00042 };
00043 static int K1 = 3;
00044 static int K2 = 1;
00045 static double psi,aug,den,sgn,upper,w,x,xmax1,xmx0,xsmall,z;
00046 static int i,m,n,nq;
00047 /*
00048      ..
00049      .. Executable Statements ..
00050 */
00051 /*
00052 ---------------------------------------------------------------------
00053      MACHINE DEPENDENT CONSTANTS ...
00054         XMAX1  = THE SMALLEST POSITIVE FLOATING POINT CONSTANT
00055                  WITH ENTIRELY INTEGER REPRESENTATION.  ALSO USED
00056                  AS NEGATIVE OF LOWER BOUND ON ACCEPTABLE NEGATIVE
00057                  ARGUMENTS AND AS THE POSITIVE ARGUMENT BEYOND WHICH
00058                  PSI MAY BE REPRESENTED AS ALOG(X).
00059         XSMALL = ABSOLUTE ARGUMENT BELOW WHICH PI*COTAN(PI*X)
00060                  MAY BE REPRESENTED BY 1/X.
00061 ---------------------------------------------------------------------
00062 */
00063     xmax1 = ipmpar(&K1);
00064     xmax1 = fifdmin1(xmax1,1.0e0/spmpar(&K2));
00065     xsmall = 1.e-9;
00066     x = *xx;
00067     aug = 0.0e0;
00068     if(x >= 0.5e0) goto S50;
00069 /*
00070 ---------------------------------------------------------------------
00071      X .LT. 0.5,  USE REFLECTION FORMULA
00072      PSI(1-X) = PSI(X) + PI * COTAN(PI*X)
00073 ---------------------------------------------------------------------
00074 */
00075     if(fabs(x) > xsmall) goto S10;
00076     if(x == 0.0e0) goto S100;
00077 /*
00078 ---------------------------------------------------------------------
00079      0 .LT. ABS(X) .LE. XSMALL.  USE 1/X AS A SUBSTITUTE
00080      FOR  PI*COTAN(PI*X)
00081 ---------------------------------------------------------------------
00082 */
00083     aug = -(1.0e0/x);
00084     goto S40;
00085 S10:
00086 /*
00087 ---------------------------------------------------------------------
00088      REDUCTION OF ARGUMENT FOR COTAN
00089 ---------------------------------------------------------------------
00090 */
00091     w = -x;
00092     sgn = piov4;
00093     if(w > 0.0e0) goto S20;
00094     w = -w;
00095     sgn = -sgn;
00096 S20:
00097 /*
00098 ---------------------------------------------------------------------
00099      MAKE AN ERROR EXIT IF X .LE. -XMAX1
00100 ---------------------------------------------------------------------
00101 */
00102     if(w >= xmax1) goto S100;
00103     nq = fifidint(w);
00104     w -= (double)nq;
00105     nq = fifidint(w*4.0e0);
00106     w = 4.0e0*(w-(double)nq*.25e0);
00107 /*
00108 ---------------------------------------------------------------------
00109      W IS NOW RELATED TO THE FRACTIONAL PART OF  4.0 * X.
00110      ADJUST ARGUMENT TO CORRESPOND TO VALUES IN FIRST
00111      QUADRANT AND DETERMINE SIGN
00112 ---------------------------------------------------------------------
00113 */
00114     n = nq/2;
00115     if(n+n != nq) w = 1.0e0-w;
00116     z = piov4*w;
00117     m = n/2;
00118     if(m+m != n) sgn = -sgn;
00119 /*
00120 ---------------------------------------------------------------------
00121      DETERMINE FINAL VALUE FOR  -PI*COTAN(PI*X)
00122 ---------------------------------------------------------------------
00123 */
00124     n = (nq+1)/2;
00125     m = n/2;
00126     m += m;
00127     if(m != n) goto S30;
00128 /*
00129 ---------------------------------------------------------------------
00130      CHECK FOR SINGULARITY
00131 ---------------------------------------------------------------------
00132 */
00133     if(z == 0.0e0) goto S100;
00134 /*
00135 ---------------------------------------------------------------------
00136      USE COS/SIN AS A SUBSTITUTE FOR COTAN, AND
00137      SIN/COS AS A SUBSTITUTE FOR TAN
00138 ---------------------------------------------------------------------
00139 */
00140     aug = sgn*(cos(z)/sin(z)*4.0e0);
00141     goto S40;
00142 S30:
00143     aug = sgn*(sin(z)/cos(z)*4.0e0);
00144 S40:
00145     x = 1.0e0-x;
00146 S50:
00147     if(x > 3.0e0) goto S70;
00148 /*
00149 ---------------------------------------------------------------------
00150      0.5 .LE. X .LE. 3.0
00151 ---------------------------------------------------------------------
00152 */
00153     den = x;
00154     upper = p1[0]*x;
00155     for(i=1; i<=5; i++) {
00156         den = (den+q1[i-1])*x;
00157         upper = (upper+p1[i+1-1])*x;
00158     }
00159     den = (upper+p1[6])/(den+q1[5]);
00160     xmx0 = x-dx0;
00161     psi = den*xmx0+aug;
00162     return psi;
00163 S70:
00164 /*
00165 ---------------------------------------------------------------------
00166      IF X .GE. XMAX1, PSI = LN(X)
00167 ---------------------------------------------------------------------
00168 */
00169     if(x >= xmax1) goto S90;
00170 /*
00171 ---------------------------------------------------------------------
00172      3.0 .LT. X .LT. XMAX1
00173 ---------------------------------------------------------------------
00174 */
00175     w = 1.0e0/(x*x);
00176     den = w;
00177     upper = p2[0]*w;
00178     for(i=1; i<=3; i++) {
00179         den = (den+q2[i-1])*w;
00180         upper = (upper+p2[i+1-1])*w;
00181     }
00182     aug = upper/(den+q2[3])-0.5e0/x+aug;
00183 S90:
00184     psi = aug+log(x);
00185     return psi;
00186 S100:
00187 /*
00188 ---------------------------------------------------------------------
00189      ERROR RETURN
00190 ---------------------------------------------------------------------
00191 */
00192     psi = 0.0e0;
00193     return psi;
00194 } /* END */