This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / g / prob.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:41  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       FUNCTION PROB(X,N)
11  
12 #include "gen/imp64.inc"
13       REAL PROB,X
14       CHARACTER NAME*(*)
15       CHARACTER*80 ERRTXT
16       PARAMETER (NAME = 'PROB')
17       PARAMETER (R1 = 1, HF = R1/2, TH = R1/3, F1 = 2*R1/9)
18       PARAMETER (C1 = 1.12837 91670 95513D0)
19       PARAMETER (NMAX = 300)
20 *                maximum chi2 per df for df >= 2., if chi2/df > chipdf prob=0.
21       PARAMETER (CHIPDF = 100.)
22       PARAMETER (XMAX = 174.673, XMAX2 = 2*XMAX)
23 #if defined(CERNLIB_IBM)
24 *
25 *     13.3 is limit of DERFC intrinsic (Wojciech Wojcik/IN2P3)
26 *
27       PARAMETER (XLIM = 13.3)
28 #endif
29 #if !defined(CERNLIB_IBM)
30       PARAMETER (XLIM = 24.)
31 #endif
32       PARAMETER (EPS = 1D-30)
33 #if defined(CERNLIB_DOUBLE)
34       GERFC(V)=DERFC(V)
35 #endif
36 #if !defined(CERNLIB_DOUBLE)
37       GERFC(V)= ERFC(V)
38 #endif
39  
40       Y=X
41       U=HF*Y
42       IF(N .LE. 0) THEN
43        H=0
44        WRITE(ERRTXT,101) N
45        CALL MTLPRT(NAME,'G100.1',ERRTXT)
46       ELSEIF(Y .LT. 0) THEN
47        H=0
48         WRITE(ERRTXT,102) X
49        CALL MTLPRT(NAME,'G100.2',ERRTXT)
50       ELSEIF(Y .EQ. 0 .OR. N/20 .GT. Y) THEN
51        H=1
52       ELSEIF(N .EQ. 1) THEN
53        W=SQRT(U)
54        IF(W .LT. XLIM) THEN
55         H=GERFC(W)
56        ELSE
57         H=0
58        ENDIF
59       ELSEIF(N .GT. NMAX) THEN
60        S=R1/N
61        T=F1*S
62        W=((Y*S)**TH-(1-T))/SQRT(2*T)
63        IF(W .LT. -XLIM) THEN
64         H=1
65        ELSEIF(W .LT. XLIM) THEN
66         H=HF*GERFC(W)
67        ELSE
68         H=0
69        ENDIF
70       ELSE
71        M=N/2
72        IF(U .LT. XMAX2 .AND. (Y/N).LE.CHIPDF ) THEN
73         S=EXP(-HF*U)
74         T=S
75         E=S
76         IF(2*M .EQ. N) THEN
77          FI=0
78          DO 1 I = 1,M-1
79          FI=FI+1
80          T=U*T/FI
81     1    S=S+T
82          H=S*E
83         ELSE
84          FI=1
85          DO 2 I=1,M-1
86          FI=FI+2
87          T=T*Y/FI
88     2    S=S+T
89          W=SQRT(U)
90          IF(W.LT.XLIM) THEN
91           H=C1*W*S*E+GERFC(W)
92          ELSE
93           H=0.
94          ENDIF
95         ENDIF
96        ELSE
97         H=0
98        ENDIF
99       ENDIF
100       IF ( H.GT. EPS ) THEN
101          PROB=H
102       ELSE
103          PROB=0.
104       ENDIF
105       RETURN
106   101 FORMAT('N = ',I6,' < 1')
107   102 FORMAT('X = ',1P,E20.10,' < 0')
108       END