*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