#include "fpreproc/library_names.h"
#include "fpreproc/fp_const.h"
*DECK DPOCH1
      REAL*8 FUNCTION DPOCH1 (A, X)
C***BEGIN PROLOGUE  DPOCH1
C***PURPOSE  Calculate a generalization of Pochhammer's symbol starting
C            from first order.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C1, C7A
C***TYPE      REAL*8 (POCH1-S, DPOCH1-D)
C***KEYWORDS  FIRST ORDER, FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate a double precision generalization of Pochhammer's symbol
C for double precision A and X for special situations that require
C especially accurate values when X is small in
C        POCH1(A,X) = (POCH(A,X)-1)/X
C                   = (GAMMA(A+X)/GAMMA(A) - 1.0)/X .
C This specification is particularly suited for stably computing
C expressions such as
C        (GAMMA(A+X)/GAMMA(A) - GAMMA(B+X)/GAMMA(B))/X
C             = POCH1(A,X) - POCH1(B,X)
C Note that POCH1(A,0.0) = PSI(A)
C
C When ABS(X) is so small that substantial cancellation will occur if
C the straightforward formula is used, we use an expansion due
C to Fields and discussed by Y. L. Luke, The Special Functions and Their
C Approximations, Vol. 1, Academic Press, 1969, page 34.
C
C The ratio POCH(A,X) = GAMMA(A+X)/GAMMA(A) is written by Luke as
C        (A+(X-1)/2)**X * polynomial in (A+(X-1)/2)**(-2) .
C In order to maintain significance in POCH1, we write for positive a
C        (A+(X-1)/2)**X = EXP(X*LOG(A+(X-1)/2)) = EXP(Q)
C                       = 1.0 + Q*EXPREL(Q) .
C Likewise the polynomial is written
C        POLY = 1.0 + X*POLY1(A,X) .
C Thus,
C        POCH1(A,X) = (POCH(A,X) - 1) / X
C                   = EXPREL(Q)*(Q/X + Q*POLY1(A,X)) + POLY1(A,X)
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCOT, DEXPRL, DPOCH, DPSI, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   890911  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900727  Added EXTERNAL statement.  (WRB)
C***END PROLOGUE  DPOCH1
C============
C idecl:  explicitize implicit INTEGER declarations:
      IMPLICIT NONE
      INTEGER incr,nterms,k,j,ndx,ii,i
C============
      REAL*8 A, X, ABSA, ABSX, ALNEPS, ALNVAR, B, BERN(20),
     1  BINV, BP, GBERN(21), GBK, PI, POLY1, Q, RHO, SINPXX, SINPX2,
     2  SQTBIG, TERM, TRIG, VAR, VAR2, D1MACH, DPSI, DEXPRL, DCOT, DPOCH
      LOGICAL FIRST
      EXTERNAL DCOT
      SAVE BERN, PI, SQTBIG, ALNEPS, FIRST
      DATA BERN  (  1) / +CONST(.833333333333333333333333333333333,-1)/
      DATA BERN  (  2) / -CONST(.138888888888888888888888888888888,-2)/
      DATA BERN  (  3) / +CONST(.330687830687830687830687830687830,-4)/
      DATA BERN  (  4) / -CONST(.826719576719576719576719576719576,-6)/
      DATA BERN  (  5) / +CONST(.208767569878680989792100903212014,-7)/
      DATA BERN  (  6) / -CONST(.528419013868749318484768220217955,-9)/
      DATA BERN  (  7) / +CONST(.133825365306846788328269809751291,-10)/
      DATA BERN  (  8) / -CONST(.338968029632258286683019539124944,-12)/
      DATA BERN  (  9) / +CONST(.858606205627784456413590545042562,-14)/
      DATA BERN  ( 10) / -CONST(.217486869855806187304151642386591,-15)/
      DATA BERN  ( 11) / +CONST(.550900282836022951520265260890225,-17)/
      DATA BERN  ( 12) / -CONST(.139544646858125233407076862640635,-18)/
      DATA BERN  ( 13) / +CONST(.353470703962946747169322997780379,-20)/
      DATA BERN  ( 14) / -CONST(.895351742703754685040261131811274,-22)/
      DATA BERN  ( 15) / +CONST(.226795245233768306031095073886816,-23)/
      DATA BERN  ( 16) / -CONST(.574472439520264523834847971943400,-24)/
      DATA BERN  ( 17) / +CONST(.145517247561486490186626486727132,-26)/
      DATA BERN  ( 18) / -CONST(.368599494066531017818178247990866,-28)/
      DATA BERN  ( 19) / +CONST(.933673425709504467203255515278562,-30)/
      DATA BERN  ( 20) / -CONST(.236502241570062993455963519636983,-31)/
      DATA PI / CONST(3.141592653589793238462643383279503,0)/
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DPOCH1
      IF (FIRST) THEN
         SQTBIG = CONST(1.0,0)/SQRT(CONST(24.0,0)*D1MACH(1))
         ALNEPS = LOG(D1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      IF (X.EQ.CONST(0.0,0)) DPOCH1 = DPSI(A)
      IF (X.EQ.CONST(0.0,0)) RETURN
C
      ABSX = ABS(X)
      ABSA = ABS(A)
      IF (ABSX.GT.CONST(0.1,0)*ABSA) GO TO 70
      IF (ABSX*LOG(MAX(ABSA,CONST(2.0,0))).GT.CONST(0.1,0)) GO TO 70
C
      BP = A
      IF (A.LT.(-CONST(0.5,0))) BP = CONST(1.0,0)- A - X
      INCR = 0
      IF (BP.LT.CONST(10.0,0)) INCR = CONST(11.0,0)- BP
      B = BP + INCR
C
      VAR = B + CONST(0.5,0)*(X-CONST(1.0,0))
      ALNVAR = LOG(VAR)
      Q = X*ALNVAR
C
      POLY1 = CONST(0.0,0)
      IF (VAR.GE.SQTBIG) GO TO 40
      VAR2 = (CONST(1.0,0)/VAR)**2
C
      RHO = CONST(0.5,0)*(X+CONST(1.0,0))
      GBERN(1) = CONST(1.0,0)
      GBERN(2) = -RHO/CONST(12.0,0)
      TERM = VAR2
      POLY1 = GBERN(2)*TERM
C
      NTERMS = -CONST(0.5,0)*ALNEPS/ALNVAR + CONST(1.0,0)
      IF (NTERMS .GT. 20) CALL XERMSG ('SLATEC', 'DPOCH1',
     +   'NTERMS IS TOO BIG, MAYBE D1MACH(3) IS BAD', 1, 2)
      IF (NTERMS.LT.2) GO TO 40
C
      DO 30 K=2,NTERMS
        GBK = CONST(0.0,0)
        DO 20 J=1,K
          NDX = K - J + 1
          GBK = GBK + BERN(NDX)*GBERN(J)
 20     CONTINUE
        GBERN(K+1) = -RHO*GBK/K
C
        TERM = TERM * (2*K-2-X)*(2*K-1-X)*VAR2
        POLY1 = POLY1 + GBERN(K+1)*TERM
 30   CONTINUE
C
 40   POLY1 = (X-CONST(1.0,0))*POLY1
      DPOCH1 = DEXPRL(Q)*(ALNVAR+Q*POLY1) + POLY1
C
      IF (INCR.EQ.0) GO TO 60
C
C WE HAVE DPOCH1(B,X), BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
C TO OBTAIN DPOCH1(BP,X).
C
      DO 50 II=1,INCR
        I = INCR - II
        BINV = CONST(1.0,0)/(BP+I)
        DPOCH1 = (DPOCH1 - BINV) / (CONST(1.0,0)+ X*BINV)
 50   CONTINUE
C
 60   IF (BP.EQ.A) RETURN
C
C WE HAVE DPOCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A REFLECTION
C FORMULA TO OBTAIN DPOCH1(A,X).
C
      SINPXX = SIN(PI*X)/X
      SINPX2 = SIN(CONST(0.5,0)*PI*X)
      TRIG = SINPXX*DCOT(PI*B) - CONST(2.0,0)*SINPX2*(SINPX2/X)
C
      DPOCH1 = TRIG + (CONST(1.0,0)+ X*TRIG)*DPOCH1
      RETURN
C
 70   DPOCH1 = (DPOCH(A,X) - CONST(1.0,0)) / X
      RETURN
C
      END
! 22Jun2000 fgtok -s r8_precision.sub "r8con.csh conversion"
! 22Jun2000 fgtok
! 22Jun2000 fgtok
