1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 2006-8 The R Core Team
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * http://www.r-project.org/Licenses/
18 */
19
20 #include "nmath.h"
21 #include "dpq.h"
22
qnt(double p,double df,double ncp,int lower_tail,int log_p)23 double qnt(double p, double df, double ncp, int lower_tail, int log_p)
24 {
25 const static double accu = 1e-13;
26 const static double Eps = 1e-11; /* must be > accu */
27
28 double ux, lx, nx, pp;
29
30 #ifdef IEEE_754
31 if (ISNAN(p) || ISNAN(df) || ISNAN(ncp))
32 return p + df + ncp;
33 #endif
34 if (!R_FINITE(df)) ML_ERR_return_NAN;
35
36 /* Was
37 * df = floor(df + 0.5);
38 * if (df < 1 || ncp < 0) ML_ERR_return_NAN;
39 */
40 if (df <= 0.0) ML_ERR_return_NAN;
41
42 if(ncp == 0.0 && df >= 1.0) return qt(p, df, lower_tail, log_p);
43
44 R_Q_P01_boundaries(p, ML_NEGINF, ML_POSINF);
45
46 p = R_DT_qIv(p);
47
48 /* Invert pnt(.) :
49 * 1. finding an upper and lower bound */
50 if(p > 1 - DBL_EPSILON) return ML_POSINF;
51 pp = fmin2(1 - DBL_EPSILON, p * (1 + Eps));
52 for(ux = fmax2(1., ncp);
53 ux < DBL_MAX && pnt(ux, df, ncp, TRUE, FALSE) < pp;
54 ux *= 2);
55 pp = p * (1 - Eps);
56 for(lx = fmin2(-1., -ncp);
57 lx > -DBL_MAX && pnt(lx, df, ncp, TRUE, FALSE) > pp;
58 lx *= 2);
59
60 /* 2. interval (lx,ux) halving : */
61 do {
62 nx = 0.5 * (lx + ux);
63 if (pnt(nx, df, ncp, TRUE, FALSE) > p) ux = nx; else lx = nx;
64 }
65 while ((ux - lx) / fabs(nx) > accu);
66
67 return 0.5 * (lx + ux);
68 }
69