This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / g / chisin.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:42  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 C     This corresponds to CHISIN,IF=DOUBLE and CHISIN64,IF=-DOUBLE
11       FUNCTION CHISIN(Q,N)
12
13 C     Computes the inverse of the
14 C     chi**2 distribution with n degrees of freedom.
15 C     Based on R.B. Goldstein,  ALGORITHM 451 Chi-Square Quantiles,
16 C     Collected Algorithms from CACM
17 C     Note the complementary definition of the integral!
18
19       CHARACTER NAME*(*)
20       CHARACTER*80 ERRTXT
21       PARAMETER (NAME = 'CHISIN')
22       DIMENSION C1(7),C2(7),C3(7),A(19)
23       PARAMETER (Z1 = 1, HF = Z1/2)
24
25       DATA C1
26
27      1/ 1.565326E-3, 1.060438E-3,-6.950356E-3,-1.323293E-2, 2.277679E-2,
28      2 -8.986007E-3,-1.513904E-2/
29       DATA C2
30      1/-1.450117E-3, 2.530010E-3, 5.169654E-3,-1.153761E-2, 1.128186E-2,
31      2  2.607083E-2,-2.237368E-1/
32       DATA C3
33      1/ 9.780499E-5,-8.426812E-4, 3.125580E-3,-8.553069E-3, 1.348028E-4,
34      2  4.713941E-1, 1.0000886/
35
36       DATA A
37      1/ 1.264616E-2,-1.425296E-2, 1.400483E-2,-5.886090E-3,-1.091214E-2,
38      2 -2.304527E-2, 3.135411E-3,-2.728484E-4,-9.699681E-3, 1.316872E-2,
39      3  2.618914E-2,-2.222222E-1, 5.406674E-5, 3.483789E-5,-7.274761E-4,
40      4  3.292181E-3,-8.729713E-3, 4.714045E-1, 1/
41
42       IF(Q .LT. 0 .OR. Q .GE. 1) THEN
43        H=0
44        WRITE(ERRTXT,101) Q
45        CALL MTLPRT(NAME,'G101.1',ERRTXT)
46       ELSEIF(N .LE. 0) THEN
47        H=0
48        WRITE(ERRTXT,102) N
49        CALL MTLPRT(NAME,'G101.2',ERRTXT)
50       ELSEIF(Q .EQ. 0) THEN
51        H=0
52       ELSEIF(N .EQ. 1) THEN
53        H=GAUSIN(HF*(1-Q))**2
54       ELSEIF(N .EQ. 2) THEN
55        H=-2*LOG(1-Q)
56       ELSE
57        F1=Z1/N
58        T=GAUSIN(Q)
59        F2=SQRT(F1)*T
60        IF(N .LT. 2+INT(4*ABS(T))) THEN
61         S1=C1(1)
62         S2=C2(1)
63         S3=C3(1)
64         DO 1 I = 2,7
65         S1=C1(I)+S1*F2
66         S2=C2(I)+S2*F2
67     1   S3=C3(I)+S3*F2
68         H=(S1*F1+S2)*F1+S3
69        ELSE
70         H=(((A(1)+A(2)*F2)*F1+(((A(3)+A(4)*F2)*F2+A(5))*F2+A(6)))*F1
71      1    +(((((A(7)+A(8)*F2)*F2+A(9))*F2+A(10))*F2+A(11))*F2+A(12)))*F1
72      2    +(((((A(13)*F2+A(14))*F2+A(15))*F2+A(16))*F2+A(17))*F2*F2
73      3    +A(18))*F2+A(19)
74        ENDIF
75        H=N*H**3
76       ENDIF
77       CHISIN=H
78       RETURN
79   101 FORMAT('ILLEGAL ARGUMENT Q = ',1P,E15.5)
80   102 FORMAT('ILLEGAL DEGREE OF FREEDOM N = ',I5)
81       END