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 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
|
DOUBLE PRECISION FUNCTION gam1(a)
C ------------------------------------------------------------------
C COMPUTATION OF 1/GAMMA(A+1) - 1 FOR -0.5 .LE. A .LE. 1.5
C ------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a
C ..
C .. Local Scalars ..
DOUBLE PRECISION bot,d,s1,s2,t,top,w
C ..
C .. Local Arrays ..
DOUBLE PRECISION p(7),q(5),r(9)
C ..
C .. Data statements ..
C -------------------
C -------------------
C -------------------
C -------------------
DATA p(1)/.577215664901533D+00/,p(2)/-.409078193005776D+00/,
+ p(3)/-.230975380857675D+00/,p(4)/.597275330452234D-01/,
+ p(5)/.766968181649490D-02/,p(6)/-.514889771323592D-02/,
+ p(7)/.589597428611429D-03/
DATA q(1)/.100000000000000D+01/,q(2)/.427569613095214D+00/,
+ q(3)/.158451672430138D+00/,q(4)/.261132021441447D-01/,
+ q(5)/.423244297896961D-02/
DATA r(1)/-.422784335098468D+00/,r(2)/-.771330383816272D+00/,
+ r(3)/-.244757765222226D+00/,r(4)/.118378989872749D+00/,
+ r(5)/.930357293360349D-03/,r(6)/-.118290993445146D-01/,
+ r(7)/.223047661158249D-02/,r(8)/.266505979058923D-03/,
+ r(9)/-.132674909766242D-03/
DATA s1/.273076135303957D+00/,s2/.559398236957378D-01/
C ..
C .. Executable Statements ..
C -------------------
t = a
d = a - 0.5D0
IF (d.GT.0.0D0) t = d - 0.5D0
IF (t.lt.0) GO TO 40
IF (t.eq.0) GO TO 10
GO TO 20
C
10 gam1 = 0.0D0
RETURN
C
20 top = (((((p(7)*t+p(6))*t+p(5))*t+p(4))*t+p(3))*t+p(2))*t + p(1)
bot = (((q(5)*t+q(4))*t+q(3))*t+q(2))*t + 1.0D0
w = top/bot
IF (d.GT.0.0D0) GO TO 30
gam1 = a*w
RETURN
30 gam1 = (t/a)* ((w-0.5D0)-0.5D0)
RETURN
C
40 top = (((((((r(9)*t+r(8))*t+r(7))*t+r(6))*t+r(5))*t+r(4))*t+r(3))*
+ t+r(2))*t + r(1)
bot = (s2*t+s1)*t + 1.0D0
w = top/bot
IF (d.GT.0.0D0) GO TO 50
gam1 = a* ((w+0.5D0)+0.5D0)
RETURN
50 gam1 = t*w/a
RETURN
END
|