1 /*
2 * Mathlib : A C Library of Special Functions
3 * Copyright (C) 1998 Ross Ihaka
4 * Copyright (C) 2000-8 The R Core Team
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, a copy is available at
18 * http://www.r-project.org/Licenses/
19 *
20 * DESCRIPTION
21 *
22 * The distribution function of the F distribution.
23 */
24
25 #include "nmath.h"
26 #include "dpq.h"
27
pf(double x,double df1,double df2,int lower_tail,int log_p)28 double pf(double x, double df1, double df2, int lower_tail, int log_p)
29 {
30 #ifdef IEEE_754
31 if (ISNAN(x) || ISNAN(df1) || ISNAN(df2))
32 return x + df2 + df1;
33 #endif
34 if (df1 <= 0. || df2 <= 0.) ML_ERR_return_NAN;
35
36 R_P_bounds_01(x, 0., ML_POSINF);
37
38 /* move to pchisq for very large values - was 'df1 > 4e5' in 2.0.x,
39 now only needed for df1 = Inf or df2 = Inf {since pbeta(0,*)=0} : */
40 if (df2 == ML_POSINF) {
41 if (df1 == ML_POSINF) {
42 if(x < 1.) return R_DT_0;
43 if(x == 1.) return (log_p ? -M_LN2 : 0.5);
44 if(x > 1.) return R_DT_1;
45 }
46
47 return pchisq(x * df1, df1, lower_tail, log_p);
48 }
49
50 if (df1 == ML_POSINF)/* was "fudge" 'df1 > 4e5' in 2.0.x */
51 return pchisq(df2 / x , df2, !lower_tail, log_p);
52
53 /* Avoid squeezing pbeta's first parameter against 1 : */
54 if (df1 * x > df2)
55 x = pbeta(df2 / (df2 + df1 * x), df2 / 2., df1 / 2.,
56 !lower_tail, log_p);
57 else
58 x = pbeta(df1 * x / (df2 + df1 * x), df1 / 2., df2 / 2.,
59 lower_tail, log_p);
60
61 return ML_VALID(x) ? x : ML_NAN;
62 }
63