1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
|
DOUBLE PRECISION FUNCTION apser(a,b,x,eps)
C-----------------------------------------------------------------------
C APSER YIELDS THE INCOMPLETE BETA RATIO I(SUB(1-X))(B,A) FOR
C A .LE. MIN(EPS,EPS*B), B*X .LE. 1, AND X .LE. 0.5. USED WHEN
C A IS VERY SMALL. USE ONLY IF ABOVE INEQUALITIES ARE SATISFIED.
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,eps,x
C ..
C .. Local Scalars ..
DOUBLE PRECISION aj,bx,c,g,j,s,t,tol
C ..
C .. External Functions ..
DOUBLE PRECISION psi
EXTERNAL psi
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,dlog
C ..
C .. Data statements ..
C--------------------
DATA g/.577215664901533D0/
C ..
C .. Executable Statements ..
C--------------------
bx = b*x
t = x - bx
IF (b*eps.GT.2.D-2) GO TO 10
c = dlog(x) + psi(b) + g + t
GO TO 20
10 c = dlog(bx) + g + t
C
20 tol = 5.0D0*eps*abs(c)
j = 1.0D0
s = 0.0D0
30 j = j + 1.0D0
t = t* (x-bx/j)
aj = t/j
s = s + aj
IF (abs(aj).GT.tol) GO TO 30
C
apser = -a* (c+s)
RETURN
END
|