1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 1998--2017 The R Core Team
4 * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
5 * based on code (C) 1979 and later Royal Statistical Society
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, a copy is available at
19 * https://www.R-project.org/Licenses/
20 *
21
22 * Reference:
23 * Cran, G. W., K. J. Martin and G. E. Thomas (1977).
24 * Remark AS R19 and Algorithm AS 109,
25 * Applied Statistics, 26(1), 111-114.
26 * Remark AS R83 (v.39, 309-310) and the correction (v.40(1) p.236)
27 * have been incorporated in this version.
28 */
29
30
31 #include <R_ext/Arith.h>
32 #include "nmath.h"
33 #include "dpq.h"
34
35 #ifdef DEBUG_qbeta
36 # define R_ifDEBUG_printf(...) REprintf(__VA_ARGS__)
37 #else
38 # define R_ifDEBUG_printf(...)
39 #endif
40
41 #define USE_LOG_X_CUTOFF -5.
42 // --- based on some testing; had = -10
43
44 #define n_NEWTON_FREE 4
45 // --- based on some testing; had = 10
46
47 #define MLOGICAL_NA -1
48 // an "NA_LOGICAL" substitute for Mathlib {only used here, for now}
49
50 //attribute_hidden
51 static void
52 qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p,
53 int swap_01, double log_q_cut, int n_N, double* qb);
54
qbeta(double alpha,double p,double q,int lower_tail,int log_p)55 double qbeta(double alpha, double p, double q, int lower_tail, int log_p)
56 {
57
58 /* test for admissibility of parameters */
59 #ifdef IEEE_754
60 if (ISNAN(p) || ISNAN(q) || ISNAN(alpha))
61 return p + q + alpha;
62 #endif
63 if(p < 0. || q < 0.) ML_ERR_return_NAN;
64 // allowing p==0 and q==0 <==> treat as one- or two-point mass
65
66 double qbet[2];// = { qbeta(), 1 - qbeta() }
67 qbeta_raw(alpha, p, q, lower_tail, log_p,
68 MLOGICAL_NA, USE_LOG_X_CUTOFF, n_NEWTON_FREE, qbet);
69 return qbet[0];
70 }
71
72 static const double
73 #ifdef IEEE_754
74 // CARE: assumes subnormal numbers, i.e., no underflow at DBL_MIN:
75 DBL_very_MIN = DBL_MIN / 4.,
76 DBL_log_v_MIN = M_LN2*(DBL_MIN_EXP - 2),
77 // Too extreme: inaccuracy in pbeta(); e.g for qbeta(0.95, 1e-9, 20):
78 // -> in pbeta() --> bgrat(..... b*z == 0 underflow, hence inaccurate pbeta()
79 /* DBL_very_MIN = 0x0.0000001p-1022, // = 2^-1050 = 2^(-1022 - 28) */
80 /* DBL_log_v_MIN = -1050. * M_LN2, // = log(DBL_very_MIN) */
81 // the most extreme -- not ok, as pbeta() then behaves strangely,
82 // e.g., for qbeta(0.95, 1e-8, 20):
83 /* DBL_very_MIN = 0x0.0000000000001p-1022, // = 2^-1074 = 2^(-1022 -52) */
84 /* DBL_log_v_MIN = -1074. * M_LN2, // = log(DBL_very_MIN) */
85
86 DBL_1__eps = 0x1.fffffffffffffp-1; // = 1 - 2^-53
87 #else // untested :
88 DBL_1__eps = 1 - DBL_EPSILON; // or rather (1 - DBL_EPSILON/2) (??)
89 #endif
90
91 /* set the exponent of acu to -2r-2 for r digits of accuracy */
92 /*---- NEW ---- -- still fails for p = 1e11, q=.5*/
93
94 #define fpu 3e-308
95 /* acu_min: Minimal value for accuracy 'acu' which will depend on (a,p);
96 acu_min >= fpu ! */
97 #define acu_min 1e-300
98 #define p_lo fpu
99 #define p_hi 1-2.22e-16
100
101 #define const1 2.30753
102 #define const2 0.27061
103 #define const3 0.99229
104 #define const4 0.04481
105
106 // Returns both qbeta() and its "mirror" 1-qbeta(). Useful notably when qbeta() ~= 1
107 attribute_hidden void
qbeta_raw(double alpha,double p,double q,int lower_tail,int log_p,int swap_01,double log_q_cut,int n_N,double * qb)108 qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p,
109 int swap_01, // {TRUE, NA, FALSE}: if NA, algorithm decides swap_tail
110 double log_q_cut, /* if == Inf: return log(qbeta(..));
111 otherwise, if finite: the bound for
112 switching to log(x)-scale; see use_log_x */
113 int n_N, // number of "unconstrained" Newton steps before switching to constrained
114 double *qb) // = qb[0:1] = { qbeta(), 1 - qbeta() }
115 {
116 Rboolean
117 swap_choose = (swap_01 == MLOGICAL_NA),
118 swap_tail,
119 log_, give_log_q = (log_q_cut == ML_POSINF),
120 use_log_x = give_log_q, // or u < log_q_cut below
121 warned = FALSE, add_N_step = TRUE;
122 int i_pb, i_inn;
123 double a, la, logbeta, g, h, pp, p_, qq, r, s, t, w, y = -1.;
124 volatile double u, xinbta;
125
126 // Assuming p >= 0, q >= 0 here ...
127
128 // Deal with boundary cases here:
129 if(alpha == R_DT_0) {
130 #define return_q_0 \
131 if(give_log_q) { qb[0] = ML_NEGINF; qb[1] = 0; } \
132 else { qb[0] = 0; qb[1] = 1; } \
133 return
134
135 return_q_0;
136 }
137 if(alpha == R_DT_1) {
138 #define return_q_1 \
139 if(give_log_q) { qb[0] = 0; qb[1] = ML_NEGINF; } \
140 else { qb[0] = 1; qb[1] = 0; } \
141 return
142
143 return_q_1;
144 }
145
146 // check alpha {*before* transformation which may all accuracy}:
147 if((log_p && alpha > 0) ||
148 (!log_p && (alpha < 0 || alpha > 1))) { // alpha is outside
149 R_ifDEBUG_printf("qbeta(alpha=%g, %g, %g, .., log_p=%d): %s%s\n",
150 alpha, p,q, log_p, "alpha not in ",
151 log_p ? "[-Inf, 0]" : "[0,1]");
152 // ML_ERR_return_NAN :
153 ML_ERROR(ME_DOMAIN, "");
154 qb[0] = qb[1] = ML_NAN; return;
155 }
156
157 // p==0, q==0, p = Inf, q = Inf <==> treat as one- or two-point mass
158 if(p == 0 || q == 0 || !R_FINITE(p) || !R_FINITE(q)) {
159 // We know 0 < T(alpha) < 1 : pbeta() is constant and trivial in {0, 1/2, 1}
160 R_ifDEBUG_printf(
161 "qbeta(%g, %g, %g, lower_t=%d, log_p=%d): (p,q)-boundary: trivial\n",
162 alpha, p,q, lower_tail, log_p);
163 if(p == 0 && q == 0) { // point mass 1/2 at each of {0,1} :
164 if(alpha < R_D_half) { return_q_0; }
165 if(alpha > R_D_half) { return_q_1; }
166 // else: alpha == "1/2"
167 #define return_q_half \
168 if(give_log_q) qb[0] = qb[1] = -M_LN2; \
169 else qb[0] = qb[1] = 0.5; \
170 return
171
172 return_q_half;
173 } else if (p == 0 || p/q == 0) { // point mass 1 at 0 - "flipped around"
174 return_q_0;
175 } else if (q == 0 || q/p == 0) { // point mass 1 at 0 - "flipped around"
176 return_q_1;
177 }
178 // else: p = q = Inf : point mass 1 at 1/2
179 return_q_half;
180 }
181
182 /* initialize */
183 p_ = R_DT_qIv(alpha);/* lower_tail prob (in any case) */
184 // Conceptually, 0 < p_ < 1 (but can be 0 or 1 because of cancellation!)
185 logbeta = lbeta(p, q);
186
187 swap_tail = (swap_choose) ? (p_ > 0.5) : swap_01;
188 // change tail; default (swap_01 = NA): afterwards 0 < a <= 1/2
189 if(swap_tail) { /* change tail, swap p <-> q :*/
190 a = R_DT_CIv(alpha); // = 1 - p_ < 1/2
191 /* la := log(a), but without numerical cancellation: */
192 la = R_DT_Clog(alpha);
193 pp = q; qq = p;
194 }
195 else {
196 a = p_;
197 la = R_DT_log(alpha);
198 pp = p; qq = q;
199 }
200
201 /* calculate the initial approximation */
202
203 /* Desired accuracy for Newton iterations (below) should depend on (a,p)
204 * This is from Remark .. on AS 109, adapted.
205 * However, it's not clear if this is "optimal" for IEEE double prec.
206
207 * acu = fmax2(acu_min, pow(10., -25. - 5./(pp * pp) - 1./(a * a)));
208
209 * NEW: 'acu' accuracy NOT for squared adjustment, but simple;
210 * ---- i.e., "new acu" = sqrt(old acu)
211 */
212 double acu = fmax2(acu_min, pow(10., -13. - 2.5/(pp * pp) - 0.5/(a * a)));
213 // try to catch "extreme left tail" early
214 double tx, u0 = (la + log(pp) + logbeta) / pp; // = log(x_0)
215 static const double
216 log_eps_c = M_LN2 * (1. - DBL_MANT_DIG);// = log(DBL_EPSILON) = -36.04..
217 r = pp*(1.-qq)/(pp+1.);
218
219 t = 0.2;
220 // FIXME: Factor 0.2 is a bit arbitrary; '1' is clearly much too much.
221
222 R_ifDEBUG_printf(
223 "qbeta(%g, %g, %g, lower_t=%d, log_p=%d):%s\n"
224 " swap_tail=%d, la=%#8g, u0=%#8g (bnd: %g (%g)) ",
225 alpha, p,q, lower_tail, log_p,
226 (log_p && (p_ == 0. || p_ == 1.)) ? (p_==0.?" p_=0":" p_=1") : "",
227 swap_tail, la, u0,
228 (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2.,
229 t*log_eps_c - log(fabs(r))
230 );
231
232 if(M_LN2 * DBL_MIN_EXP < u0 && // cannot allow exp(u0) = 0 ==> exp(u1) = exp(u0) = 0
233 u0 < -0.01 && // (must: u0 < 0, but too close to 0 <==> x = exp(u0) = 0.99..)
234 // qq <= 2 && // <--- "arbitrary"
235 // u0 < t*log_eps_c - log(fabs(r)) &&
236 u0 < (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2.)
237 {
238 // TODO: maybe jump here from below, when initial u "fails" ?
239 // L_tail_u:
240 // MM's one-step correction (cheaper than 1 Newton!)
241 r = r*exp(u0);// = r*x0
242 if(r > -1.) {
243 u = u0 - log1p(r)/pp;
244 R_ifDEBUG_printf("u1-u0=%9.3g --> choosing u = u1\n", u-u0);
245 } else {
246 u = u0;
247 R_ifDEBUG_printf("cannot cheaply improve u0\n");
248 }
249 tx = xinbta = exp(u);
250 use_log_x = TRUE; // or (u < log_q_cut) ??
251 goto L_Newton;
252 }
253
254 // y := y_\alpha in AS 64 := Hastings(1955) approximation of qnorm(1 - a) :
255 r = sqrt(-2 * la);
256 y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r);
257
258 if (pp > 1 && qq > 1) { // use Carter(1947), see AS 109, remark '5.'
259 r = (y * y - 3.) / 6.;
260 s = 1. / (pp + pp - 1.);
261 t = 1. / (qq + qq - 1.);
262 h = 2. / (s + t);
263 w = y * sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h));
264 R_ifDEBUG_printf("p,q > 1 => w=%g", w);
265 if(w > 300) { // exp(w+w) is huge or overflows
266 t = w+w + log(qq) - log(pp); // = argument of log1pexp(.)
267 u = // log(xinbta) = - log1p(qq/pp * exp(w+w)) = -log(1 + exp(t))
268 (t <= 18) ? -log1p(exp(t)) : -t - exp(-t);
269 xinbta = exp(u);
270 } else {
271 xinbta = pp / (pp + qq * exp(w + w));
272 u = // log(xinbta)
273 - log1p(qq/pp * exp(w+w));
274 }
275 } else { // use the original AS 64 proposal, Scheffé-Tukey (1944) and Wilson-Hilferty
276 r = qq + qq;
277 /* A slightly more stable version of t := \chi^2_{alpha} of AS 64
278 * t = 1. / (9. * qq); t = r * R_pow_di(1. - t + y * sqrt(t), 3); */
279 t = 1. / (3. * sqrt(qq));
280 t = r * R_pow_di(1. + t*(-t + y), 3);// = \chi^2_{alpha} of AS 64
281 s = 4. * pp + r - 2.;// 4p + 2q - 2 = numerator of new t = (...) / chi^2
282 R_ifDEBUG_printf("min(p,q) <= 1: t=%g", t);
283 if (t == 0 || (t < 0. && s >= t)) { // cannot use chisq approx
284 // x0 = 1 - { (1-a)*q*B(p,q) } ^{1/q} {AS 65}
285 // xinbta = 1. - exp((log(1-a)+ log(qq) + logbeta) / qq);
286 double l1ma;/* := log(1-a), directly from alpha (as 'la' above):
287 * FIXME: not worth it? log1p(-a) always the same ?? */
288 if(swap_tail)
289 l1ma = R_DT_log(alpha);
290 else
291 l1ma = R_DT_Clog(alpha);
292 R_ifDEBUG_printf(" t <= 0 : log1p(-a)=%.15g, better l1ma=%.15g\n", log1p(-a), l1ma);
293 double xx = (l1ma + log(qq) + logbeta) / qq;
294 if(xx <= 0.) {
295 xinbta = -expm1(xx);
296 u = R_Log1_Exp (xx);// = log(xinbta) = log(1 - exp(...A...))
297 } else { // xx > 0 ==> 1 - e^xx < 0 .. is nonsense
298 R_ifDEBUG_printf(" xx=%g > 0: xinbta:= 1-e^xx < 0\n", xx);
299 xinbta = 0; u = ML_NEGINF; /// FIXME can do better?
300 }
301 } else {
302 t = s / t;
303 R_ifDEBUG_printf(" t > 0 or s < t < 0: new t = %g ( > 1 ?)\n", t);
304 if (t <= 1.) { // cannot use chisq, either
305 u = (la + log(pp) + logbeta) / pp;
306 xinbta = exp(u);
307 } else { // (1+x0)/(1-x0) = t, solved for x0 :
308 xinbta = 1. - 2. / (t + 1.);
309 u = log1p(-2. / (t + 1.));
310 }
311 }
312 }
313
314 // Problem: If initial u is completely wrong, we make a wrong decision here
315 if(swap_choose &&
316 (( swap_tail && u >= -exp( log_q_cut)) || // ==> "swap back"
317 (!swap_tail && u >= -exp(4*log_q_cut) && pp / qq < 1000.) // ==> "swap now"
318 )) {
319 // "revert swap" -- and use_log_x
320 swap_tail = !swap_tail;
321 R_ifDEBUG_printf(" u = %g (e^u = xinbta = %.16g) ==> ", u, xinbta);
322 if(swap_tail) { // "swap now" (much less easily)
323 a = R_DT_CIv(alpha); // needed ?
324 la = R_DT_Clog(alpha);
325 pp = q; qq = p;
326 }
327 else { // swap back :
328 a = p_;
329 la = R_DT_log(alpha);
330 pp = p; qq = q;
331 }
332 R_ifDEBUG_printf("\"%s\"; la = %g\n",
333 (swap_tail ? "swap now" : "swap back"), la);
334 // we could redo computations above, but this should be stable
335 u = R_Log1_Exp(u);
336 xinbta = exp(u);
337
338 /* Careful: "swap now" should not fail if
339 1) the above initial xinbta is "completely wrong"
340 2) The correction step can go outside (u_n > 0 ==> e^u > 1 is illegal)
341 e.g., for qbeta(0.2066, 0.143891, 0.05)
342 */
343 } else R_ifDEBUG_printf("\n");
344
345 if(!use_log_x)
346 use_log_x = (u < log_q_cut);// <==> xinbta = e^u < exp(log_q_cut)
347 Rboolean
348 bad_u = !R_FINITE(u),
349 bad_init = bad_u || xinbta > p_hi;
350
351 R_ifDEBUG_printf(" -> u = %g, e^u = xinbta = %.16g, (Newton acu=%g%s%s%s)\n",
352 u, xinbta, acu, (bad_u ? ", ** bad u **" : ""),
353 ((bad_init && !bad_u) ? ", ** bad_init **" : ""),
354 (use_log_x ? ", on u = LOG(x) SCALE" : ""));
355
356 double u_n = 1.; // -Wall
357 tx = xinbta; // keeping "original initial x" (for now)
358
359 if(bad_u || u < log_q_cut) {
360 /* e.g.
361 qbeta(0.21, .001, 0.05)
362 try "left border" quickly, i.e.,
363 try at smallest positive number: */
364 w = pbeta_raw(DBL_very_MIN, pp, qq, TRUE, log_p);
365 if(w > (log_p ? la : a)) {
366 R_ifDEBUG_printf(
367 " quantile is left of %g; \"convergence\"\n", DBL_very_MIN);
368 if(log_p || fabs(w - a) < fabs(0 - a)) { // DBL_very_MIN is better than 0
369 tx = DBL_very_MIN;
370 u_n = DBL_log_v_MIN;// = log(DBL_very_MIN)
371 } else {
372 tx = 0.;
373 u_n = ML_NEGINF;
374 }
375 use_log_x = log_p; add_N_step = FALSE; goto L_return;
376 }
377 else {
378 R_ifDEBUG_printf(" pbeta(%g, *) = %g <= %g (= %s) --> continuing\n",
379 DBL_log_v_MIN, w, (log_p ? la : a), (log_p ? "la" : "a"));
380 if(u < DBL_log_v_MIN) {
381 u = DBL_log_v_MIN;// = log(DBL_very_MIN)
382 xinbta = DBL_very_MIN;
383 }
384 }
385 }
386
387
388 /* Sometimes the approximation is negative (and == 0 is also not "ok") */
389 if (bad_init && !(use_log_x && tx > 0)) {
390 if(u == ML_NEGINF) {
391 R_ifDEBUG_printf(" u = -Inf;");
392 u = M_LN2 * DBL_MIN_EXP;
393 xinbta = DBL_MIN;
394 } else {
395 R_ifDEBUG_printf(" bad_init: u=%g, xinbta=%g;", u,xinbta);
396 xinbta = (xinbta > 1.1) // i.e. "way off"
397 ? 0.5 // otherwise, keep the respective boundary:
398 : ((xinbta < p_lo) ? exp(u) : p_hi);
399 if(bad_u)
400 u = log(xinbta);
401 // otherwise: not changing "potentially better" u than the above
402 }
403 R_ifDEBUG_printf(" -> (partly)new u=%g, xinbta=%g\n", u,xinbta);
404 }
405
406 L_Newton:
407 /* --------------------------------------------------------------------
408
409 * Solve for x by a modified Newton-Raphson method, using pbeta_raw()
410 */
411 r = 1 - pp;
412 t = 1 - qq;
413 double wprev = 0., prev = 1., adj = 1.; // -Wall
414
415 if(use_log_x) { // find log(xinbta) -- work in u := log(x) scale
416 // if(bad_init && tx > 0) xinbta = tx;// may have been better
417
418 for (i_pb=0; i_pb < 1000; i_pb++) {
419 // using log_p == TRUE unconditionally here
420 /* FIXME: if exp(u) = xinbta underflows to 0,
421 * want different formula pbeta_log(u, ..) */
422 y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, TRUE);
423
424 /* w := Newton step size for L(u) = log F(e^u) =!= 0; u := log(x)
425 * = (L(.) - la) / L'(.); L'(u)= (F'(e^u) * e^u ) / F(e^u)
426 * = (L(.) - la)*F(.) / {F'(e^u) * e^u } =
427 * = (L(.) - la) * e^L(.) * e^{-log F'(e^u) - u}
428 * = ( y - la) * e^{ y - u -log F'(e^u)}
429 and -log F'(x)= -log f(x) = - -logbeta + (1-p) log(x) + (1-q) log(1-x)
430 = logbeta + (1-p) u + (1-q) log(1-e^u)
431 */
432 w = (y == ML_NEGINF) // y = -Inf well possible: we are on log scale!
433 ? 0. : (y - la) * exp(y - u + logbeta + r * u + t * R_Log1_Exp(u));
434 if(!R_FINITE(w))
435 break;
436 if (i_pb >= n_N && w * wprev <= 0.)
437 prev = fmax2(fabs(adj),fpu);
438 R_ifDEBUG_printf(
439 "N(i=%2d): u=%#20.16g, pb(e^u)=%#15.9g, w=%#15.9g, %s prev=%g,",
440 i_pb, u, y, w,
441 (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev);
442 g = 1;
443 for (i_inn=0; i_inn < 1000; i_inn++) {
444 adj = g * w;
445 // safe guard (here, from the very beginning)
446 if (fabs(adj) < prev) {
447 u_n = u - adj; // u_{n+1} = u_n - g*w
448 if (u_n <= 0.) { // <==> 0 < xinbta := e^u <= 1
449 if (prev <= acu || fabs(w) <= acu) {
450 R_ifDEBUG_printf(
451 " it{in}=%d, -adj=%g, %s <= acu ==> convergence\n",
452 i_inn, -adj, (prev <= acu) ? "prev" : "|w|");
453 goto L_converged;
454 }
455 // if (u_n != ML_NEGINF && u_n != 1)
456 break;
457 }
458 }
459 g /= 3;
460 }
461 // (cancellation in (u_n -u) => may differ from adj:
462 double D = fmin2(fabs(adj), fabs(u_n - u));
463 /* R_ifDEBUG_printf(" delta(u)=%g\n", u_n - u); */
464 R_ifDEBUG_printf(" it{in}=%d, delta(u)=%9.3g, D/|.|=%.3g\n",
465 i_inn, u_n - u, D/fabs(u_n + u));
466 if (D <= 4e-16 * fabs(u_n + u))
467 goto L_converged;
468 u = u_n;
469 xinbta = exp(u);
470 wprev = w;
471 } // for(i )
472
473 } else { // "normal scale" Newton
474
475 for (i_pb=0; i_pb < 1000; i_pb++) {
476 y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p);
477 // delta{y} : d_y = y - (log_p ? la : a);
478 #ifdef IEEE_754
479 if(!R_FINITE(y) && !(log_p && y == ML_NEGINF))// y = -Inf is ok if(log_p)
480 #else
481 if (errno)
482 #endif
483 { // ML_ERR_return_NAN :
484 ML_ERROR(ME_DOMAIN, "");
485 qb[0] = qb[1] = ML_NAN; return;
486 }
487
488
489 /* w := Newton step size (F(.) - a) / F'(.) or,
490 * -- log: (lF - la) / (F' / F) = exp(lF) * (lF - la) / F'
491 */
492 w = log_p
493 ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta))
494 : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta));
495 if (i_pb >= n_N && w * wprev <= 0.)
496 prev = fmax2(fabs(adj),fpu);
497 R_ifDEBUG_printf(
498 "N(i=%2d): x0=%#17.15g, pb(x0)=%#15.9g, w=%#15.9g, %s prev=%g,",
499 i_pb, xinbta, y, w,
500 (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev);
501 g = 1;
502 for (i_inn=0; i_inn < 1000;i_inn++) {
503 adj = g * w;
504 // take full Newton steps at the beginning; only then safe guard:
505 if (i_pb < n_N || fabs(adj) < prev) {
506 tx = xinbta - adj; // x_{n+1} = x_n - g*w
507 if (0. <= tx && tx <= 1.) {
508 if (prev <= acu || fabs(w) <= acu) {
509 R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g, %s <= acu ==> convergence\n",
510 i_inn, -adj, (prev <= acu) ? "prev" : "|w|");
511 goto L_converged;
512 }
513 if (tx != 0. && tx != 1)
514 break;
515 }
516 }
517 g /= 3;
518 }
519 R_ifDEBUG_printf(" it{in}=%d, delta(x)=%g\n", i_inn, tx - xinbta);
520 if (fabs(tx - xinbta) <= 4e-16 * (tx + xinbta)) // "<=" : (.) == 0
521 goto L_converged;
522 xinbta = tx;
523 if(tx == 0) // "we have lost"
524 break;
525 wprev = w;
526 } // for( i_pb ..)
527
528 } // end{else : normal scale Newton}
529
530 /*-- NOT converged: Iteration count --*/
531 warned = TRUE;
532 ML_ERROR(ME_PRECISION, "qbeta");
533
534 L_converged:
535 log_ = log_p || use_log_x; // only for printing
536 R_ifDEBUG_printf(" %s: Final delta(y) = %g%s\n",
537 warned ? "_NO_ convergence" : "converged",
538 y - (log_ ? la : a), (log_ ? " (log_)" : ""));
539 if((log_ && y == ML_NEGINF) || (!log_ && y == 0)) {
540 // stuck at left, try if smallest positive number is "better"
541 w = pbeta_raw(DBL_very_MIN, pp, qq, TRUE, log_);
542 if(log_ || fabs(w - a) <= fabs(y - a)) {
543 tx = DBL_very_MIN;
544 u_n = DBL_log_v_MIN;// = log(DBL_very_MIN)
545 }
546 add_N_step = FALSE; // not trying to do better anymore
547 }
548 else if(!warned && (log_ ? fabs(y - la) > 3 : fabs(y - a) > 1e-4)) {
549 if(!(log_ && y == ML_NEGINF &&
550 // e.g. qbeta(-1e-10, .2, .03, log=TRUE) cannot get accurate ==> do NOT warn
551 pbeta_raw(DBL_1__eps, // = 1 - eps
552 pp, qq, TRUE, TRUE) > la + 2))
553 MATHLIB_WARNING2( // low accuracy for more platform independent output:
554 "qbeta(a, *) =: x0 with |pbeta(x0,*%s) - alpha| = %.5g is not accurate",
555 (log_ ? ", log_" : ""), fabs(y - (log_ ? la : a)));
556 }
557 L_return:
558 if(give_log_q) { // ==> use_log_x , too
559 if(!use_log_x) // (see if claim above is true)
560 MATHLIB_WARNING(
561 "qbeta() L_return, u_n=%g; give_log_q=TRUE but use_log_x=FALSE -- please report!",
562 u_n);
563 double r = R_Log1_Exp(u_n);
564 if(swap_tail) {
565 qb[0] = r; qb[1] = u_n;
566 } else {
567 qb[0] = u_n; qb[1] = r;
568 }
569 } else {
570 if(use_log_x) {
571 if(add_N_step) {
572 /* add one last Newton step on original x scale, e.g., for
573 qbeta(2^-98, 0.125, 2^-96) */
574 xinbta = exp(u_n);
575 y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p);
576 w = log_p
577 ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta))
578 : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta));
579 tx = xinbta - w;
580 R_ifDEBUG_printf(" Final Newton correction(non-log scale):\n"
581 // \n xinbta=%.16g
582 " xinbta=%.16g, y=%g, w=-Delta(x)=%g. \n=> new x=%.16g\n",
583 xinbta, y, w, tx);
584 } else {
585 if(swap_tail) {
586 qb[0] = -expm1(u_n); qb[1] = exp (u_n);
587 } else {
588 qb[0] = exp (u_n); qb[1] = -expm1(u_n);
589 }
590 return;
591 }
592 }
593 if(swap_tail) {
594 qb[0] = 1 - tx; qb[1] = tx;
595 } else {
596 qb[0] = tx; qb[1] = 1 - tx;
597 }
598 }
599 return;
600 }
601