OpenLibm/slatec/dpoch1.f

*DECK DPOCH1
      DOUBLE PRECISION 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      DOUBLE PRECISION (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
      DOUBLE PRECISION 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) / +.8333333333 3333333333 3333333333 333 D-1    /
      DATA BERN  (  2) / -.1388888888 8888888888 8888888888 888 D-2    /
      DATA BERN  (  3) / +.3306878306 8783068783 0687830687 830 D-4    /
      DATA BERN  (  4) / -.8267195767 1957671957 6719576719 576 D-6    /
      DATA BERN  (  5) / +.2087675698 7868098979 2100903212 014 D-7    /
      DATA BERN  (  6) / -.5284190138 6874931848 4768220217 955 D-9    /
      DATA BERN  (  7) / +.1338253653 0684678832 8269809751 291 D-10   /
      DATA BERN  (  8) / -.3389680296 3225828668 3019539124 944 D-12   /
      DATA BERN  (  9) / +.8586062056 2778445641 3590545042 562 D-14   /
      DATA BERN  ( 10) / -.2174868698 5580618730 4151642386 591 D-15   /
      DATA BERN  ( 11) / +.5509002828 3602295152 0265260890 225 D-17   /
      DATA BERN  ( 12) / -.1395446468 5812523340 7076862640 635 D-18   /
      DATA BERN  ( 13) / +.3534707039 6294674716 9322997780 379 D-20   /
      DATA BERN  ( 14) / -.8953517427 0375468504 0261131811 274 D-22   /
      DATA BERN  ( 15) / +.2267952452 3376830603 1095073886 816 D-23   /
      DATA BERN  ( 16) / -.5744724395 2026452383 4847971943 400 D-24   /
      DATA BERN  ( 17) / +.1455172475 6148649018 6626486727 132 D-26   /
      DATA BERN  ( 18) / -.3685994940 6653101781 8178247990 866 D-28   /
      DATA BERN  ( 19) / +.9336734257 0950446720 3255515278 562 D-30   /
      DATA BERN  ( 20) / -.2365022415 7006299345 5963519636 983 D-31   /
      DATA PI / 3.1415926535 8979323846 2643383279 503 D0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DPOCH1
      IF (FIRST) THEN
         SQTBIG = 1.0D0/SQRT(24.0D0*D1MACH(1))
         ALNEPS = LOG(D1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      IF (X.EQ.0.0D0) DPOCH1 = DPSI(A)
      IF (X.EQ.0.0D0) RETURN
C
      ABSX = ABS(X)
      ABSA = ABS(A)
      IF (ABSX.GT.0.1D0*ABSA) GO TO 70
      IF (ABSX*LOG(MAX(ABSA,2.0D0)).GT.0.1D0) GO TO 70
C
      BP = A
      IF (A.LT.(-0.5D0)) BP = 1.0D0 - A - X
      INCR = 0
      IF (BP.LT.10.0D0) INCR = 11.0D0 - BP
      B = BP + INCR
C
      VAR = B + 0.5D0*(X-1.0D0)
      ALNVAR = LOG(VAR)
      Q = X*ALNVAR
C
      POLY1 = 0.0D0
      IF (VAR.GE.SQTBIG) GO TO 40
      VAR2 = (1.0D0/VAR)**2
C
      RHO = 0.5D0*(X+1.0D0)
      GBERN(1) = 1.0D0
      GBERN(2) = -RHO/12.0D0
      TERM = VAR2
      POLY1 = GBERN(2)*TERM
C
      NTERMS = -0.5D0*ALNEPS/ALNVAR + 1.0D0
      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 = 0.0D0
        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-1.0D0)*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 = 1.0D0/(BP+I)
        DPOCH1 = (DPOCH1 - BINV) / (1.0D0 + 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(0.5D0*PI*X)
      TRIG = SINPXX*DCOT(PI*B) - 2.0D0*SINPX2*(SINPX2/X)
C
      DPOCH1 = TRIG + (1.0D0 + X*TRIG)*DPOCH1
      RETURN
C
 70   DPOCH1 = (DPOCH(A,X) - 1.0D0) / X
      RETURN
C
      END