FUNCTION PSI(XX)
C----------------------------------------------------------------------
C
C This function program evaluates the logarithmic derivative of the
C   gamma function, 
C
C      psi(x) = d/dx (gamma(x)) / gamma(x) = d/dx (ln gamma(x))
C
C   for real x, where either
C
C          -xmax1 < x < -xmin (x not a negative integer), or
C            xmin < x.
C
C   The calling sequence for this function is 
C
C                  Y = PSI(X)
C
C   The main computation uses rational Chebyshev approximations
C   published in Math. Comp. 27, 123-127 (1973) by Cody, Strecok and
C   Thacher.  This transportable program is patterned after the
C   machine-dependent FUNPACK program PSI(X), but cannot match that
C   version for efficiency or accuracy.  This version uses rational
C   approximations that are theoretically accurate to 20 significant
C   decimal digits.  The accuracy achieved depends on the arithmetic
C   system, the compiler, the intrinsic functions, and proper selection
C   of the machine-dependent constants.
C
C*******************************************************************
C*******************************************************************
C
C Explanation of machine-dependent constants
C
C   XINF   = largest positive machine number
C   XMAX1  = beta ** (p-1), where beta is the radix for the
C            floating-point system, and p is the number of base-beta
C            digits in the floating-point significand.  This is an
C            upper bound on non-integral floating-point numbers, and
C            the negative of the lower bound on acceptable negative
C            arguments for PSI.  If rounding is necessary, round this
C            value down.
C   XMIN1  = the smallest in magnitude acceptable argument.  We
C            recommend XMIN1 = MAX(1/XINF,xmin) rounded up, where
C            xmin is the smallest positive floating-point number.
C   XSMALL = absolute argument below which  PI*COTAN(PI*X)  may be
C            represented by 1/X.  We recommend XSMALL < sqrt(3 eps)/pi,
C            where eps is the smallest positive number such that
C            1+eps > 1. 
C   XLARGE = argument beyond which PSI(X) may be represented by
C            LOG(X).  The solution to the equation
C               x*ln(x) = beta ** p
C            is a safe value.
C
C     Approximate values for some important machines are
C
C                        beta  p     eps     xmin       XINF  
C
C  CDC 7600      (S.P.)    2  48  7.11E-15  3.13E-294  1.26E+322
C  CRAY-1        (S.P.)    2  48  7.11E-15  4.58E-2467 5.45E+2465
C  IEEE (IBM/XT,
C    SUN, etc.)  (S.P.)    2  24  1.19E-07  1.18E-38   3.40E+38
C  IEEE (IBM/XT,
C    SUN, etc.)  (D.P.)    2  53  1.11D-16  2.23E-308  1.79D+308
C  IBM 3033      (D.P.)   16  14  1.11D-16  5.40D-79   7.23D+75
C  SUN 3/160     (D.P.)    2  53  1.11D-16  2.23D-308  1.79D+308
C  VAX 11/780    (S.P.)    2  24  5.96E-08  2.94E-39   1.70E+38
C                (D.P.)    2  56  1.39D-17  2.94D-39   1.70D+38
C   (G Format)   (D.P.)    2  53  1.11D-16  5.57D-309  8.98D+307
C
C                         XMIN1      XMAX1     XSMALL    XLARGE
C
C  CDC 7600      (S.P.)  3.13E-294  1.40E+14  4.64E-08  9.42E+12
C  CRAY-1        (S.P.)  1.84E-2466 1.40E+14  4.64E-08  9.42E+12
C  IEEE (IBM/XT,
C    SUN, etc.)  (S.P.)  1.18E-38   8.38E+06  1.90E-04  1.20E+06
C  IEEE (IBM/XT,
C    SUN, etc.)  (D.P.)  2.23D-308  4.50D+15  5.80D-09  2.71D+14
C  IBM 3033      (D.P.)  1.39D-76   4.50D+15  5.80D-09  2.05D+15
C  SUN 3/160     (D.P.)  2.23D-308  4.50D+15  5.80D-09  2.71D+14
C  VAX 11/780    (S.P.)  5.89E-39   8.38E+06  1.35E-04  1.20E+06
C                (D.P.)  5.89D-39   3.60D+16  2.05D-09  2.05D+15
C   (G Format)   (D.P.)  1.12D-308  4.50D+15  5.80D-09  2.71D+14
C
C*******************************************************************
C*******************************************************************
C
C Error Returns
C
C  The program returns XINF for  X < -XMAX1, for X zero or a negative
C    integer, or when X lies in (-XMIN1, 0), and returns -XINF
C    when X lies in (0, XMIN1).
C
C Intrinsic functions required are:
C
C     ABS, AINT, DBLE, INT, LOG, REAL, TAN
C
C
C  Author: W. J. Cody
C          Mathematics and Computer Science Division 
C          Argonne National Laboratory
C          Argonne, IL 60439
C
C  Latest modification: June 8, 1988
C
C----------------------------------------------------------------------
      INTEGER I,N,NQ
CS    REAL
CD    DOUBLE PRECISION
     1   AUG,CONV,DEN,PSI,FOUR,FOURTH,HALF,ONE,P1,P2,PIOV4,Q1,Q2,
     2   SGN,THREE,XLARGE,UPPER,W,X,XINF,XMAX1,XMIN1,XSMALL,X01,
     3   X01D,X02,XX,Z,ZERO
      DIMENSION P1(9),P2(7),Q1(8),Q2(6)
C----------------------------------------------------------------------
C  Mathematical constants.  PIOV4 = pi / 4
C----------------------------------------------------------------------
CS    DATA ZERO,FOURTH,HALF,ONE/0.0E0,0.25E0,0.5E0,1.0E0/
CS    DATA THREE,FOUR/3.0E0,4.0E0/,PIOV4/7.8539816339744830962E-01/
CD    DATA ZERO,FOURTH,HALF,ONE/0.0D0,0.25D0,0.5D0,1.0D0/
CD    DATA THREE,FOUR/3.0D0,4.0D0/,PIOV4/7.8539816339744830962D-01/
C----------------------------------------------------------------------
C  Machine-dependent constants
C----------------------------------------------------------------------
CS    DATA XINF/1.70E+38/, XMIN1/5.89E-39/, XMAX1/8.38E+06/,
CS   1     XSMALL/1.35E-04/, XLARGE/1.20E+06/
CD    DATA XINF/1.70D+38/, XMIN1/5.89D-39/, XMAX1/3.60D+16/,
CD   1     XSMALL/2.05D-09/, XLARGE/2.04D+15/
C----------------------------------------------------------------------
C  Zero of psi(x)
C----------------------------------------------------------------------
CS    DATA X01/187.0E0/,X01D/128.0E0/,X02/6.9464496836234126266E-04/
CD    DATA X01/187.0D0/,X01D/128.0D0/,X02/6.9464496836234126266D-04/
C----------------------------------------------------------------------
C  Coefficients for approximation to  psi(x)/(x-x0)  over [0.5, 3.0]
C----------------------------------------------------------------------
CS    DATA P1/4.5104681245762934160E-03,5.4932855833000385356E+00,
CS   1        3.7646693175929276856E+02,7.9525490849151998065E+03,
CS   2        7.1451595818951933210E+04,3.0655976301987365674E+05,
CS   3        6.3606997788964458797E+05,5.8041312783537569993E+05,
CS   4        1.6585695029761022321E+05/
CS    DATA Q1/9.6141654774222358525E+01,2.6287715790581193330E+03,
CS   1        2.9862497022250277920E+04,1.6206566091533671639E+05,
CS   2        4.3487880712768329037E+05,5.4256384537269993733E+05,
CS   3        2.4242185002017985252E+05,6.4155223783576225996E-08/
CD    DATA P1/4.5104681245762934160D-03,5.4932855833000385356D+00,
CD   1        3.7646693175929276856D+02,7.9525490849151998065D+03,
CD   2        7.1451595818951933210D+04,3.0655976301987365674D+05,
CD   3        6.3606997788964458797D+05,5.8041312783537569993D+05,
CD   4        1.6585695029761022321D+05/
CD    DATA Q1/9.6141654774222358525D+01,2.6287715790581193330D+03,
CD   1        2.9862497022250277920D+04,1.6206566091533671639D+05,
CD   2        4.3487880712768329037D+05,5.4256384537269993733D+05,
CD   3        2.4242185002017985252D+05,6.4155223783576225996D-08/
C----------------------------------------------------------------------
C  Coefficients for approximation to  psi(x) - ln(x) + 1/(2x) 
C     for  x > 3.0
C----------------------------------------------------------------------
CS    DATA P2/-2.7103228277757834192E+00,-1.5166271776896121383E+01,
CS   1        -1.9784554148719218667E+01,-8.8100958828312219821E+00,
CS   2        -1.4479614616899842986E+00,-7.3689600332394549911E-02,
CS   3        -6.5135387732718171306E-21/
CS    DATA Q2/ 4.4992760373789365846E+01, 2.0240955312679931159E+02,
CS   1         2.4736979003315290057E+02, 1.0742543875702278326E+02,
CS   2         1.7463965060678569906E+01, 8.8427520398873480342E-01/
CD    DATA P2/-2.7103228277757834192D+00,-1.5166271776896121383D+01,
CD   1        -1.9784554148719218667D+01,-8.8100958828312219821D+00,
CD   2        -1.4479614616899842986D+00,-7.3689600332394549911D-02,
CD   3        -6.5135387732718171306D-21/
CD    DATA Q2/ 4.4992760373789365846D+01, 2.0240955312679931159D+02,
CD   1         2.4736979003315290057D+02, 1.0742543875702278326D+02,
CD   2         1.7463965060678569906D+01, 8.8427520398873480342D-01/
C----------------------------------------------------------------------
CS    CONV(I) = REAL(I)
CD    CONV(I) = DBLE(I)
      X = XX
      W = ABS(X)
      AUG = ZERO
C----------------------------------------------------------------------
C  Check for valid arguments, then branch to appropriate algorithm
C----------------------------------------------------------------------
      IF ((-X .GE. XMAX1) .OR. (W .LT. XMIN1)) THEN
            GO TO 410
         ELSE IF (X .GE. HALF) THEN
            GO TO 200
C----------------------------------------------------------------------
C  X < 0.5, use reflection formula: psi(1-x) = psi(x) + pi * cot(pi*x)
C     Use 1/X for PI*COTAN(PI*X)  when  XMIN1 < |X| <= XSMALL.  
C----------------------------------------------------------------------
         ELSE IF (W .LE. XSMALL) THEN
            AUG = -ONE / X
            GO TO 150
      END IF
C----------------------------------------------------------------------
C  Argument reduction for cot
C----------------------------------------------------------------------
  100 IF (X .LT. ZERO) THEN
            SGN = PIOV4
         ELSE
            SGN = -PIOV4
      END IF
      W = W - AINT(W)
      NQ = INT(W * FOUR)
      W = FOUR * (W - CONV(NQ) * FOURTH)
C----------------------------------------------------------------------
C  W is now related to the fractional part of  4.0 * X.
C     Adjust argument to correspond to values in the first
C     quadrant and determine the sign.
C----------------------------------------------------------------------
      N = NQ / 2
      IF ((N+N) .NE. NQ) W = ONE - W
      Z = PIOV4 * W
      IF (MOD(N,2) .NE. 0) SGN = - SGN
C----------------------------------------------------------------------
C  determine the final value for  -pi * cotan(pi*x)
C----------------------------------------------------------------------
      N = (NQ + 1) / 2
      IF (MOD(N,2) .EQ. 0) THEN
C----------------------------------------------------------------------
C  Check for singularity
C----------------------------------------------------------------------
            IF (Z .EQ. ZERO) GO TO 410
            AUG = SGN * (FOUR / TAN(Z))
         ELSE
            AUG = SGN * (FOUR * TAN(Z))
      END IF
  150 X = ONE - X
  200 IF (X .GT. THREE) GO TO 300
C----------------------------------------------------------------------
C  0.5 <= X <= 3.0
C----------------------------------------------------------------------
      DEN = X
      UPPER = P1(1) * X
      DO 210 I = 1, 7
         DEN = (DEN + Q1(I)) * X
         UPPER = (UPPER + P1(I+1)) * X
  210 CONTINUE
      DEN = (UPPER + P1(9)) / (DEN + Q1(8))
      X = (X-X01/X01D) - X02
      PSI = DEN * X + AUG
      GO TO 500
C----------------------------------------------------------------------
C  3.0 < X 
C----------------------------------------------------------------------
  300 IF (X .LT. XLARGE) THEN
         W = ONE / (X * X)
         DEN = W
         UPPER = P2(1) * W
         DO 310 I = 1, 5
            DEN = (DEN + Q2(I)) * W
            UPPER = (UPPER + P2(I+1)) * W
  310    CONTINUE
         AUG = (UPPER + P2(7)) / (DEN + Q2(6)) - HALF / X + AUG
      END IF
      PSI = AUG + LOG(X)
      GO TO 500
C----------------------------------------------------------------------
C  Error return
C----------------------------------------------------------------------
  410 PSI = XINF
      IF (X .GT. ZERO) PSI = -XINF
  500 RETURN
C---------- Last card of PSI ----------
      END