* Bug #7296 fixed - Enabling %nan, %inf and -%inf for the cdf* functions
[scilab.git] / scilab / modules / statistics / src / dcdflib / cdffnc.f
1       SUBROUTINE cdffnc(which,p,q,f,dfn,dfd,phonc,status,bound)
2 C**********************************************************************
3 C
4 C      SUBROUTINE CDFFNC( WHICH, P, Q, F, DFN, DFD, PNONC, STATUS, BOUND
5 C               Cumulative Distribution Function
6 C               Non-central F distribution
7 C
8 C
9 C                              Function
10 C
11 C
12 C     Calculates any one parameter of the Non-central F
13 C     distribution given values for the others.
14 C
15 C
16 C                              Arguments
17 C
18 C
19 C     WHICH --> Integer indicating which of the next five argument
20 C               values is to be calculated from the others.
21 C               Legal range: 1..5
22 C               iwhich = 1 : Calculate P and Q from F,DFN,DFD and PNONC
23 C               iwhich = 2 : Calculate F from P,Q,DFN,DFD and PNONC
24 C               iwhich = 3 : Calculate DFN from P,Q,F,DFD and PNONC
25 C               iwhich = 4 : Calculate DFD from P,Q,F,DFN and PNONC
26 C               iwhich = 5 : Calculate PNONC from P,Q,F,DFN and DFD
27 C                    INTEGER WHICH
28 C
29 C       P <--> The integral from 0 to F of the non-central f-density.
30 C              Input range: [0,1-1E-16).
31 C                    DOUBLE PRECISION P
32 C
33 C       Q <--> 1-P.
34 C            Q is not used by this subroutine and is only included
35 C            for similarity with other cdf* routines.
36 C                    DOUBLE PRECISION Q
37 C
38 C       F <--> Upper limit of integration of the non-central f-density.
39 C              Input range: [0, +infinity).
40 C              Search range: [0,1E300]
41 C                    DOUBLE PRECISION F
42 C
43 C     DFN < --> Degrees of freedom of the numerator sum of squares.
44 C               Input range: (0, +infinity).
45 C               Search range: [ 1E-300, 1E300]
46 C                    DOUBLE PRECISION DFN
47 C
48 C     DFD < --> Degrees of freedom of the denominator sum of squares.
49 C               Must be in range: (0, +infinity).
50 C               Input range: (0, +infinity).
51 C               Search range: [ 1E-300, 1E300]
52 C                    DOUBLE PRECISION DFD
53 C
54 C     PNONC <-> The non-centrality parameter
55 C               Input range: [0,infinity)
56 C               Search range: [0,1E4]
57 C                    DOUBLE PRECISION PHONC
58 C
59 C     STATUS <-- 0 if calculation completed correctly
60 C               -I if input parameter number I is out of range
61 C                1 if answer appears to be lower than lowest
62 C                  search bound
63 C                2 if answer appears to be higher than greatest
64 C                  search bound
65 C                3 if P + Q .ne. 1
66 C                    INTEGER STATUS
67 C
68 C     BOUND <-- Undefined if STATUS is 0
69 C
70 C               Bound exceeded by parameter number I if STATUS
71 C               is negative.
72 C
73 C               Lower search bound if STATUS is 1.
74 C
75 C               Upper search bound if STATUS is 2.
76 C
77 C
78 C                              Method
79 C
80 C
81 C     Formula  26.6.20   of   Abramowitz   and   Stegun,  Handbook  of
82 C     Mathematical  Functions (1966) is used to compute the cumulative
83 C     distribution function.
84 C
85 C     Computation of other parameters involve a seach for a value that
86 C     produces  the desired  value  of P.   The search relies  on  the
87 C     monotinicity of P with the other parameter.
88 C
89 C                            WARNING
90 C
91 C     The computation time  required for this  routine is proportional
92 C     to the noncentrality  parameter  (PNONC).  Very large  values of
93 C     this parameter can consume immense  computer resources.  This is
94 C     why the search range is bounded by 10,000.
95 C
96 C                              WARNING
97 C
98 C     The  value  of the  cumulative  noncentral F distribution is not
99 C     necessarily monotone in either degrees  of freedom.  There  thus
100 C     may be two values that provide a given  CDF value.  This routine
101 C     assumes monotonicity  and will find  an arbitrary one of the two
102 C     values.
103 C
104 C**********************************************************************
105 C     .. Parameters ..
106       DOUBLE PRECISION tent4
107       PARAMETER (tent4=1.0D4)
108       DOUBLE PRECISION tol
109       PARAMETER (tol=1.0D-13)
110       DOUBLE PRECISION atol
111       PARAMETER (atol=1.0D-50)
112       DOUBLE PRECISION zero,one,inf
113       PARAMETER (zero=1.0D-300,one=1.0D0-1.0D-16,inf=1.0D300)
114 C     ..
115 C     .. Scalar Arguments ..
116       DOUBLE PRECISION bound,dfd,dfn,f,p,q,phonc
117       INTEGER status,which
118 C     ..
119 C     .. Local Scalars ..
120       DOUBLE PRECISION fx,cum,ccum
121       LOGICAL qhi,qleft
122 C     ..
123 C     .. External Functions ..
124       INTEGER vfinite
125 C     ..
126 C     .. External Subroutines ..
127       EXTERNAL dinvr,dstinv,cumfnc
128 C     ..
129 C     .. Executable Statements ..
130 C
131 C     Check arguments
132 C
133       IF (.NOT. ((which.LT.1).OR. (which.GT.5))) GO TO 30
134       IF (.NOT. (which.LT.1)) GO TO 10
135       bound = 1.0D0
136       GO TO 20
137
138    10 bound = 5.0D0
139    20 status = -1
140       RETURN
141
142    30 IF (which.EQ.1) GO TO 70
143 C
144 C     P
145 C
146       IF (ISANAN(p).EQ.1) THEN
147          CALL RETURNANANFORTRAN(dfn)
148          CALL RETURNANANFORTRAN(dfd)
149          CALL RETURNANANFORTRAN(f)
150          CALL RETURNANANFORTRAN(phonc)
151          RETURN
152       ENDIF
153       IF (.NOT. ((p.LT.0.0D0).OR. (p.GT.one))) GO TO 60
154       IF (.NOT. (p.LT.0.0D0)) GO TO 40
155       bound = 0.0D0
156       GO TO 50
157
158    40 bound = one
159    50 status = -2
160       RETURN
161
162    60 CONTINUE
163    70 IF (which.EQ.2) GO TO 90
164 C
165 C     F
166 C
167       IF (ISANAN(f).EQ.1) THEN
168          CALL RETURNANANFORTRAN(p)
169          CALL RETURNANANFORTRAN(q)
170          CALL RETURNANANFORTRAN(dfn)
171          CALL RETURNANANFORTRAN(dfd)
172          CALL RETURNANANFORTRAN(phonc)
173          RETURN
174       ENDIF
175       IF (vfinite(1,f).EQ.0) then
176          IF (which.EQ.1) then
177             IF (f.GT.0) then
178                p = 1
179                q = 0
180                RETURN
181             ENDIF
182          ELSE
183             f = SIGN(1D300,f)
184          ENDIF
185       ENDIF
186       IF (.NOT. (f.LT.0.0D0)) GO TO 80
187       bound = 0.0D0
188       status = -4
189       RETURN
190
191    80 CONTINUE
192    90 IF (which.EQ.3) GO TO 110
193 C
194 C     DFN
195 C
196       IF (ISANAN(dfn).EQ.1) THEN
197          CALL RETURNANANFORTRAN(p)
198          CALL RETURNANANFORTRAN(q)
199          CALL RETURNANANFORTRAN(f)
200          CALL RETURNANANFORTRAN(dfd)
201          CALL RETURNANANFORTRAN(phonc)
202          RETURN
203       ENDIF
204       IF (vfinite(1,dfn).EQ.0) dfn = SIGN(inf,dfn)
205       IF (.NOT. (dfn.LE.0.0D0)) GO TO 100
206       bound = 0.0D0
207       status = -5
208       RETURN
209
210   100 CONTINUE
211   110 IF (which.EQ.4) GO TO 130
212 C
213 C     DFD
214 C
215       IF (ISANAN(dfd).EQ.1) THEN
216          CALL RETURNANANFORTRAN(p)
217          CALL RETURNANANFORTRAN(q)
218          CALL RETURNANANFORTRAN(dfn)
219          CALL RETURNANANFORTRAN(f)
220          CALL RETURNANANFORTRAN(phonc)
221          RETURN
222       ENDIF
223       IF (vfinite(1,dfd).EQ.0) dfd = SIGN(inf,dfd)
224       IF (.NOT. (dfd.LE.0.0D0)) GO TO 120
225       bound = 0.0D0
226       status = -6
227       RETURN
228
229   120 CONTINUE
230   130 IF (which.EQ.5) GO TO 150
231 C
232 C     PHONC
233 C
234       IF (ISANAN(phonc).EQ.1) THEN
235          CALL RETURNANANFORTRAN(p)
236          CALL RETURNANANFORTRAN(q)
237          CALL RETURNANANFORTRAN(dfn)
238          CALL RETURNANANFORTRAN(dfd)
239          CALL RETURNANANFORTRAN(f)
240          RETURN
241       ENDIF
242       IF (vfinite(1,phonc).EQ.0) phonc = SIGN(inf,phonc)
243       IF (.NOT. (phonc.LT.0.0D0)) GO TO 140
244       bound = 0.0D0
245       status = -7
246       RETURN
247
248   140 CONTINUE
249 C
250 C     Calculate ANSWERS
251 C
252   150 IF ((1).EQ. (which)) THEN
253 C
254 C     Calculating P
255 C
256           CALL cumfnc(f,dfn,dfd,phonc,p,q)
257           status = 0
258
259       ELSE IF ((2).EQ. (which)) THEN
260 C
261 C     Calculating F
262 C
263           f = 5.0D0
264           CALL dstinv(0.0D0,inf,0.5D0,0.5D0,5.0D0,atol,tol)
265           status = 0
266           CALL dinvr(status,f,fx,qleft,qhi)
267   160     IF (.NOT. (status.EQ.1)) GO TO 170
268           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
269           fx = cum - p
270           CALL dinvr(status,f,fx,qleft,qhi)
271           GO TO 160
272
273   170     IF (.NOT. (status.EQ.-1)) GO TO 200
274           IF (.NOT. (qleft)) GO TO 180
275           status = 1
276           bound = 0.0D0
277           GO TO 190
278
279   180     status = 2
280           bound = inf
281   190     CONTINUE
282   200     CONTINUE
283
284       ELSE IF ((3).EQ. (which)) THEN
285 C
286 C     Calculating DFN
287 C
288           dfn = 5.0D0
289           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
290           status = 0
291           CALL dinvr(status,dfn,fx,qleft,qhi)
292   210     IF (.NOT. (status.EQ.1)) GO TO 220
293           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
294           fx = cum - p
295           CALL dinvr(status,dfn,fx,qleft,qhi)
296           GO TO 210
297
298   220     IF (.NOT. (status.EQ.-1)) GO TO 250
299           IF (.NOT. (qleft)) GO TO 230
300           status = 1
301           bound = zero
302           GO TO 240
303
304   230     status = 2
305           bound = inf
306   240     CONTINUE
307   250     CONTINUE
308
309       ELSE IF ((4).EQ. (which)) THEN
310 C
311 C     Calculating DFD
312 C
313           dfd = 5.0D0
314           CALL dstinv(zero,inf,0.5D0,0.5D0,5.0D0,atol,tol)
315           status = 0
316           CALL dinvr(status,dfd,fx,qleft,qhi)
317   260     IF (.NOT. (status.EQ.1)) GO TO 270
318           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
319           fx = cum - p
320           CALL dinvr(status,dfd,fx,qleft,qhi)
321           GO TO 260
322
323   270     IF (.NOT. (status.EQ.-1)) GO TO 300
324           IF (.NOT. (qleft)) GO TO 280
325           status = 1
326           bound = zero
327           GO TO 290
328
329   280     status = 2
330           bound = inf
331   290     CONTINUE
332   300     CONTINUE
333
334       ELSE IF ((5).EQ. (which)) THEN
335 C
336 C     Calculating PHONC
337 C
338           phonc = 5.0D0
339           CALL dstinv(0.0D0,tent4,0.5D0,0.5D0,5.0D0,atol,tol)
340           status = 0
341           CALL dinvr(status,phonc,fx,qleft,qhi)
342   310     IF (.NOT. (status.EQ.1)) GO TO 320
343           CALL cumfnc(f,dfn,dfd,phonc,cum,ccum)
344           fx = cum - p
345           CALL dinvr(status,phonc,fx,qleft,qhi)
346           GO TO 310
347
348   320     IF (.NOT. (status.EQ.-1)) GO TO 350
349           IF (.NOT. (qleft)) GO TO 330
350           status = 1
351           bound = 0.0D0
352           GO TO 340
353
354   330     status = 2
355           bound = tent4
356   340     CONTINUE
357   350 END IF
358
359       RETURN
360
361       END