relibc_openlibm/slatec/poch1.f

*DECK POCH1
      FUNCTION POCH1 (A, X)
C***BEGIN PROLOGUE  POCH1
C***PURPOSE  Calculate a generalization of Pochhammer's symbol starting
C            from first order.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C1, C7A
C***TYPE      SINGLE PRECISION (POCH1-S, DPOCH1-D)
C***KEYWORDS  FIRST ORDER, FNLIB, POCHHAMMER, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate a generalization of Pochhammer's symbol for special
C situations that require 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  COT, EXPREL, POCH, PSI, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770801  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  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  POCH1
      DIMENSION BERN(9), GBERN(10)
      LOGICAL FIRST
      EXTERNAL COT
      SAVE BERN, PI, SQTBIG, ALNEPS, FIRST
      DATA BERN( 1) /   .8333333333 3333333E-01 /
      DATA BERN( 2) /  -.1388888888 8888889E-02 /
      DATA BERN( 3) /   .3306878306 8783069E-04 /
      DATA BERN( 4) /  -.8267195767 1957672E-06 /
      DATA BERN( 5) /   .2087675698 7868099E-07 /
      DATA BERN( 6) /  -.5284190138 6874932E-09 /
      DATA BERN( 7) /   .1338253653 0684679E-10 /
      DATA BERN( 8) /  -.3389680296 3225829E-12 /
      DATA BERN( 9) /   .8586062056 2778446E-14 /
      DATA PI / 3.1415926535 8979324 E0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  POCH1
      IF (FIRST) THEN
         SQTBIG = 1.0/SQRT(24.0*R1MACH(1))
         ALNEPS = LOG(R1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      IF (X.EQ.0.0) POCH1 = PSI(A)
      IF (X.EQ.0.0) RETURN
C
      ABSX = ABS(X)
      ABSA = ABS(A)
      IF (ABSX.GT.0.1*ABSA) GO TO 70
      IF (ABSX*LOG(MAX(ABSA,2.0)).GT.0.1) GO TO 70
C
      BP = A
      IF (A.LT.(-0.5)) BP = 1.0 - A - X
      INCR = 0
      IF (BP.LT.10.0) INCR = 11.0 - BP
      B = BP + INCR
C
      VAR = B + 0.5*(X-1.0)
      ALNVAR = LOG(VAR)
      Q = X*ALNVAR
C
      POLY1 = 0.0
      IF (VAR.GE.SQTBIG) GO TO 40
      VAR2 = (1.0/VAR)**2
C
      RHO = 0.5*(X+1.0)
      GBERN(1) = 1.0
      GBERN(2) = -RHO/12.0
      TERM = VAR2
      POLY1 = GBERN(2)*TERM
C
      NTERMS = -0.5*ALNEPS/ALNVAR + 1.0
      IF (NTERMS .GT. 9) CALL XERMSG ('SLATEC', 'POCH1',
     +   'NTERMS IS TOO BIG, MAYBE R1MACH(3) IS BAD', 1, 2)
      IF (NTERMS.LT.2) GO TO 40
C
      DO 30 K=2,NTERMS
        GBK = 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-1.0)*POLY1
      POCH1 = EXPREL(Q)*(ALNVAR + Q*POLY1) + POLY1
C
      IF (INCR.EQ.0) GO TO 60
C
C WE HAVE POCH1(B,X).  BUT BP IS SMALL, SO WE USE BACKWARDS RECURSION
C TO OBTAIN POCH1(BP,X).
C
      DO 50 II=1,INCR
        I = INCR - II
        BINV = 1.0/(BP+I)
        POCH1 = (POCH1-BINV)/(1.0+X*BINV)
 50   CONTINUE
C
 60   IF (BP.EQ.A) RETURN
C
C WE HAVE POCH1(BP,X), BUT A IS LT -0.5.  WE THEREFORE USE A REFLECTION
C FORMULA TO OBTAIN POCH1(A,X).
C
      SINPXX = SIN(PI*X)/X
      SINPX2 = SIN(0.5*PI*X)
      TRIG = SINPXX*COT(PI*B) - 2.0*SINPX2*(SINPX2/X)
C
      POCH1 = TRIG + (1.0 + X*TRIG) * POCH1
      RETURN
C
 70   POCH1 = (POCH(A,X) - 1.0) / X
      RETURN
C
      END