← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/dcd/grat1.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/dcd/grat1.f
 
1
      SUBROUTINE grat1(a,x,r,p,q,eps)
 
2
C     .. Scalar Arguments ..
 
3
      DOUBLE PRECISION a,eps,p,q,r,x
 
4
C     ..
 
5
C     .. Local Scalars ..
 
6
      DOUBLE PRECISION a2n,a2nm1,am0,an,an0,b2n,b2nm1,c,cma,g,h,j,l,sum,
 
7
     +                 t,tol,w,z
 
8
C     ..
 
9
C     .. External Functions ..
 
10
      DOUBLE PRECISION erf,erfc1,gam1,rexp
 
11
      EXTERNAL erf,erfc1,gam1,rexp
 
12
C     ..
 
13
C     .. Intrinsic Functions ..
 
14
      INTRINSIC abs,dlog,exp,sqrt
 
15
C     ..
 
16
C     .. Executable Statements ..
 
17
C-----------------------------------------------------------------------
 
18
C        EVALUATION OF THE INCOMPLETE GAMMA RATIO FUNCTIONS
 
19
C                      P(A,X) AND Q(A,X)
 
20
C
 
21
C     IT IS ASSUMED THAT A .LE. 1.  EPS IS THE TOLERANCE TO BE USED.
 
22
C     THE INPUT ARGUMENT R HAS THE VALUE E**(-X)*X**A/GAMMA(A).
 
23
C-----------------------------------------------------------------------
 
24
      IF (a*x.EQ.0.0D0) GO TO 120
 
25
      IF (a.EQ.0.5D0) GO TO 100
 
26
      IF (x.LT.1.1D0) GO TO 10
 
27
      GO TO 60
 
28
C
 
29
C             TAYLOR SERIES FOR P(A,X)/X**A
 
30
C
 
31
   10 an = 3.0D0
 
32
      c = x
 
33
      sum = x/ (a+3.0D0)
 
34
      tol = 0.1D0*eps/ (a+1.0D0)
 
35
   20 an = an + 1.0D0
 
36
      c = -c* (x/an)
 
37
      t = c/ (a+an)
 
38
      sum = sum + t
 
39
      IF (abs(t).GT.tol) GO TO 20
 
40
      j = a*x* ((sum/6.0D0-0.5D0/ (a+2.0D0))*x+1.0D0/ (a+1.0D0))
 
41
C
 
42
      z = a*dlog(x)
 
43
      h = gam1(a)
 
44
      g = 1.0D0 + h
 
45
      IF (x.LT.0.25D0) GO TO 30
 
46
      IF (a.LT.x/2.59D0) GO TO 50
 
47
      GO TO 40
 
48
 
 
49
   30 IF (z.GT.-.13394D0) GO TO 50
 
50
C
 
51
   40 w = exp(z)
 
52
      p = w*g* (0.5D0+ (0.5D0-j))
 
53
      q = 0.5D0 + (0.5D0-p)
 
54
      RETURN
 
55
C
 
56
   50 l = rexp(z)
 
57
      w = 0.5D0 + (0.5D0+l)
 
58
      q = (w*j-l)*g - h
 
59
      IF (q.LT.0.0D0) GO TO 90
 
60
      p = 0.5D0 + (0.5D0-q)
 
61
      RETURN
 
62
C
 
63
C              CONTINUED FRACTION EXPANSION
 
64
C
 
65
   60 a2nm1 = 1.0D0
 
66
      a2n = 1.0D0
 
67
      b2nm1 = x
 
68
      b2n = x + (1.0D0-a)
 
69
      c = 1.0D0
 
70
   70 a2nm1 = x*a2n + c*a2nm1
 
71
      b2nm1 = x*b2n + c*b2nm1
 
72
      am0 = a2nm1/b2nm1
 
73
      c = c + 1.0D0
 
74
      cma = c - a
 
75
      a2n = a2nm1 + cma*a2n
 
76
      b2n = b2nm1 + cma*b2n
 
77
      an0 = a2n/b2n
 
78
      IF (abs(an0-am0).GE.eps*an0) GO TO 70
 
79
      q = r*an0
 
80
      p = 0.5D0 + (0.5D0-q)
 
81
      RETURN
 
82
C
 
83
C                SPECIAL CASES
 
84
C
 
85
   80 p = 0.0D0
 
86
      q = 1.0D0
 
87
      RETURN
 
88
C
 
89
   90 p = 1.0D0
 
90
      q = 0.0D0
 
91
      RETURN
 
92
C
 
93
  100 IF (x.GE.0.25D0) GO TO 110
 
94
      p = erf(sqrt(x))
 
95
      q = 0.5D0 + (0.5D0-p)
 
96
      RETURN
 
97
 
 
98
  110 q = erfc1(0,sqrt(x))
 
99
      p = 0.5D0 + (0.5D0-q)
 
100
      RETURN
 
101
C
 
102
  120 IF (x.LE.a) GO TO 80
 
103
      GO TO 90
 
104
 
 
105
      END