From be5aaa48f9c40361e4f4235bf39b0826f3633197 Mon Sep 17 00:00:00 2001 From: Simon Byrne Date: Wed, 25 Jan 2023 10:33:48 -0800 Subject: [PATCH] update to 4.2.2 (#43) --- RVERSION | 2 +- include/R_ext/RS.h | 37 ++--- include/R_ext/Random.h | 4 +- include/Rmath.h | 24 +-- src/bd0.c | 321 ++++++++++++++++++++++++++++++++++++++--- src/cospi.c | 27 ++-- src/dbinom.c | 7 +- src/dnbinom.c | 40 +++-- src/dnorm.c | 2 +- src/dnt.c | 2 +- src/dpois.c | 18 ++- src/dpq.h | 4 +- src/fprec.c | 10 +- src/fround.c | 6 +- src/i1mach.c | 2 +- src/lgamma.c | 6 +- src/lgammacor.c | 4 +- src/log1p.c | 2 +- src/nmath.h | 8 +- src/pbeta.c | 11 +- src/pgamma.c | 16 +- src/plogis.c | 9 +- src/pnbeta.c | 10 +- src/pnchisq.c | 28 ++-- src/pnorm.c | 39 ++--- src/polygamma.c | 121 +++++++++------- src/qDiscrete_search.h | 219 ++++++++++++++++++++++++++++ src/qbeta.c | 206 ++++++++++++++++---------- src/qbinom.c | 102 ++++--------- src/qf.c | 6 +- src/qnbinom.c | 97 ++++--------- src/qnbinom_mu.c | 83 +++++++++++ src/qnchisq.c | 11 +- src/qnorm.c | 38 +++-- src/qpois.c | 79 +++------- src/qt.c | 27 +++- src/rmultinom.c | 2 +- src/stirlerr.c | 7 +- 38 files changed, 1107 insertions(+), 530 deletions(-) create mode 100644 src/qDiscrete_search.h create mode 100644 src/qnbinom_mu.c diff --git a/RVERSION b/RVERSION index c5106e6..af8c8ec 100644 --- a/RVERSION +++ b/RVERSION @@ -1 +1 @@ -4.0.4 +4.2.2 diff --git a/include/R_ext/RS.h b/include/R_ext/RS.h index 4fc82ff..50f3c8b 100644 --- a/include/R_ext/RS.h +++ b/include/R_ext/RS.h @@ -1,6 +1,6 @@ /* * R : A Computer Language for Statistical Data Analysis - * Copyright (C) 1999-2016 The R Core Team. + * Copyright (C) 1999-2022 The R Core Team. * * This header file is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -21,7 +21,7 @@ * https://www.R-project.org/Licenses/ */ -/* Included by R.h: API */ +/* Included by R.h: mainly API */ #ifndef R_RS_H #define R_RS_H @@ -42,25 +42,6 @@ extern "C" { #endif -/* S Like Error Handling */ - -#include /* for error and warning */ - -#ifndef STRICT_R_HEADERS - -#define R_PROBLEM_BUFSIZE 4096 -/* Parentheses added for FC4 with gcc4 and -D_FORTIFY_SOURCE=2 */ -#define PROBLEM {char R_problem_buf[R_PROBLEM_BUFSIZE];(sprintf)(R_problem_buf, -#define MESSAGE {char R_problem_buf[R_PROBLEM_BUFSIZE];(sprintf)(R_problem_buf, -#define ERROR ),error(R_problem_buf);} -#define RECOVER(x) ),error(R_problem_buf);} -#define WARNING(x) ),warning(R_problem_buf);} -#define LOCAL_EVALUATOR /**/ -#define NULL_ENTRY /**/ -#define WARN WARNING(NULL) - -#endif - /* S Like Memory Management */ extern void *R_chk_calloc(R_SIZE_T, R_SIZE_T); @@ -68,21 +49,25 @@ extern void *R_chk_realloc(void *, R_SIZE_T); extern void R_chk_free(void *); #ifndef STRICT_R_HEADERS -/* S-PLUS 3.x but not 5.x NULLs the pointer in the following */ +/* S-PLUS 3.x but not 5.x NULLed the pointer in Free */ #define Calloc(n, t) (t *) R_chk_calloc( (R_SIZE_T) (n), sizeof(t) ) #define Realloc(p,n,t) (t *) R_chk_realloc( (void *)(p), (R_SIZE_T)((n) * sizeof(t)) ) #define Free(p) (R_chk_free( (void *)(p) ), (p) = NULL) #endif + #define R_Calloc(n, t) (t *) R_chk_calloc( (R_SIZE_T) (n), sizeof(t) ) #define R_Realloc(p,n,t) (t *) R_chk_realloc( (void *)(p), (R_SIZE_T)((n) * sizeof(t)) ) #define R_Free(p) (R_chk_free( (void *)(p) ), (p) = NULL) +/* undocumented until 4.1.2: widely used. */ #define Memcpy(p,q,n) memcpy( p, q, (R_SIZE_T)(n) * sizeof(*p) ) -/* added for 3.0.0 */ +/* added for 3.0.0 but undocumented until 4.1.2. + Used by a couple of packages. */ #define Memzero(p,n) memset(p, 0, (R_SIZE_T)(n) * sizeof(*p)) -#define CallocCharBuf(n) (char *) R_chk_calloc((R_SIZE_T) ((n)+1), sizeof(char)) +/* Added in R 2.6.0 */ +#define CallocCharBuf(n) (char *) R_chk_calloc(((R_SIZE_T)(n))+1, sizeof(char)) /* S Like Fortran Interface */ /* These may not be adequate everywhere. Convex had _ prepending common @@ -98,9 +83,11 @@ extern void R_chk_free(void *); #define F77_COM(x) F77_CALL(x) #define F77_COMDECL(x) F77_CALL(x) -#ifndef NO_CALL_R +/* Depreacated in R 2.15.0, non-API +#if !defined(NO_CALL_R) && defined(DECLARE_LEGACY_CALL_R) void call_R(char*, long, void**, char**, long*, char**, long, char**); #endif +*/ #ifdef __cplusplus } diff --git a/include/R_ext/Random.h b/include/R_ext/Random.h index c9d5d03..7dc13f9 100644 --- a/include/R_ext/Random.h +++ b/include/R_ext/Random.h @@ -1,6 +1,6 @@ /* * R : A Computer Language for Statistical Data Analysis - * Copyright (C) 1998-2019 The R Core Team + * Copyright (C) 1998-2022 The R Core Team * * This header file is free software; you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by @@ -58,7 +58,7 @@ typedef enum { ROUNDING, REJECTION } Sampletype; -Sampletype R_sample_kind(); +Sampletype R_sample_kind(void); void GetRNGstate(void); void PutRNGstate(void); diff --git a/include/Rmath.h b/include/Rmath.h index 08dd95d..609eb9f 100644 --- a/include/Rmath.h +++ b/include/Rmath.h @@ -1,6 +1,6 @@ /* -*- C -*- * Mathlib : A C Library of Special Functions - * Copyright (C) 1998-2018 The R Core Team + * Copyright (C) 1998-2022 The R Core Team * Copyright (C) 2004 The R Foundation * * This program is free software; you can redistribute it and/or modify @@ -26,6 +26,7 @@ and nothing else. It is part of the API and supports 'standalone Rmath'. + Some entries possibly are not yet documented in 'Writing R Extensions'. */ #ifndef RMATH_H @@ -35,12 +36,9 @@ #ifndef __STDC_WANT_IEC_60559_FUNCS_EXT__ # define __STDC_WANT_IEC_60559_FUNCS_EXT__ 1 #endif + #if defined(__cplusplus) && !defined(DO_NOT_USE_CXX_HEADERS) # include -// See comment in R.h -# ifdef __SUNPRO_CC -using namespace std; -# endif #else # include #endif @@ -254,6 +252,7 @@ double Rlog1p(double); #define lgammafn Rf_lgammafn #define lgammafn_sign Rf_lgammafn_sign #define lgamma1p Rf_lgamma1p +#define log1mexp Rf_log1mexp #define log1pexp Rf_log1pexp #define log1pmx Rf_log1pmx #define logspace_add Rf_logspace_add @@ -390,9 +389,11 @@ double pgamma(double, double, double, int, int); double qgamma(double, double, double, int, int); double rgamma(double, double); -double log1pmx(double); +double log1pmx(double); /* Accurate log(1+x) - x, {care for small x} */ double log1pexp(double); // <-- ../nmath/plogis.c -double lgamma1p(double); +double log1mexp(double); +double lgamma1p(double);/* accurate log(gamma(x+1)), small x (0 < x < 0.5) */ + double logspace_add(double, double); double logspace_sub(double, double); double logspace_sum(const double *, int); @@ -447,7 +448,7 @@ double pbinom(double, double, double, int, int); double qbinom(double, double, double, int, int); double rbinom(double, double); - /* Multnomial Distribution */ + /* Multinomial Distribution */ void rmultinom(int, double*, int, int*); @@ -543,13 +544,14 @@ double dwilcox(double, double, double, int); double pwilcox(double, double, double, int, int); double qwilcox(double, double, double, int, int); double rwilcox(double, double); - +void wilcox_free(void); /* Wilcoxon Signed Rank Distribution */ double dsignrank(double, double, int); double psignrank(double, double, int, int); double qsignrank(double, double, int, int); double rsignrank(double); +void signrank_free(void); /* Gamma and Related Functions */ double gammafn(double); @@ -592,9 +594,6 @@ double fround(double, double); double fsign(double, double); double ftrunc(double); -double log1pmx(double); /* Accurate log(1+x) - x, {care for small x} */ -double lgamma1p(double);/* accurate log(gamma(x+1)), small x (0 < x < 0.5) */ - /* More accurate cos(pi*x), sin(pi*x), tan(pi*x) These declarations might clash with system headers if someone had @@ -607,6 +606,7 @@ double cospi(double); double sinpi(double); double tanpi(double); #endif +double Rtanpi(double); /* our own in any case */ /* Compute the log of a sum or difference from logs of terms, i.e., * diff --git a/src/bd0.c b/src/bd0.c index 713b5b5..5b9d36d 100644 --- a/src/bd0.c +++ b/src/bd0.c @@ -1,10 +1,11 @@ /* - * AUTHOR - * Catherine Loader, catherine@research.bell-labs.com. - * October 23, 2000. + * AUTHORS + * Catherine Loader, catherine@research.bell-labs.com, October 23, 2000. [ bd0() ] + * Morten Welinder, see Bugzilla PR#15628, 2014 [ebd0() ] * - * Merge in to R: - * Copyright (C) 2000-2014 The R Core Team + * Merge in to R (and much more): + + * Copyright (C) 2000-2022 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -31,31 +32,315 @@ * for all x and M=np. In particular for x/np close to 1, direct * evaluation fails, and evaluation is based on the Taylor series * of log((1+v)/(1-v)) with v = (x-M)/(x+M) = (x-np)/(x+np). + * + * Martyn Plummer had the nice idea to use log1p() and Martin Maechler + * emphasized the extra need to control cancellation. + * + * MP: t := (x-M)/M ( <==> 1+t = x/M ==> + * + * bd0 = M*[ x/M * log(x/M) + 1 - (x/M) ] = M*[ (1+t)*log1p(t) + 1 - (1+t) ] + * = M*[ (1+t)*log1p(t) - t ] =: M * p1log1pm(t) =: M * p1l1(t) + * MM: The above is very nice, as the "simple" p1l1() function would be useful + * to have available in a fast numerical stable way more generally. */ #include "nmath.h" double attribute_hidden bd0(double x, double np) { - double ej, s, s1, v; - int j; - if(!R_FINITE(x) || !R_FINITE(np) || np == 0.0) ML_WARN_return_NAN; if (fabs(x-np) < 0.1*(x+np)) { - v = (x-np)/(x+np); // might underflow to 0 - s = (x-np)*v;/* s using v -- change by MM */ + double + v = (x-np)/(x+np), // might underflow to 0 + s = (x-np)*v; if(fabs(s) < DBL_MIN) return s; - ej = 2*x*v; - v = v*v; - for (j = 1; j < 1000; j++) { /* Taylor series; 1000: no infinite loop + double ej = 2*x*v; + v *= v; // "v = v^2" + for (int j = 1; j < 1000; j++) { /* Taylor series; 1000: no infinite loop as |v| < .1, v^2000 is "zero" */ - ej *= v;// = v^(2j+1) - s1 = s+ej/((j<<1)+1); - if (s1 == s) /* last term was effectively 0 */ - return s1 ; - s = s1; + ej *= v;// = 2 x v^(2j+1) + double s_ = s; + s += ej/((j<<1)+1); + if (s == s_) { /* last term was effectively 0 */ +#ifdef DEBUG_bd0 + REprintf("bd0(%g, %g): T.series w/ %d terms -> bd0=%g\n", x, np, j, s); +#endif + return s; + } } + MATHLIB_WARNING4("bd0(%g, %g): T.series failed to converge in 1000 it.; s=%g, ej/(2j+1)=%g\n", + x, np, s, ej/((1000<<1)+1)); } /* else: | x - np | is not too small */ return(x*log(x/np)+np-x); } + + +// ebd0(): R Bugzilla PR#15628 -- proposed accuracy improvement by Morten Welinder + +/* + * A table of logs for scaling purposes. Each value has four parts with + * 23 bits in each. That means each part can be multiplied by a double + * with at most 30 bits set and not have any rounding error. Note, that + * the first entry is log(2). + * + * Entry i is associated with the value r = 0.5 + i / 256.0. The + * argument to log is p/q where q=1024 and p=floor(q / r + 0.5). + * Thus r*p/q is close to 1. + */ +static const float bd0_scale[128 + 1][4] = { + { +0x1.62e430p-1, -0x1.05c610p-29, -0x1.950d88p-54, +0x1.d9cc02p-79 }, /* 128: log(2048/1024.) */ + { +0x1.5ee02cp-1, -0x1.6dbe98p-25, -0x1.51e540p-50, +0x1.2bfa48p-74 }, /* 129: log(2032/1024.) */ + { +0x1.5ad404p-1, +0x1.86b3e4p-26, +0x1.9f6534p-50, +0x1.54be04p-74 }, /* 130: log(2016/1024.) */ + { +0x1.570124p-1, -0x1.9ed750p-25, -0x1.f37dd0p-51, +0x1.10b770p-77 }, /* 131: log(2001/1024.) */ + { +0x1.5326e4p-1, -0x1.9b9874p-25, -0x1.378194p-49, +0x1.56feb2p-74 }, /* 132: log(1986/1024.) */ + { +0x1.4f4528p-1, +0x1.aca70cp-28, +0x1.103e74p-53, +0x1.9c410ap-81 }, /* 133: log(1971/1024.) */ + { +0x1.4b5bd8p-1, -0x1.6a91d8p-25, -0x1.8e43d0p-50, -0x1.afba9ep-77 }, /* 134: log(1956/1024.) */ + { +0x1.47ae54p-1, -0x1.abb51cp-25, +0x1.19b798p-51, +0x1.45e09cp-76 }, /* 135: log(1942/1024.) */ + { +0x1.43fa00p-1, -0x1.d06318p-25, -0x1.8858d8p-49, -0x1.1927c4p-75 }, /* 136: log(1928/1024.) */ + { +0x1.3ffa40p-1, +0x1.1a427cp-25, +0x1.151640p-53, -0x1.4f5606p-77 }, /* 137: log(1913/1024.) */ + { +0x1.3c7c80p-1, -0x1.19bf48p-34, +0x1.05fc94p-58, -0x1.c096fcp-82 }, /* 138: log(1900/1024.) */ + { +0x1.38b320p-1, +0x1.6b5778p-25, +0x1.be38d0p-50, -0x1.075e96p-74 }, /* 139: log(1886/1024.) */ + { +0x1.34e288p-1, +0x1.d9ce1cp-25, +0x1.316eb8p-49, +0x1.2d885cp-73 }, /* 140: log(1872/1024.) */ + { +0x1.315124p-1, +0x1.c2fc60p-29, -0x1.4396fcp-53, +0x1.acf376p-78 }, /* 141: log(1859/1024.) */ + { +0x1.2db954p-1, +0x1.720de4p-25, -0x1.d39b04p-49, -0x1.f11176p-76 }, /* 142: log(1846/1024.) */ + { +0x1.2a1b08p-1, -0x1.562494p-25, +0x1.a7863cp-49, +0x1.85dd64p-73 }, /* 143: log(1833/1024.) */ + { +0x1.267620p-1, +0x1.3430e0p-29, -0x1.96a958p-56, +0x1.f8e636p-82 }, /* 144: log(1820/1024.) */ + { +0x1.23130cp-1, +0x1.7bebf4p-25, +0x1.416f1cp-52, -0x1.78dd36p-77 }, /* 145: log(1808/1024.) */ + { +0x1.1faa34p-1, +0x1.70e128p-26, +0x1.81817cp-50, -0x1.c2179cp-76 }, /* 146: log(1796/1024.) */ + { +0x1.1bf204p-1, +0x1.3a9620p-28, +0x1.2f94c0p-52, +0x1.9096c0p-76 }, /* 147: log(1783/1024.) */ + { +0x1.187ce4p-1, -0x1.077870p-27, +0x1.655a80p-51, +0x1.eaafd6p-78 }, /* 148: log(1771/1024.) */ + { +0x1.1501c0p-1, -0x1.406cacp-25, -0x1.e72290p-49, +0x1.5dd800p-73 }, /* 149: log(1759/1024.) */ + { +0x1.11cb80p-1, +0x1.787cd0p-25, -0x1.efdc78p-51, -0x1.5380cep-77 }, /* 150: log(1748/1024.) */ + { +0x1.0e4498p-1, +0x1.747324p-27, -0x1.024548p-51, +0x1.77a5a6p-75 }, /* 151: log(1736/1024.) */ + { +0x1.0b036cp-1, +0x1.690c74p-25, +0x1.5d0cc4p-50, -0x1.c0e23cp-76 }, /* 152: log(1725/1024.) */ + { +0x1.077070p-1, -0x1.a769bcp-27, +0x1.452234p-52, +0x1.6ba668p-76 }, /* 153: log(1713/1024.) */ + { +0x1.04240cp-1, -0x1.a686acp-27, -0x1.ef46b0p-52, -0x1.5ce10cp-76 }, /* 154: log(1702/1024.) */ + { +0x1.00d22cp-1, +0x1.fc0e10p-25, +0x1.6ee034p-50, -0x1.19a2ccp-74 }, /* 155: log(1691/1024.) */ + { +0x1.faf588p-2, +0x1.ef1e64p-27, -0x1.26504cp-54, -0x1.b15792p-82 }, /* 156: log(1680/1024.) */ + { +0x1.f4d87cp-2, +0x1.d7b980p-26, -0x1.a114d8p-50, +0x1.9758c6p-75 }, /* 157: log(1670/1024.) */ + { +0x1.ee1414p-2, +0x1.2ec060p-26, +0x1.dc00fcp-52, +0x1.f8833cp-76 }, /* 158: log(1659/1024.) */ + { +0x1.e7e32cp-2, -0x1.ac796cp-27, -0x1.a68818p-54, +0x1.235d02p-78 }, /* 159: log(1649/1024.) */ + { +0x1.e108a0p-2, -0x1.768ba4p-28, -0x1.f050a8p-52, +0x1.00d632p-82 }, /* 160: log(1638/1024.) */ + { +0x1.dac354p-2, -0x1.d3a6acp-30, +0x1.18734cp-57, -0x1.f97902p-83 }, /* 161: log(1628/1024.) */ + { +0x1.d47424p-2, +0x1.7dbbacp-31, -0x1.d5ada4p-56, +0x1.56fcaap-81 }, /* 162: log(1618/1024.) */ + { +0x1.ce1af0p-2, +0x1.70be7cp-27, +0x1.6f6fa4p-51, +0x1.7955a2p-75 }, /* 163: log(1608/1024.) */ + { +0x1.c7b798p-2, +0x1.ec36ecp-26, -0x1.07e294p-50, -0x1.ca183cp-75 }, /* 164: log(1598/1024.) */ + { +0x1.c1ef04p-2, +0x1.c1dfd4p-26, +0x1.888eecp-50, -0x1.fd6b86p-75 }, /* 165: log(1589/1024.) */ + { +0x1.bb7810p-2, +0x1.478bfcp-26, +0x1.245b8cp-50, +0x1.ea9d52p-74 }, /* 166: log(1579/1024.) */ + { +0x1.b59da0p-2, -0x1.882b08p-27, +0x1.31573cp-53, -0x1.8c249ap-77 }, /* 167: log(1570/1024.) */ + { +0x1.af1294p-2, -0x1.b710f4p-27, +0x1.622670p-51, +0x1.128578p-76 }, /* 168: log(1560/1024.) */ + { +0x1.a925d4p-2, -0x1.0ae750p-27, +0x1.574ed4p-51, +0x1.084996p-75 }, /* 169: log(1551/1024.) */ + { +0x1.a33040p-2, +0x1.027d30p-29, +0x1.b9a550p-53, -0x1.b2e38ap-78 }, /* 170: log(1542/1024.) */ + { +0x1.9d31c0p-2, -0x1.5ec12cp-26, -0x1.5245e0p-52, +0x1.2522d0p-79 }, /* 171: log(1533/1024.) */ + { +0x1.972a34p-2, +0x1.135158p-30, +0x1.a5c09cp-56, +0x1.24b70ep-80 }, /* 172: log(1524/1024.) */ + { +0x1.911984p-2, +0x1.0995d4p-26, +0x1.3bfb5cp-50, +0x1.2c9dd6p-75 }, /* 173: log(1515/1024.) */ + { +0x1.8bad98p-2, -0x1.1d6144p-29, +0x1.5b9208p-53, +0x1.1ec158p-77 }, /* 174: log(1507/1024.) */ + { +0x1.858b58p-2, -0x1.1b4678p-27, +0x1.56cab4p-53, -0x1.2fdc0cp-78 }, /* 175: log(1498/1024.) */ + { +0x1.7f5fa0p-2, +0x1.3aaf48p-27, +0x1.461964p-51, +0x1.4ae476p-75 }, /* 176: log(1489/1024.) */ + { +0x1.79db68p-2, -0x1.7e5054p-26, +0x1.673750p-51, -0x1.a11f7ap-76 }, /* 177: log(1481/1024.) */ + { +0x1.744f88p-2, -0x1.cc0e18p-26, -0x1.1e9d18p-50, -0x1.6c06bcp-78 }, /* 178: log(1473/1024.) */ + { +0x1.6e08ecp-2, -0x1.5d45e0p-26, -0x1.c73ec8p-50, +0x1.318d72p-74 }, /* 179: log(1464/1024.) */ + { +0x1.686c80p-2, +0x1.e9b14cp-26, -0x1.13bbd4p-50, -0x1.efeb1cp-78 }, /* 180: log(1456/1024.) */ + { +0x1.62c830p-2, -0x1.a8c70cp-27, -0x1.5a1214p-51, -0x1.bab3fcp-79 }, /* 181: log(1448/1024.) */ + { +0x1.5d1bdcp-2, -0x1.4fec6cp-31, +0x1.423638p-56, +0x1.ee3feep-83 }, /* 182: log(1440/1024.) */ + { +0x1.576770p-2, +0x1.7455a8p-26, -0x1.3ab654p-50, -0x1.26be4cp-75 }, /* 183: log(1432/1024.) */ + { +0x1.5262e0p-2, -0x1.146778p-26, -0x1.b9f708p-52, -0x1.294018p-77 }, /* 184: log(1425/1024.) */ + { +0x1.4c9f08p-2, +0x1.e152c4p-26, -0x1.dde710p-53, +0x1.fd2208p-77 }, /* 185: log(1417/1024.) */ + { +0x1.46d2d8p-2, +0x1.c28058p-26, -0x1.936284p-50, +0x1.9fdd68p-74 }, /* 186: log(1409/1024.) */ + { +0x1.41b940p-2, +0x1.cce0c0p-26, -0x1.1a4050p-50, +0x1.bc0376p-76 }, /* 187: log(1402/1024.) */ + { +0x1.3bdd24p-2, +0x1.d6296cp-27, +0x1.425b48p-51, -0x1.cddb2cp-77 }, /* 188: log(1394/1024.) */ + { +0x1.36b578p-2, -0x1.287ddcp-27, -0x1.2d0f4cp-51, +0x1.38447ep-75 }, /* 189: log(1387/1024.) */ + { +0x1.31871cp-2, +0x1.2a8830p-27, +0x1.3eae54p-52, -0x1.898136p-77 }, /* 190: log(1380/1024.) */ + { +0x1.2b9304p-2, -0x1.51d8b8p-28, +0x1.27694cp-52, -0x1.fd852ap-76 }, /* 191: log(1372/1024.) */ + { +0x1.265620p-2, -0x1.d98f3cp-27, +0x1.a44338p-51, -0x1.56e85ep-78 }, /* 192: log(1365/1024.) */ + { +0x1.211254p-2, +0x1.986160p-26, +0x1.73c5d0p-51, +0x1.4a861ep-75 }, /* 193: log(1358/1024.) */ + { +0x1.1bc794p-2, +0x1.fa3918p-27, +0x1.879c5cp-51, +0x1.16107cp-78 }, /* 194: log(1351/1024.) */ + { +0x1.1675ccp-2, -0x1.4545a0p-26, +0x1.c07398p-51, +0x1.f55c42p-76 }, /* 195: log(1344/1024.) */ + { +0x1.111ce4p-2, +0x1.f72670p-37, -0x1.b84b5cp-61, +0x1.a4a4dcp-85 }, /* 196: log(1337/1024.) */ + { +0x1.0c81d4p-2, +0x1.0c150cp-27, +0x1.218600p-51, -0x1.d17312p-76 }, /* 197: log(1331/1024.) */ + { +0x1.071b84p-2, +0x1.fcd590p-26, +0x1.a3a2e0p-51, +0x1.fe5ef8p-76 }, /* 198: log(1324/1024.) */ + { +0x1.01ade4p-2, -0x1.bb1844p-28, +0x1.db3cccp-52, +0x1.1f56fcp-77 }, /* 199: log(1317/1024.) */ + { +0x1.fa01c4p-3, -0x1.12a0d0p-29, -0x1.f71fb0p-54, +0x1.e287a4p-78 }, /* 200: log(1311/1024.) */ + { +0x1.ef0adcp-3, +0x1.7b8b28p-28, -0x1.35bce4p-52, -0x1.abc8f8p-79 }, /* 201: log(1304/1024.) */ + { +0x1.e598ecp-3, +0x1.5a87e4p-27, -0x1.134bd0p-51, +0x1.c2cebep-76 }, /* 202: log(1298/1024.) */ + { +0x1.da85d8p-3, -0x1.df31b0p-27, +0x1.94c16cp-57, +0x1.8fd7eap-82 }, /* 203: log(1291/1024.) */ + { +0x1.d0fb80p-3, -0x1.bb5434p-28, -0x1.ea5640p-52, -0x1.8ceca4p-77 }, /* 204: log(1285/1024.) */ + { +0x1.c765b8p-3, +0x1.e4d68cp-27, +0x1.5b59b4p-51, +0x1.76f6c4p-76 }, /* 205: log(1279/1024.) */ + { +0x1.bdc46cp-3, -0x1.1cbb50p-27, +0x1.2da010p-51, +0x1.eb282cp-75 }, /* 206: log(1273/1024.) */ + { +0x1.b27980p-3, -0x1.1b9ce0p-27, +0x1.7756f8p-52, +0x1.2ff572p-76 }, /* 207: log(1266/1024.) */ + { +0x1.a8bed0p-3, -0x1.bbe874p-30, +0x1.85cf20p-56, +0x1.b9cf18p-80 }, /* 208: log(1260/1024.) */ + { +0x1.9ef83cp-3, +0x1.2769a4p-27, -0x1.85bda0p-52, +0x1.8c8018p-79 }, /* 209: log(1254/1024.) */ + { +0x1.9525a8p-3, +0x1.cf456cp-27, -0x1.7137d8p-52, -0x1.f158e8p-76 }, /* 210: log(1248/1024.) */ + { +0x1.8b46f8p-3, +0x1.11b12cp-30, +0x1.9f2104p-54, -0x1.22836ep-78 }, /* 211: log(1242/1024.) */ + { +0x1.83040cp-3, +0x1.2379e4p-28, +0x1.b71c70p-52, -0x1.990cdep-76 }, /* 212: log(1237/1024.) */ + { +0x1.790ed4p-3, +0x1.dc4c68p-28, -0x1.910ac8p-52, +0x1.dd1bd6p-76 }, /* 213: log(1231/1024.) */ + { +0x1.6f0d28p-3, +0x1.5cad68p-28, +0x1.737c94p-52, -0x1.9184bap-77 }, /* 214: log(1225/1024.) */ + { +0x1.64fee8p-3, +0x1.04bf88p-28, +0x1.6fca28p-52, +0x1.8884a8p-76 }, /* 215: log(1219/1024.) */ + { +0x1.5c9400p-3, +0x1.d65cb0p-29, -0x1.b2919cp-53, +0x1.b99bcep-77 }, /* 216: log(1214/1024.) */ + { +0x1.526e60p-3, -0x1.c5e4bcp-27, -0x1.0ba380p-52, +0x1.d6e3ccp-79 }, /* 217: log(1208/1024.) */ + { +0x1.483bccp-3, +0x1.9cdc7cp-28, -0x1.5ad8dcp-54, -0x1.392d3cp-83 }, /* 218: log(1202/1024.) */ + { +0x1.3fb25cp-3, -0x1.a6ad74p-27, +0x1.5be6b4p-52, -0x1.4e0114p-77 }, /* 219: log(1197/1024.) */ + { +0x1.371fc4p-3, -0x1.fe1708p-27, -0x1.78864cp-52, -0x1.27543ap-76 }, /* 220: log(1192/1024.) */ + { +0x1.2cca10p-3, -0x1.4141b4p-28, -0x1.ef191cp-52, +0x1.00ee08p-76 }, /* 221: log(1186/1024.) */ + { +0x1.242310p-3, +0x1.3ba510p-27, -0x1.d003c8p-51, +0x1.162640p-76 }, /* 222: log(1181/1024.) */ + { +0x1.1b72acp-3, +0x1.52f67cp-27, -0x1.fd6fa0p-51, +0x1.1a3966p-77 }, /* 223: log(1176/1024.) */ + { +0x1.10f8e4p-3, +0x1.129cd8p-30, +0x1.31ef30p-55, +0x1.a73e38p-79 }, /* 224: log(1170/1024.) */ + { +0x1.08338cp-3, -0x1.005d7cp-27, -0x1.661a9cp-51, +0x1.1f138ap-79 }, /* 225: log(1165/1024.) */ + { +0x1.fec914p-4, -0x1.c482a8p-29, -0x1.55746cp-54, +0x1.99f932p-80 }, /* 226: log(1160/1024.) */ + { +0x1.ed1794p-4, +0x1.d06f00p-29, +0x1.75e45cp-53, -0x1.d0483ep-78 }, /* 227: log(1155/1024.) */ + { +0x1.db5270p-4, +0x1.87d928p-32, -0x1.0f52a4p-57, +0x1.81f4a6p-84 }, /* 228: log(1150/1024.) */ + { +0x1.c97978p-4, +0x1.af1d24p-29, -0x1.0977d0p-60, -0x1.8839d0p-84 }, /* 229: log(1145/1024.) */ + { +0x1.b78c84p-4, -0x1.44f124p-28, -0x1.ef7bc4p-52, +0x1.9e0650p-78 }, /* 230: log(1140/1024.) */ + { +0x1.a58b60p-4, +0x1.856464p-29, +0x1.c651d0p-55, +0x1.b06b0cp-79 }, /* 231: log(1135/1024.) */ + { +0x1.9375e4p-4, +0x1.5595ecp-28, +0x1.dc3738p-52, +0x1.86c89ap-81 }, /* 232: log(1130/1024.) */ + { +0x1.814be4p-4, -0x1.c073fcp-28, -0x1.371f88p-53, -0x1.5f4080p-77 }, /* 233: log(1125/1024.) */ + { +0x1.6f0d28p-4, +0x1.5cad68p-29, +0x1.737c94p-53, -0x1.9184bap-78 }, /* 234: log(1120/1024.) */ + { +0x1.60658cp-4, -0x1.6c8af4p-28, +0x1.d8ef74p-55, +0x1.c4f792p-80 }, /* 235: log(1116/1024.) */ + { +0x1.4e0110p-4, +0x1.146b5cp-29, +0x1.73f7ccp-54, -0x1.d28db8p-79 }, /* 236: log(1111/1024.) */ + { +0x1.3b8758p-4, +0x1.8b1b70p-28, -0x1.20aca4p-52, -0x1.651894p-76 }, /* 237: log(1106/1024.) */ + { +0x1.28f834p-4, +0x1.43b6a4p-30, -0x1.452af8p-55, +0x1.976892p-80 }, /* 238: log(1101/1024.) */ + { +0x1.1a0fbcp-4, -0x1.e4075cp-28, +0x1.1fe618p-52, +0x1.9d6dc2p-77 }, /* 239: log(1097/1024.) */ + { +0x1.075984p-4, -0x1.4ce370p-29, -0x1.d9fc98p-53, +0x1.4ccf12p-77 }, /* 240: log(1092/1024.) */ + { +0x1.f0a30cp-5, +0x1.162a68p-37, -0x1.e83368p-61, -0x1.d222a6p-86 }, /* 241: log(1088/1024.) */ + { +0x1.cae730p-5, -0x1.1a8f7cp-31, -0x1.5f9014p-55, +0x1.2720c0p-79 }, /* 242: log(1083/1024.) */ + { +0x1.ac9724p-5, -0x1.e8ee08p-29, +0x1.a7de04p-54, -0x1.9bba74p-78 }, /* 243: log(1079/1024.) */ + { +0x1.868a84p-5, -0x1.ef8128p-30, +0x1.dc5eccp-54, -0x1.58d250p-79 }, /* 244: log(1074/1024.) */ + { +0x1.67f950p-5, -0x1.ed684cp-30, -0x1.f060c0p-55, -0x1.b1294cp-80 }, /* 245: log(1070/1024.) */ + { +0x1.494accp-5, +0x1.a6c890p-32, -0x1.c3ad48p-56, -0x1.6dc66cp-84 }, /* 246: log(1066/1024.) */ + { +0x1.22c71cp-5, -0x1.8abe2cp-32, -0x1.7e7078p-56, -0x1.ddc3dcp-86 }, /* 247: log(1061/1024.) */ + { +0x1.03d5d8p-5, +0x1.79cfbcp-31, -0x1.da7c4cp-58, +0x1.4e7582p-83 }, /* 248: log(1057/1024.) */ + { +0x1.c98d18p-6, +0x1.a01904p-31, -0x1.854164p-55, +0x1.883c36p-79 }, /* 249: log(1053/1024.) */ + { +0x1.8b31fcp-6, -0x1.356500p-30, +0x1.c3ab48p-55, +0x1.b69bdap-80 }, /* 250: log(1049/1024.) */ + { +0x1.3cea44p-6, +0x1.a352bcp-33, -0x1.8865acp-57, -0x1.48159cp-81 }, /* 251: log(1044/1024.) */ + { +0x1.fc0a8cp-7, -0x1.e07f84p-32, +0x1.e7cf6cp-58, +0x1.3a69c0p-82 }, /* 252: log(1040/1024.) */ + { +0x1.7dc474p-7, +0x1.f810a8p-31, -0x1.245b5cp-56, -0x1.a1f4f8p-80 }, /* 253: log(1036/1024.) */ + { +0x1.fe02a8p-8, -0x1.4ef988p-32, +0x1.1f86ecp-57, +0x1.20723cp-81 }, /* 254: log(1032/1024.) */ + { +0x1.ff00acp-9, -0x1.d4ef44p-33, +0x1.2821acp-63, +0x1.5a6d32p-87 }, /* 255: log(1028/1024.) */ + { 0, 0, 0, 0 } /* log(1024/1024) = log(1) = 0 */ +}; + + +/* + * Compute x * log (x / M) + (M - x) + * aka -x * log1pmx ((M - x) / x) + * + * Deliver the result back in two parts, *yh and *yl. + */ +void attribute_hidden ebd0(double x, double M, double *yh, double *yl) +{ + const int Sb = 10; + const double S = 1u << Sb; // = 2^10 = 1024 + const int N = 128; // == ? == G_N_ELEMENTS(bd0_scale) - 1; <<<< FIXME: + + *yl = *yh = 0; + + if (x == M) return; + if (x == 0) { *yh = M; return; } + if (M == 0) { *yh = ML_POSINF; return; } + + if (M/x == ML_POSINF) { *yh = M; return; }// as when (x == 0) + + int e; + // NB: M/x overflow handled above; underflow should be handled by fg = Inf + double r = frexp (M / x, &e); // => r in [0.5, 1) and 'e' (int) such that M/x = r * 2^e + + // prevent later overflow + if (M_LN2 * ((double) -e) > 1. + DBL_MAX / x) { *yh = ML_POSINF; return; } + + int i = (int) floor ((r - 0.5) * (2 * N) + 0.5); + // now, 0 <= i <= N + double f = floor (S / (0.5 + i / (2.0 * N)) + 0.5); + double fg = ldexp (f, -(e + Sb)); // ldexp(f, E) := f * 2^E +#ifdef DEBUG_bd0 + REprintf("ebd0(x=%g, M=%g): M/x = (r=%.15g) * 2^(e=%d); i=%d,\n f=%g, fg=f*2^-(e+%d)=%g\n", + x, M, r,e, i, f, Sb, fg); + if (fg == ML_POSINF) { + REprintf(" --> fg = +Inf --> return( +Inf )\n"); + *yh = fg; return; + } + REprintf(" bd0_sc[0][0..3]= ("); for(int j=0; j < 4; j++) REprintf("%g ", bd0_scale[0][j]); REprintf(")\n"); + REprintf("i -> bd0_sc[i][0..3]= ("); for(int j=0; j < 4; j++) REprintf("%g ", bd0_scale[i][j]); REprintf(")\n"); + REprintf( " small(?) (M*fg-x)/x = (M*fg)/x - 1 = %.16g\n", (M*fg-x)/x); +#else + if (fg == ML_POSINF) { + *yh = fg; return; + } +#endif + /* We now have (M * fg / x) close to 1. */ + + /* + * We need to compute this: + * (x/M)^x * exp(M-x) = + * (M/x)^-x * exp(M-x) = + * (M*fg/x)^-x * (fg)^x * exp(M-x) = + * (M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg) + * + * In log terms: + * log((x/M)^x * exp(M-x)) = + * log((M*fg/x)^-x * (fg)^x * exp(M*fg-x) * exp(M-M*fg)) = + * log((M*fg/x)^-x * exp(M*fg-x)) + x*log(fg) + (M-M*fg) = + * -x*log1pmx((M*fg-x)/x) + x*log(fg) + M - M*fg = + * + * Note, that fg has at most 10 bits. If M and x are suitably + * "nice" -- such as being integers or half-integers -- then + * we can compute M*fg as well as x * bd0_scale[.][.] without + * rounding errors. + */ + +#define ADD1(d_) do { \ + volatile double d = (d_); \ + double d1 = floor (d + 0.5); \ + double d2 = d - d1;/* in [-.5,.5) */ \ + *yh += d1; \ + *yl += d2; \ + } while(0) + +#ifdef DEBUG_bd0 + { + double log1__ = log1pmx((M * fg - x) / x), + xl = -x * log1__; + REprintf(" 1a. before adding -x * log1pmx(.) = -x * %g = %g\n", log1__, xl); + ADD1(xl); + REprintf(" 1. after A.(-x*l..): yl,yh = (%13g, %13g); yl+yh= %g\n", + *yl, *yh, (*yl)+(*yh)); + } + if(fg == 1) { + REprintf("___ fg = 1 ___ skipping further steps\n"); + return; + } + // else [ fg != 1 ] + REprintf(" 2: A(x*b[i,j]) and A(-x*e*b[0,j]), j=1:4:\n"); + for (int j = 0; j < 4; j++) { + ADD1( x * bd0_scale[i][j]); // handles x*log(fg*2^e) + REprintf(" j=%d: (%13g, %13g);", j, *yl, *yh); + ADD1(-x * bd0_scale[0][j] * e); // handles x*log(1/ 2^e) + REprintf(" (%13g, %13g); yl+yh= %g\n", *yl, *yh, (*yl)+(*yh)); + if(!R_FINITE(*yh)) { + REprintf(" non-finite yh --> return((yh=Inf, yl=0))\n"); + *yh = ML_POSINF; *yl = 0; return; + } + } +#else + ADD1(-x * log1pmx ((M * fg - x) / x)); + if(fg == 1) return; + // else (fg != 1) : + for (int j = 0; j < 4; j++) { + ADD1( x * bd0_scale[i][j]); // handles x*log(fg*2^e) + ADD1(-x * bd0_scale[0][j] * e); // handles x*log(1/ 2^e) + // ^^^ at end prevents overflow in ebd0(1e307, 1e300) + if(!R_FINITE(*yh)) { *yh = ML_POSINF; *yl = 0; return; } + } +#endif + + ADD1(M); +#ifdef DEBUG_bd0 + REprintf(" 3. after ADD1(M): yl,yh = (%13g, %13g); yl+yh= %g\n", *yl, *yh, (*yl)+(*yh)); +#endif + ADD1(-M * fg); +#ifdef DEBUG_bd0 + REprintf(" 4. after ADD1(- M*fg): yl,yh = (%13g, %13g); yl+yh= %g\n\n", *yl, *yh, (*yl)+(*yh)); +#endif +} + +#undef ADD1 diff --git a/src/cospi.c b/src/cospi.c index d55d895..e12c627 100644 --- a/src/cospi.c +++ b/src/cospi.c @@ -1,6 +1,6 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 2013-2016 The R Core Team + * Copyright (C) 2013-2022 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -68,13 +68,9 @@ double sinpi(double x) { } #endif -// tan(pi * x) -- exact when x = k/2 for all integer k -#if defined(HAVE_TANPI) || defined(HAVE___TANPI) -// for use in arithmetic.c, half-values documented to give NaN +// tan(pi * x) -- exact when x = k/4 for all integer k and half-values give NaN +// ----------- e.g. used in ../main/arithmetic.c : double Rtanpi(double x) -#else -double tanpi(double x) -#endif { #ifdef IEEE_754 if (ISNAN(x)) return x; @@ -82,13 +78,26 @@ double tanpi(double x) if(!R_FINITE(x)) ML_WARN_return_NAN; x = fmod(x, 1.); // tan(pi(x + k)) == tan(pi x) for all integer k - // map (-1,1) --> (-1/2, 1/2] : + // map (-1,1] --> (-1/2, 1/2] : if(x <= -0.5) x++; else if(x > 0.5) x--; - return (x == 0.) ? 0. : ((x == 0.5) ? ML_NAN : tan(M_PI * x)); + return (x == 0.) ? 0. : + ((x == 0.5 ) ? ML_NAN : + ((x == 0.25) ? 1. : + ((x == -0.25) ? -1. : + tan(M_PI * x) + ))); } +#if defined(HAVE_TANPI) || defined(HAVE___TANPI) +#else +double tanpi(double x) { + return Rtanpi(x); +} +#endif + #if !defined(HAVE_TANPI) && defined(HAVE___TANPI) double tanpi(double x) { return __tanpi(x); } +/* #else tanpi() defined from C standard math lib */ #endif diff --git a/src/dbinom.c b/src/dbinom.c index a5963db..d2f1cd9 100644 --- a/src/dbinom.c +++ b/src/dbinom.c @@ -4,7 +4,7 @@ * October 23, 2000. * * Merge in to R and further tweaks : - * Copyright (C) 2000-2015 The R Core Team + * Copyright (C) 2000-2020 The R Core Team * Copyright (C) 2008 The R Foundation * * This program is free software; you can redistribute it and/or modify @@ -43,11 +43,10 @@ double dbinom_raw(double x, double n, double p, double q, int give_log) { - double lf, lc; - if (p == 0) return((x == 0) ? R_D__1 : R_D__0); if (q == 0) return((x == n) ? R_D__1 : R_D__0); + double lc; if (x == 0) { if(n == 0) return R_D__1; lc = (p < 0.1) ? -bd0(n,n*q) - n*p : n*log(q); @@ -67,7 +66,7 @@ double dbinom_raw(double x, double n, double p, double q, int give_log) /* Upto R 2.7.1: * lf = log(M_2PI) + log(x) + log(n-x) - log(n); * -- following is much better for x << n : */ - lf = M_LN_2PI + log(x) + log1p(- x/n); + double lf = M_LN_2PI + log(x) + log1p(- x/n); return R_D_exp(lc - 0.5*lf); } diff --git a/src/dnbinom.c b/src/dnbinom.c index 2ecd831..ada687b 100644 --- a/src/dnbinom.c +++ b/src/dnbinom.c @@ -5,8 +5,8 @@ * * dnbinom_mu(): Martin Maechler, June 2008 * - * Merge in to R: - * Copyright (C) 2000--2016, The R Core Team + * Merge in to R and improvements notably for |x| << size : + * Copyright (C) 2000--2021, The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -37,8 +37,6 @@ double dnbinom(double x, double size, double prob, int give_log) { - double ans, p; - #ifdef IEEE_754 if (ISNAN(x) || ISNAN(size) || ISNAN(prob)) return x + size + prob; @@ -47,14 +45,26 @@ double dnbinom(double x, double size, double prob, int give_log) if (prob <= 0 || prob > 1 || size < 0) ML_WARN_return_NAN; R_D_nonint_check(x); if (x < 0 || !R_FINITE(x)) return R_D__0; - /* limiting case as size approaches zero is point mass at zero */ - if (x == 0 && size==0) return R_D__1; x = R_forceint(x); + if(x == 0) { + /* limiting case as size approaches zero is point mass at zero */ + if(size == 0) return R_D__1; + // size > 0: P(x, ..) = pr^n : + return(give_log ? size*log(prob) : pow(prob, size)); + } if(!R_FINITE(size)) size = DBL_MAX; - ans = dbinom_raw(size, x+size, prob, 1-prob, give_log); - p = ((double)size)/(size+x); - return((give_log) ? log(p) + ans : p * ans); + if(x < 1e-10 * size) { // instead of dbinom_raw(), use 2 terms of Abramowitz & Stegun (6.1.47) + return R_D_exp(size * log(prob) + x * (log(size) + log1p(-prob)) + - lgamma1p(x) + log1p(x*(x-1)/(2*size))); + } else { + /* log( size/(size+x) ) is much less accurate than log1p(- x/(size+x)) + for |x| << size (and actually when x < size): */ + double p = give_log ? (x < size ? log1p(-x/(size+x)) : log(size/(size+x))) + : size/(size+x), + ans = dbinom_raw(size, x+size, prob, 1-prob, give_log); + return((give_log) ? p + ans : p * ans); + } } double dnbinom_mu(double x, double size, double mu, int give_log) @@ -77,6 +87,7 @@ double dnbinom_mu(double x, double size, double mu, int give_log) */ if (x == 0 && size == 0) return R_D__1; x = R_forceint(x); + // FIXME use also for size "almost" Inf because that gives NaN ??? if(!R_FINITE(size)) // limit case: Poisson return(dpois_raw(x, mu, give_log)); @@ -85,13 +96,16 @@ double dnbinom_mu(double x, double size, double mu, int give_log) if(x < 1e-10 * size) { /* don't use dbinom_raw() but MM's formula: */ /* FIXME --- 1e-8 shows problem; rather use algdiv() from ./toms708.c */ double p = (size < mu ? log(size/(1 + size/mu)) : log(mu / (1 + mu/size))); - return R_D_exp(x * p - mu - lgamma(x+1) + + return R_D_exp(x * p - mu - lgamma1p(x) + log1p(x*(x-1)/(2*size))); } else { /* no unnecessary cancellation inside dbinom_raw, when - * x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy */ - double p = ((double)size)/(size+x), + x_ = size and n_ = x+size are so close that n_ - x_ loses accuracy + but log( size/(size+x) ) is much less accurate than log1p(- x/(size+x)) + for |x| << size (and actually when x < size): */ + double p = give_log ? (x < size ? log1p(-x/(size+x)) : log(size/(size+x))) + : size/(size+x), ans = dbinom_raw(size, x+size, size/(size+mu), mu/(size+mu), give_log); - return((give_log) ? log(p) + ans : p * ans); + return((give_log) ? p + ans : p * ans); } } diff --git a/src/dnorm.c b/src/dnorm.c index c855db8..c32b9c7 100644 --- a/src/dnorm.c +++ b/src/dnorm.c @@ -62,7 +62,7 @@ double dnorm4(double x, double mu, double sigma, int give_log) * x*x may lose upto about two digits accuracy for "large" x * Morten Welinder's proposal for PR#15620 - * https://bugs.r-project.org/bugzilla/show_bug.cgi?id=15620 + * https://bugs.r-project.org/show_bug.cgi?id=15620 * -- 1 -- No hoop jumping when we underflow to zero anyway: diff --git a/src/dnt.c b/src/dnt.c index 53f5877..0a9dcbf 100644 --- a/src/dnt.c +++ b/src/dnt.c @@ -1,6 +1,6 @@ /* * AUTHOR - * Claus Ekstrøm, ekstrom@dina.kvl.dk + * Claus Ekstrøm, ekstrom@dina.kvl.dk * July 15, 2003. * * Merge in to R: diff --git a/src/dpois.c b/src/dpois.c index 7425e6a..bc9b163 100644 --- a/src/dpois.c +++ b/src/dpois.c @@ -4,7 +4,7 @@ * October 23, 2000. * * Merge in to R: - * Copyright (C) 2000-2016 The R Core Team + * Copyright (C) 2000-2021 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -35,6 +35,10 @@ #include "nmath.h" #include "dpq.h" +#define M_SQRT_2PI 2.50662827463100050241576528481104525301 /* sqrt(2*pi) */ +// sqrt(2 * Rmpfr::Const("pi", 128)) +#define x_LRG 2.86111748575702815380240589208115399625e+307 /* = 2^1023 / pi */ + // called also from dgamma.c, pgamma.c, dnbeta.c, dnbinom.c, dnchisq.c : double dpois_raw(double x, double lambda, int give_log) { @@ -51,7 +55,17 @@ double dpois_raw(double x, double lambda, int give_log) // else return(R_D_exp(-lambda + x*log(lambda) -lgammafn(x+1))); } - return(R_D_fexp( M_2PI*x, -stirlerr(x)-bd0(x,lambda) )); + // R <= 4.0.x had return(R_D_fexp( M_2PI*x, -stirlerr(x)-bd0(x,lambda) )); + double yh, yl; + ebd0 (x, lambda, &yh, &yl); + yl += stirlerr(x); + Rboolean Lrg_x = (x >= x_LRG); //really large x <==> 2*pi*x overflows + double r = Lrg_x + ? M_SQRT_2PI * sqrt(x) // sqrt(.): avoid overflow for very large x + : M_2PI * x; + return give_log + ? -yl - yh - (Lrg_x ? log(r) : 0.5 * log(r)) + : exp(-yl) * exp(-yh) / (Lrg_x ? r : sqrt(r)); } double dpois(double x, double lambda, int give_log) diff --git a/src/dpq.h b/src/dpq.h index a03dfe1..88affab 100644 --- a/src/dpq.h +++ b/src/dpq.h @@ -1,6 +1,6 @@ /* * R : A Computer Language for Statistical Data Analysis - * Copyright (C) 2000--2015 The R Core Team + * Copyright (C) 2000--2021 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -114,6 +114,8 @@ /* additions for density functions (C.Loader) */ #define R_D_fexp(f,x) (give_log ? -0.5*log(f)+(x) : exp(x)/sqrt(f)) +// version working with rf := sqrt(f) [avoiding overflow in computation of f in the caller] +#define R_D_rtxp(rf,x) (give_log ? -log(rf)+(x) : exp(x)/(rf)) /* [neg]ative or [non int]eger : */ #define R_D_negInonint(x) (x < 0. || R_nonint(x)) diff --git a/src/fprec.c b/src/fprec.c index ebb6f67..f6b3bc2 100644 --- a/src/fprec.c +++ b/src/fprec.c @@ -1,7 +1,7 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2019 The R Core Team * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2018 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -54,13 +54,13 @@ --MM-- */ - +// R's signif(x, digits) via Math2(args, fprec) in ../main/arithmetic.c : double fprec(double x, double digits) { double l10, pow10, sgn, p10, P10; int e10, e2, do_round, dig; - /* Max.expon. of 10 (=308.2547) */ - const static int max10e = (int) (DBL_MAX_EXP * M_LOG10_2); + // Max.expon. of 10 (w/o denormalizing or overflow; = R's trunc( log10(.Machine$double.xmax) ) + const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE") if (ISNAN(x) || ISNAN(digits)) return x + digits; @@ -88,7 +88,7 @@ double fprec(double x, double digits) if(e10 > max10e) { /* numbers less than 10^(dig-1) * 1e-308 */ p10 = R_pow_di(10., e10-max10e); e10 = max10e; - } + } if(e10 > 0) { /* Try always to have pow >= 1 and so exactly representable */ pow10 = R_pow_di(10., e10); diff --git a/src/fround.c b/src/fround.c index 5aa8631..36c7c0b 100644 --- a/src/fround.c +++ b/src/fround.c @@ -32,9 +32,9 @@ #include "nmath.h" double fround(double x, double digits) { -#define MAX_DIGITS (DBL_MAX_10_EXP + DBL_DIG) - /* was DBL_MAX_10_EXP (= 308, IEEE) till R 3.6.x; before, - was (DBL_DIG - 1) till R 0.99 */ +#define MAX_DIGITS (DBL_MAX_10_EXP + DBL_DIG) /* typically = 308+15 = 323 + * was DBL_MAX_10_EXP (= 308, IEEE) till R 3.6.x; before, + * was (DBL_DIG - 1) till R 0.99 */ const static int max10e = (int) DBL_MAX_10_EXP; // == 308 ("IEEE") /* Note that large digits make sense for very small numbers */ diff --git a/src/i1mach.c b/src/i1mach.c index 5cad06c..c939c80 100644 --- a/src/i1mach.c +++ b/src/i1mach.c @@ -43,7 +43,7 @@ attribute_hidden int Rf_i1mach(int i) case 12: return FLT_MIN_EXP; case 13: return FLT_MAX_EXP; - case 14: return DBL_MANT_DIG; + case 14: return DBL_MANT_DIG; // 53 case 15: return DBL_MIN_EXP; case 16: return DBL_MAX_EXP; diff --git a/src/lgamma.c b/src/lgamma.c index 775ed75..8c06fbc 100644 --- a/src/lgamma.c +++ b/src/lgamma.c @@ -1,6 +1,6 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 2000-2018 The R Core Team + * Copyright (C) 2000-2020 The R Core Team * Copyright (C) 1998 Ross Ihaka * * This program is free software; you can redistribute it and/or modify @@ -37,6 +37,8 @@ * The accuracy of this routine compares (very) favourably * with those of the Sun Microsystems portable mathematical * library. + * + * ./toms708.c has gamln() */ #include "nmath.h" @@ -112,7 +114,7 @@ double lgammafn_sign(double x, int *sgn) if(fabs((x - trunc(x - 0.5)) * ans / x) < dxrel) { /* The answer is less than half precision because - * the argument is too near a negative integer. */ + * the argument is too near a negative integer; e.g. for lgamma(1e-7 - 11) */ ML_WARNING(ME_PRECISION, "lgamma"); } diff --git a/src/lgammacor.c b/src/lgammacor.c index 5e51b3b..f0addc5 100644 --- a/src/lgammacor.c +++ b/src/lgammacor.c @@ -1,7 +1,7 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2021 The R Core Team * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2001 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -45,7 +45,7 @@ double attribute_hidden lgammacor(double x) { - const static double algmcs[15] = { + const static double algmcs[15] = { // below, nalgm = 5 ==> only the first 5 are used! +.1666389480451863247205729650822e+0, -.1384948176067563840732986059135e-4, +.9810825646924729426157171547487e-8, diff --git a/src/log1p.c b/src/log1p.c index 38861e6..a87e5bf 100644 --- a/src/log1p.c +++ b/src/log1p.c @@ -101,7 +101,7 @@ double Rlog1p(double x) static double xmin = 0.0; if (xmin == 0.0) xmin = -1 + sqrt(DBL_EPSILON);/*was sqrt(d1mach(4)); */ - if (nlnrel == 0) /* initialize chebychev coefficients */ + if (nlnrel == 0) /* initialize chebyshev coefficients */ nlnrel = chebyshev_init(alnrcs, 43, DBL_EPSILON/20);/*was .1*d1mach(3)*/ #else # define nlnrel 22 diff --git a/src/nmath.h b/src/nmath.h index 2f86a67..5f723ac 100644 --- a/src/nmath.h +++ b/src/nmath.h @@ -1,6 +1,6 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 1998-2020 The R Core Team + * Copyright (C) 1998-2022 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -58,7 +58,7 @@ double Rf_gamma_cody(double); #else # define R_forceint(x) round(x) #endif -//R >= 3.1.0: # define R_nonint(x) (fabs((x) - R_forceint(x)) > 1e-7) +//R >= 3.1.0: previously, was defined as (fabs((x) - R_forceint(x)) > 1e-7) # define R_nonint(x) (fabs((x) - R_forceint(x)) > 1e-7*fmax2(1., fabs(x))) #ifndef MATHLIB_STANDALONE @@ -70,6 +70,7 @@ double Rf_gamma_cody(double); # define MATHLIB_WARNING3(fmt,x,x2,x3) warning(fmt,x,x2,x3) # define MATHLIB_WARNING4(fmt,x,x2,x3,x4) warning(fmt,x,x2,x3,x4) # define MATHLIB_WARNING5(fmt,x,x2,x3,x4,x5) warning(fmt,x,x2,x3,x4,x5) +# define MATHLIB_WARNING6(fmt,x,x2,x3,x4,x5,x6) warning(fmt,x,x2,x3,x4,x5,x6) #include #define ML_POSINF R_PosInf @@ -107,6 +108,7 @@ void R_CheckUserInterrupt(void); #define MATHLIB_WARNING3(fmt,x,x2,x3) printf(fmt,x,x2,x3) #define MATHLIB_WARNING4(fmt,x,x2,x3,x4) printf(fmt,x,x2,x3,x4) #define MATHLIB_WARNING5(fmt,x,x2,x3,x4,x5) printf(fmt,x,x2,x3,x4,x5) +#define MATHLIB_WARNING6(fmt,x,x2,x3,x4,x5,x6) printf(fmt,x,x2,x3,x4,x5,x6) #define ISNAN(x) (isnan(x)!=0) // Arith.h defines it @@ -187,6 +189,7 @@ int R_finite(double); /* always remap internal functions */ #define bd0 Rf_bd0 +#define ebd0 Rf_ebd0 #define chebyshev_eval Rf_chebyshev_eval #define chebyshev_init Rf_chebyshev_init #define gammalims Rf_gammalims @@ -213,6 +216,7 @@ double attribute_hidden stirlerr(double); /* Stirling expansion "error" */ double attribute_hidden lfastchoose(double, double); double attribute_hidden bd0(double, double); +void attribute_hidden ebd0(double, double, double*, double*); double attribute_hidden pnchisq_raw(double, double, double, double, double, int, Rboolean, Rboolean); diff --git a/src/pbeta.c b/src/pbeta.c index b62f265..98a98f5 100644 --- a/src/pbeta.c +++ b/src/pbeta.c @@ -1,6 +1,6 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 2006 The R Core Team + * Copyright (C) 2022 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -52,6 +52,9 @@ double pbeta_raw(double x, double a, double b, int lower_tail, int log_p) // else, remaining case: a = b = Inf : point mass 1 at 1/2 if (x < 0.5) return R_DT_0; else return R_DT_1; } + if (x >= 1) // may happen when called from qbeta() + return R_DT_1; + // Now: 0 < a < Inf; 0 < b < Inf double x1 = 0.5 - x + 0.5, w, wc; @@ -61,8 +64,8 @@ double pbeta_raw(double x, double a, double b, int lower_tail, int log_p) //==== // ierr in {10,14} <==> bgrat() error code ierr-10 in 1:4; for 1 and 4, warned *there* if(ierr && ierr != 11 && ierr != 14) - MATHLIB_WARNING4(_("pbeta_raw(%g, a=%g, b=%g, ..) -> bratio() gave error code %d"), - x, a,b, ierr); + MATHLIB_WARNING6(_("pbeta_raw(%g, a=%g, b=%g, lower=%d, log=%d) -> bratio() gave error code %d"), + x, a,b, lower_tail, log_p, ierr); return lower_tail ? w : wc; } /* pbeta_raw() */ @@ -77,8 +80,6 @@ double pbeta(double x, double a, double b, int lower_tail, int log_p) if (x <= 0) return R_DT_0; - if (x >= 1) - return R_DT_1; return pbeta_raw(x, a, b, lower_tail, log_p); } diff --git a/src/pgamma.c b/src/pgamma.c index 01e88d1..4dc8b58 100644 --- a/src/pgamma.c +++ b/src/pgamma.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2006-2019 The R Core Team * Copyright (C) 2005-6 Morten Welinder * Copyright (C) 2005-10 The R Foundation - * Copyright (C) 2006-2015 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -67,7 +67,7 @@ static const double M_cutoff = M_LN2 * DBL_MAX_EXP / DBL_EPSILON;/*=3.196577e18* /* Continued fraction for calculation of * 1/i + x/(i+d) + x^2/(i+2*d) + x^3/(i+3*d) + ... = sum_{k=0}^Inf x^k/(i+k*d) * - * auxilary in log1pmx() and lgamma1p() + * auxiliary in log1pmx() and lgamma1p() */ static double logcf (double x, double i, double d, @@ -145,6 +145,9 @@ double log1pmx (double x) /* Compute log(gamma(a+1)) accurately also for small a (0 < a < 0.5). */ double lgamma1p (double a) { + if (fabs (a) >= 0.5) + return lgammafn (a + 1); + const double eulers_const = 0.5772156649015328606065120900824024; /* coeffs[i] holds (zeta(i+2)-1)/(i+2) , i = 0:(N-1), N = 40 : */ @@ -194,11 +197,6 @@ double lgamma1p (double a) const double c = 0.2273736845824652515226821577978691e-12;/* zeta(N+2)-1 */ const double tol_logcf = 1e-14; - double lgam; - int i; - - if (fabs (a) >= 0.5) - return lgammafn (a + 1); /* Abramowitz & Stegun 6.1.33 : for |x| < 2, * <==> log(gamma(1+x)) = -(log(1+x) - x) - gamma*x + x^2 * \sum_{n=0}^\infty c_n (-x)^n @@ -207,8 +205,8 @@ double lgamma1p (double a) * Here, another convergence acceleration trick is used to compute * lgam(x) := sum_{n=0..Inf} c_n (-x)^n */ - lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf); - for (i = N - 1; i >= 0; i--) + double lgam = c * logcf(-a / 2, N + 2, 1, tol_logcf); + for (int i = N - 1; i >= 0; i--) lgam = coeffs[i] - a * lgam; return (a * lgam - eulers_const) * a - log1pmx (a); diff --git a/src/plogis.c b/src/plogis.c index 5d00cd3..df679fe 100644 --- a/src/plogis.c +++ b/src/plogis.c @@ -1,7 +1,7 @@ /* * R : A Computer Language for Statistical Data Analysis + * Copyright (C) 2000--2020 The R Core Team * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka - * Copyright (C) 2000 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -25,13 +25,16 @@ curve(log1p(exp(x)) - x, 33.1, 33.5, n=2^10) curve(x+exp(-x) - log1p(exp(x)), 15, 25, n=2^11) */ -double Rf_log1pexp(double x) { +double log1pexp(double x) { if(x <= 18.) return log1p(exp(x)); if(x > 33.3) return x; // else: 18.0 < x <= 33.3 : return x + exp(-x); } +// API. For now, continue using macro R_Log1_Exp() in our own code. +double log1mexp(double x) { return R_Log1_Exp(-x); } + double plogis(double x, double location, double scale, int lower_tail, int log_p) { @@ -47,7 +50,7 @@ double plogis(double x, double location, double scale, if(log_p) { // log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x)) - return -Rf_log1pexp(lower_tail ? -x : x); + return -log1pexp(lower_tail ? -x : x); } else { return 1 / (1 + exp(lower_tail ? -x : x)); } diff --git a/src/pnbeta.c b/src/pnbeta.c index 24a9d36..94c1681 100644 --- a/src/pnbeta.c +++ b/src/pnbeta.c @@ -30,7 +30,7 @@ pnbeta_raw(double x, double o_x, double a, double b, double ncp) const int itrmax = 10000; /* 100 is not enough for pf(ncp=200) see PR#11277 */ - double a0, lbeta, c, errbd, x0, temp, tmp_c; + double a0, lBeta, c, errbd, x0, temp, tmp_c; int ierr; LDOUBLE ans, ax, gx, q, sumq; @@ -46,15 +46,15 @@ pnbeta_raw(double x, double o_x, double a, double b, double ncp) x0 = floor(fmax2(c - 7. * sqrt(c), 0.)); a0 = a + x0; - lbeta = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b); + lBeta = lbeta(a0, b); // = lgammafn(a0) + lgammafn(b) - lgammafn(a0 + b); /* temp = pbeta_raw(x, a0, b, TRUE, FALSE), but using (x, o_x): */ bratio(a0, b, x, o_x, &temp, &tmp_c, &ierr, FALSE); gx = exp(a0 * log(x) + b * (x < .5 ? log1p(-x) : log(o_x)) - - lbeta - log(a0)); - if (a0 > a) + - lBeta - log(a0)); + if (a0 > a) // x0 >= 1 (and *not* x0 << a) q = exp(-c + x0 * log(c) - lgammafn(x0 + 1.)); - else + else // a0 = a <== x0 << a q = exp(-c); sumq = 1. - q; diff --git a/src/pnchisq.c b/src/pnchisq.c index ae85eb4..672e1c4 100644 --- a/src/pnchisq.c +++ b/src/pnchisq.c @@ -23,7 +23,7 @@ /*----------- DEBUGGING ------------- * * make CFLAGS='-DDEBUG_pnch ....' -(cd `R-devel RHOME`/src/nmath; gcc -I. -I../../src/include -I../../../R/src/include -I/usr/local/include -DHAVE_CONFIG_H -fopenmp -g -O0 -pedantic -Wall --std=gnu99 -DDEBUG_pnch -DDEBUG_q -Wcast-align -Wclobbered -c ../../../R/src/nmath/pnchisq.c -o pnchisq.o ) +(cd `R-devel RHOME`/src/nmath; gcc -I. -I../../src/include -I../../../R/src/include -I/usr/local/include -DHAVE_CONFIG_H -fopenmp -g -O0 -pedantic -Wall -DDEBUG_pnch -DDEBUG_q -c ../../../R/src/nmath/pnchisq.c -o pnchisq.o ) * -- Feb.6, 2000 (R pre0.99); M.Maechler: still have * bad precision & non-convergence in some cases (x ~= f, both LARGE) @@ -76,7 +76,7 @@ double pnchisq(double x, double df, double ncp, int lower_tail, int log_p) * <---> "in principle" this check should happen there, not here */ if (!log_p || ans < -1e-8) return ans; - else { // log_p && ans >= -1e-8 + else { // log_p (==> ans <= 0) && -1e-8 <= ans <= 0 // prob. = exp(ans) is near one: we can do better using the other tail #ifdef DEBUG_pnch REprintf(" pnchisq_raw(*, log_p): ans=%g => 2nd call, other tail\n", ans); @@ -125,13 +125,13 @@ pnchisq_raw(double x, double f, double theta /* = ncp */, log(x) < M_LN2 + 2/f*(lgamma(f/2. + 1) + _dbl_min_exp)) { // all pchisq(x, f+2*i, lower_tail, FALSE), i=0,...,110 would underflow to 0. // ==> work in log scale - double lambda = 0.5 * theta; - double sum, sum2, pr = -lambda; + double lambda = 0.5 * theta; // < 40 + double sum, sum2, pr = -lambda, log_lam = log(lambda); sum = sum2 = ML_NEGINF; /* we need to renormalize here: the result could be very close to 1 */ - for(i = 0; i < 110; pr += log(lambda) - log(++i)) { + for(i = 0; i < 110; pr += log_lam - log(++i)) { sum2 = logspace_add(sum2, pr); - sum = logspace_add(sum, pr + pchisq(x, f+2*i, lower_tail, TRUE)); + sum = logspace_add(sum , pr + pchisq(x, f+2*i, lower_tail, TRUE)); if (sum2 >= -1e-15) /*<=> EXP(sum2) >= 1-1e-15 */ break; } ans = sum - sum2; @@ -171,9 +171,7 @@ pnchisq_raw(double x, double f, double theta /* = ncp */, lam = .5 * theta; // = lambda = ncp/2 lamSml = (-lam < _dbl_min_exp); if(lamSml) { - /* MATHLIB_ERROR( - "non centrality parameter (= %g) too large for current algorithm", -p theta) */ + // originally error: "non centrality parameter too large for current algorithm" u = 0; lu = -lam;/* == ln(u) */ l_lam = log(lam); @@ -231,7 +229,9 @@ p theta) */ for (n = 1, f_2n = f + 2., f_x_2n += 2.; n <= itrmax ; n++, f_2n += 2, f_x_2n += 2) { #ifdef DEBUG_pnch_n - REprintf("\n _OL_: n=%d",n); + if(n % 1000 == 0) + REprintf("\n _OL_: n=%d, f_x_2n = %g", n); + else REprintf(n % 100 == 0 ? ".\n" : "."); #endif #ifndef MATHLIB_STANDALONE if(n % 1000 == 0) R_CheckUserInterrupt(); @@ -243,7 +243,9 @@ p theta) */ bound = (double) (t * x / f_x_2n); #ifdef DEBUG_pnch_n - REprintf("\n L10: n=%d; term= %g; bound= %g",n,term,bound); + if(n % 1000 == 0) + REprintf("\n L10: n=%d; term, ans = %g, %g; bound= %g", + n, term, ans, bound); #endif is_r = FALSE; /* convergence only if BOTH absolute and relative error < 'bnd' */ @@ -251,9 +253,9 @@ p theta) */ (is_r = (term <= reltol * ans)))) { #ifdef DEBUG_pnch - REprintf("BREAK out of for(n = 1 ..): n=%d; bound= %g %s, rel.err= %g %s\n", + REprintf("BREAK from for(n=1 ..): n=%d; bound= %g %s; term=%g, rel.err= %g %s\n", n, - bound, (is_b ? "<= errmax" : ""), + bound, (is_b ? "<= errmax" : ""), term, term/ans, (is_r ? "<= reltol" : "")); #endif break; /* out completely */ diff --git a/src/pnorm.c b/src/pnorm.c index 3578645..9815a01 100644 --- a/src/pnorm.c +++ b/src/pnorm.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2013 The R Core Team + * Copyright (C) 2000-2020 The R Core Team * Copyright (C) 2003 The R Foundation + * Copyright (C) 1998 Ross Ihaka * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -86,7 +86,6 @@ double pnorm5(double x, double mu, double sigma, int lower_tail, int log_p) return(lower_tail ? p : cp); } -#define SIXTEN 16 /* Cutoff allowing exact "*" and "/" */ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p) { @@ -193,18 +192,20 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p) } temp = (xnum + c[7]) / (xden + d[7]); -#define do_del(X) \ - xsq = trunc(X * SIXTEN) / SIXTEN; \ - del = (X - xsq) * (X + xsq); \ - if(log_p) { \ - *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + log(temp); \ - if((lower && x > 0.) || (upper && x <= 0.)) \ - *ccum = log1p(-exp(-xsq * xsq * 0.5) * \ - exp(-del * 0.5) * temp); \ - } \ - else { \ - *cum = exp(-xsq * xsq * 0.5) * exp(-del * 0.5) * temp; \ - *ccum = 1.0 - *cum; \ +#define d_2(_x_) ldexp(_x_, -1) // == (_x_ / 2 ) "perfectly" + +#define do_del(X) \ + xsq = ldexp(trunc(ldexp(X, 4)), -4); \ + del = (X - xsq) * (X + xsq); \ + if(log_p) { \ + *cum = (-xsq * d_2(xsq)) -d_2(del) + log(temp); \ + if((lower && x > 0.) || (upper && x <= 0.)) \ + *ccum = log1p(-exp(-xsq * d_2(xsq)) * \ + exp(-d_2(del)) * temp); \ + } \ + else { \ + *cum = exp(-xsq * d_2(xsq)) * exp(-d_2(del)) * temp;\ + *ccum = 1.0 - *cum; \ } #define swap_tail \ @@ -225,8 +226,8 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p) * Note that we do want symmetry(0), lower/upper -> hence use y */ else if((log_p && y < 1e170) /* avoid underflow below */ - /* ^^^^^ MM FIXME: can speedup for log_p and much larger |x| ! - * Then, make use of Abramowitz & Stegun, 26.2.13, something like + /* ^^^^^ MM FIXME: could speedup for log_p and |x| >> 5.657 ! + * Then, make use of Abramowitz & Stegun, 26.2.13, p.932, something like xsq = x*x; @@ -239,8 +240,8 @@ void pnorm_both(double x, double *cum, double *ccum, int i_tail, int log_p) swap_tail; - [Yes, but xsq might be infinite.] - + Yes, but xsq might be infinite.; + well, actually x = -1.34..e154 = -sqrt(DBL_MAX) already overflows x^2 */ || (lower && -37.5193 < x && x < 8.2924) || (upper && -8.2924 < x && x < 37.5193) diff --git a/src/polygamma.c b/src/polygamma.c index f082b92..41dd05b 100644 --- a/src/polygamma.c +++ b/src/polygamma.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2015 The R Core Team + * Copyright (C) 2000-2022 The R Core Team * Copyright (C) 2004-2009 The R Foundation + * Copyright (C) 1998 Ross Ihaka * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -145,6 +145,32 @@ #define n_max (100) +// Compute d_n(x) = (d/dx)^n cot(x) ; cot(x) := cos(x) / sin(x) +double d_n_cot(double x, int n) +{ + if (n == 0) + return cos(x)/sin(x); + else if (n == 1) // -1/sin^2 + return -1/R_pow_di(sin(x), 2); + else if (n == 2) // 2 cos / sin^3 + return 2*cos(x)/R_pow_di(sin(x), 3); + else if (n == 3) { // (-4 cos^2 - 2) / sin^4 ; num.= -2(2 cos^2 +1) = -2(3 - 2 sin^2) + // tt = -2*(2*R_pow_di(cos(x), 2) + 1.)/R_pow_di(sin(x), 4); + double sin2 = R_pow_di(sin(x), 2); // = sin^2 + return -2*(3 - 2*sin2)/R_pow_di(sin2, 2); + } + else if (n == 4) { // 8 cos (cos^2 + 2) / sin^5 + double co = cos(x); + return 8*co*(R_pow_di(co, 2) + 2) / R_pow_di(sin(x), 5); + } + else if (n == 5) { // (-16 cos^4 - 88 cos^2 - 16)/ sin^6 + double co2 = R_pow_di(cos(x), 2); // cos^2 + return -8*(2*R_pow_di(co2, 2) + 11*co2 + 2) / R_pow_di(sin(x), 6); + } + else + return ML_NAN; +} + /* From R, currently only used for kode = 1, m = 1 : */ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) { @@ -173,21 +199,24 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) -1.92965793419400681e+16 }; - int i, j, k, mm, mx, nn, np, nx, fn; - double arg, den, elim, eps, fln, fx, rln, rxsq, - r1m4, r1m5, s, slope, t, ta, tk, tol, tols, tss, tst, - tt, t1, t2, wdtol, xdmln, xdmy, xinc, xln = 0.0 /* -Wall */, - xm, xmin, xq, yint; + int i, j, k, nn, np, nx, fn; + double arg, den, eps, fx, rxsq, + s, t, ta, tk, tol, tols, tss, tst, + tt, t1, t2, xdmln, xdmy, xinc, + xm, xmin, xq; double trm[23], trmr[n_max + 1]; *ierr = 0; + *nz = 0; if (n < 0 || kode < 1 || kode > 2 || m < 1) { *ierr = 1; return; } if (x <= 0.) { - /* use Abramowitz & Stegun 6.4.7 "Reflection Formula" - * psi(k, x) = (-1)^k psi(k, 1-x) - pi^{n+1} (d/dx)^n cot(x) + /* use Abramowitz & Stegun 6.4.7 "Reflection Formula", p.260 + * psi(n, 1-x) + (-1)^(n+1) psi(n, x) = (-1)^n pi (d/dx)^n cot(pi*x) + * psi(n, x) = (-1)^n psi(n, 1-x) - pi^{n+1} d_n(pi*x), + where d_n(x) := (d/dx)^n cot(x) */ if (x == round(x)) { /* non-positive integer : +Inf or NaN depends on n */ @@ -201,44 +230,31 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) * for j = 0:(m-1) , k = n + j */ - /* Cheat for now: only work for m = 1, n in {0,1,2,3} : */ - if(m > 1 || n > 3) {/* doesn't happen for digamma() .. pentagamma() */ - /* not yet implemented */ + /* For now: only work for n in {0,1,..,5} : */ + if(n > 5) { + /* not yet implemented for x < 0 and n >= 6 */ *ierr = 4; return; } + // tt := d_n(pi * x) x *= M_PI; /* pi * x */ - if (n == 0) - tt = cos(x)/sin(x); - else if (n == 1) - tt = -1/R_pow_di(sin(x), 2); - else if (n == 2) - tt = 2*cos(x)/R_pow_di(sin(x), 3); - else if (n == 3) - tt = -2*(2*R_pow_di(cos(x), 2) + 1.)/R_pow_di(sin(x), 4); - else /* can not happen! */ - tt = ML_NAN; - /* end cheat */ - - s = (n % 2) ? -1. : 1.;/* s = (-1)^n */ - /* t := pi^(n+1) * d_n(x) / gamma(n+1) , where - * d_n(x) := (d/dx)^n cot(x)*/ + + // t := pi^(n+1) * d_n(x) / gamma(n+1) t1 = t2 = s = 1.; for(k=0, j=k-n; j < m; k++, j++, s = -s) { /* k == n+j , s = (-1)^k */ t1 *= M_PI;/* t1 == pi^(k+1) */ if(k >= 2) t2 *= k;/* t2 == k! == gamma(k+1) */ - if(j >= 0) /* by cheat above, tt === d_k(x) */ - ans[j] = s*(ans[j] + t1/t2 * tt); + if(j >= 0) /* now using d_k(x) */ + ans[j] = s*(ans[j] + t1/t2 * d_n_cot(x, k)); } - if (n == 0 && kode == 2) /* unused from R, but "wrong": xln === 0 :*/ - ans[0] += xln; + /* if (n == 0 && kode == 2) -- nonsense for x < 0 ! + * ans[0] += log(x); */ return; } /* x <= 0 */ /* else : x > 0 */ - *nz = 0; - xln = log(x); + double xln = log(x); if(kode == 1 && m == 1) {/* the R case --- for very large x: */ double lrg = 1/(2. * DBL_EPSILON); if(n == 0 && x * xln > lrg) { @@ -250,24 +266,29 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) return; } } - mm = m; nx = imin2(-Rf_i1mach(15), Rf_i1mach(16));/* = 1021 */ - r1m5 = Rf_d1mach(5); - r1m4 = Rf_d1mach(4) * 0.5; - wdtol = fmax2(r1m4, 0.5e-18); /* 1.11e-16 */ + + const double + r1m5 = Rf_d1mach(5), // = M_LOG10_2 = log10(2) = 0.30103.. + r1m4 = Rf_d1mach(4) * 0.5, // = DBL_EPSILON * 0.5 = 2^-53 = 1.110223e-16 + wdtol = fmax2(r1m4, 0.5e-18); /* = 2^-53 = 1.11e-16 */ /* elim = approximate exponential over and underflow limit */ - elim = 2.302 * (nx * r1m5 - 3.0);/* = 700.6174... */ + double elim = 2.302 * (nx * r1m5 - 3.0);/* = 700.6174... */ + const double + rln = fmin2(r1m5 * Rf_i1mach(14), 18.06);// = 0.30103 * 53 = 15.95.. ~= #{decimals} + double fln = fmax2(rln, 3.0) - 3.0; // = 12.95.. + const double yint = 3.50 + 0.40 * fln, // = 8.6818.. + slope = 0.21 + fln * (0.0006038 * fln + 0.008677); // = 0.4237.. + int mm = m; for(;;) { nn = n + mm - 1; fn = nn; t = (fn + 1) * xln; /* overflow and underflow test for small and large x */ - if (fabs(t) > elim) { if (t <= 0.0) { - *nz = 0; *ierr = 2; return; } @@ -284,24 +305,15 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) return; } - /* compute xmin and the number of terms of the series, fln+1 */ - - rln = r1m5 * Rf_i1mach(14); - rln = fmin2(rln, 18.06); - fln = fmax2(rln, 3.0) - 3.0; - yint = 3.50 + 0.40 * fln; - slope = 0.21 + fln * (0.0006038 * fln + 0.008677); + /* compute xmin and the number of terms of the series, fln+1 */ xm = yint + slope * fn; - mx = (int)xm + 1; + int mx = (int)xm + 1; xmin = mx; if (n != 0) { xm = -2.302 * rln - fmin2(0.0, xln); - arg = xm / n; - arg = fmin2(0.0, arg); + arg = fmin2(0., xm / n); eps = exp(arg); - xm = 1.0 - eps; - if (fabs(arg) < 1.0e-3) - xm = -arg; + xm = (fabs(arg) < 1.0e-3) ? -arg : 1. - eps; fln = x * xm / eps; xm = xmin - x; if (xm > 7.0 && fln < 15.0) @@ -326,7 +338,7 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) if (tk <= elim) /* for all but large x */ goto L10; } - nz++; /* underflow */ + nz++; /* nz := #{underflows} */ mm--; ans[mm] = 0.; if (mm == 0) @@ -397,7 +409,6 @@ void dpsifn(double x, int n, int kode, int m, double *ans, int *nz, int *ierr) nx = (int)xinc; np = nn + 1; if (nx > n_max) { - *nz = 0; *ierr = 3; return; } diff --git a/src/qDiscrete_search.h b/src/qDiscrete_search.h new file mode 100644 index 0000000..0bdbcb4 --- /dev/null +++ b/src/qDiscrete_search.h @@ -0,0 +1,219 @@ +/* + * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2021 The R Core Team + * Copyright (C) 2005-2021 The R Foundation + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, a copy is available at + * https://www.R-project.org/Licenses/ + */ + +/* This is #included from ./qnbinom.c , ......... +*/ + +#define PST_0(a, b) a ## b +#define PASTE(a, b) PST_0(a, b) + +#define CHR_0(x) #x +#define AS_CHAR(x) CHR_0(x) + +#define _pDIST_ PASTE(p, _thisDIST_) +#define _qDIST_ PASTE(q, _thisDIST_) +/** + For a discrete distribution on the integers, + For P(x) := (x, ), find p-quantile y(p) :<==> P(y) < p <= P(y) + + @param y current guess + @param *z := (y, ..) + @param p target probability + @param parameters of the respective distribution function + @param incr increment, integer valued >= 1. + @param lower_tail "logical" int; if 1 (true), have lower tail probabilities; otherwise upper tail + @param log "logical" int; if 1 (true) the probabilities are log() + + @return root 'y' and z[0] = (y, ..) + */ +#define DO_SEARCH_FUN(...) \ + do_search(double y, double *z, double p, __VA_ARGS__, \ + double incr, int lower_tail, int log_p) + +// this is used in the caller file qnbinom.c etc : +#define DO_SEARCH_(Y_, incr_, ...) do_search(Y_, &z, p, __VA_ARGS__, incr_, lower_tail, log_p) + +#define P_DIST(Y_, ...) _pDIST_(Y_, __VA_ARGS__, lower_tail, log_p) + +#ifdef MATHLIB_STANDALONE +# define MAYBE_R_CheckUserInterrupt() +#else +# define MAYBE_R_CheckUserInterrupt() R_CheckUserInterrupt() +#endif + +static double DO_SEARCH_FUN(_dist_PARS_DECL_) +{ + Rboolean left = (lower_tail ? (*z >= p) : (*z < p)); + R_DBG_printf(" do_search(y=%g, z=%15.10g %s p = %15.10g) --> search to %s", + y, *z, (lower_tail ? ">=" : "<"), p, (left ? "left" : "right")); + if(incr > 1) + R_DBG_printf(", incr = %.0f\n", incr); + else R_DBG_printf("\n"); + + if(left) { // (lower_tail, *z >= p) or (upper tail, *z < p): search to the __left__ + for(int iter = 0; ; iter++) { + double newz = -1.; // -Wall +#ifndef MATHLIB_STANDALONE + if(iter % 10000 == 0) R_CheckUserInterrupt();// have seen inf.loops +#endif + if(y > 0) + newz = P_DIST(y - incr, _dist_PARS_); + else if(y < 0) + y = 0; + // note that newz may be NaN because of remaining border line bugs in _pDIST_() {bug from pbeta()} + if(y == 0 || ISNAN(newz) || (lower_tail ? (newz < p) : (newz >= p))) { + R_DBG_printf(" new y=%.15g, " AS_CHAR(_pDIST_) "(y-incr,*) %s;" + " ==> search() returns previous z=%g after %d iter.\n", y, + ISNAN(newz) ? "is NaN" : (lower_tail ? "< p" : ">= p"), *z, iter); + return y; // and previous *z + } + y = fmax2(0, y - incr); + *z = newz; + } + } + else { // (lower_tail, *z < p) or (upper tail, *z >= p): search to the __right__ + for(int iter = 0; ; iter++) { +#ifndef MATHLIB_STANDALONE + if(iter % 10000 == 0) R_CheckUserInterrupt(); +#endif + y += incr; +#ifdef _dist_MAX_y + if(y < _dist_MAX_y) + *z = P_DIST(y, _dist_PARS_); + else if(y > _dist_MAX_y) + y = _dist_MAX_y; +#else + *z = P_DIST(y, _dist_PARS_); +#endif + + if( +#ifdef _dist_MAX_y + y == _dist_MAX_y || +#endif + ISNAN(*z) || (lower_tail ? (*z >= p) : (*z < p))) + { + R_DBG_printf(" new y=%.15g, z=%g = " AS_CHAR(_pDIST_) "(y,*) %s;" + " ==> search() returns after %d iter.\n", y, *z, + ISNAN(*z) ? "is NaN" : (lower_tail ? ">= p" : "< p"), iter); + return y; + } + } + } +} // do_search() + + +/* +* Note : "same" code in qbinom.c, qnbinom.c __FIXME__ NOT YET for qpois() ?? +* FIXME: This is far from optimal [cancellation for p ~= 1, etc]: +*/ +#define q_DISCRETE_01_CHECKS() do { \ + double p_n; /* p in the "normal" (lower_tail=TRUE, log.p=FALSE) scale */ \ + if(!lower_tail || log_p) { \ + p_n = R_DT_qIv(p); /* need check again (cancellation!): */ \ + R_DBG_printf(" upper tail or log_p: (p_n, 1-p_n) = (%.15g, %.15g)\n", p_n, 1-p_n); \ + /* _OLD_ (R <= 4.0.x): */ \ + /* if (p_n == 0) return 0; */ \ + /* if (p_n == 1) return ML_POSINF; */ \ + /* end{_OLD_} */ \ + if (p_n == 0) R_DBG_printf("** p_n == 0: _NO LONGER_ returning 0"); \ + if (p_n == 1) R_DBG_printf("** p_n == 1: _NO LONGER_ returning +Inf"); \ + } else \ + p_n = p; \ + /* temporary hack : __New: On 2020-08-26, only give message but do *NOT* early return: */ \ + /* seen (only @sfs, not 'nb-mm' (???): infinite loop for chkQnbinom(1e-19, k.max=1e4) */ \ + if (p_n + 1.01*DBL_EPSILON >= 1.) { \ + R_DBG_printf("p_n + 1.01*eps >= 1 ; (1-p_n = %g): for now *NO LONGER* returning Inf\n", \ + 1-p_n); \ + /* return ML_POSINF; */ \ + } \ +} while(0) + + +#ifdef _dist_MAX_y +# define q_DISCR_CHECK_BOUNDARY(Y_) do { \ + if(Y_ > _dist_MAX_y) /* way off */ Y_ = _dist_MAX_y; \ + else if(Y_ < 0) Y_ = 0.; \ +} while(0) +#else +# define q_DISCR_CHECK_BOUNDARY(Y_) if(Y_ < 0) Y_ = 0. + /* e.g., for qnbinom(0.5, mu = 3, size = 1e-10) */ +#endif + +#define q_DISCRETE_BODY() do { \ + /* y := approx.value (Cornish-Fisher expansion) : */ \ + double \ + z = qnorm(p, 0., 1., lower_tail, log_p), \ + y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6)); \ + R_DBG_printf(" Cornish-Fisher: initial z = qnorm(p, l.t, log)= %g, y = %g;\n", z,y); \ + \ + q_DISCR_CHECK_BOUNDARY(y); \ + \ + z = P_DIST(y, _dist_PARS_); \ + \ + /* Algorithmic "tuning parameters", used to be hardwired; changed for speed &| precision */ \ + const double \ + _pf_n_ = 8, /* was hardwired to 64 */ \ + _pf_L_ = 2, /* was hardwired to 64 */ \ + _yLarge_ = 4096, /* was hardwired to 1e5 */ \ + _incF_ = (1./64),/* was hardwired to 0.001 (= 1./1000 ) */ \ + _iShrink_ = 8, /* was hardwired to 100 */ \ + _relTol_ = 1e-15,/* was hardwired to 1e-15 */ \ + _xf_ = 4; /* extra factor, *must* be >= 1 (new anyway) */ \ + \ + R_DBG_printf(" algo. tuning: fuzz factors _pf_{n,L}: {%.0f, %.0f}; yLarge = %g\n" \ + " large case: _incF_=%g, _iShrink_=%g; _relTol_=%g, _xf_=%g\n", \ + _pf_n_, _pf_L_, _yLarge_, _incF_, _iShrink_, _relTol_, _xf_); \ + \ + /* fuzz to ensure left continuity: do not loose too much (=> error in upper tail) */ \ + if(log_p) { /* <==> p \in [-Inf, 0] different adjustment: "other sign" */ \ + double e = _pf_L_ * DBL_EPSILON; \ + if(lower_tail && p > - DBL_MAX) /* prevent underflow to -Inf */ \ + p *= 1 + e; \ + else /* if(p < - DBL_MIN) // not too close to -0 */ \ + p *= 1 - e; \ + \ + } else { /* not log scale */ \ + double e = _pf_n_ * DBL_EPSILON; \ + if(lower_tail) \ + p *= 1 - e; \ + else if(1 - p > _xf_*e) /* otherwise get p > 1 */ \ + p *= 1 + e; \ + } \ + R_DBG_printf(" new z := " AS_CHAR(_pDIST_) "(y, *) = %.11g; left-cont. fuzz => p = %.11g\n", z, p); \ + \ + /* If the C-F value y is not too large a simple search is OK */ \ + if(y < _yLarge_) return DO_SEARCH_(y, 1, _dist_PARS_); \ + /* Otherwise be a bit cleverer in the search: use larger increments, notably initially: */ \ + { /* y >= _yLarge_ */ \ + double oldincr, incr = floor(y * _incF_); \ + int qIt = 0; \ + \ + R_DBG_printf(" large y: --> use larger increments than 1: incr=%.0f\n", incr); \ + do { \ + oldincr = incr; \ + y = DO_SEARCH_(y, incr, _dist_PARS_); /* also updating *z */ \ + if(++qIt % 10000 == 0) MAYBE_R_CheckUserInterrupt(); \ + incr = fmax2(1, floor(incr / _iShrink_)); \ + } while(oldincr > 1 && incr > y * _relTol_); \ + R_DBG_printf(" \\--> oldincr=%.0f, after %d \"outer\" search() iterations\n", \ + oldincr, qIt); \ + return y; \ + } \ +} while(0) diff --git a/src/qbeta.c b/src/qbeta.c index 2352c1b..5ce8798 100644 --- a/src/qbeta.c +++ b/src/qbeta.c @@ -1,6 +1,6 @@ /* * R : A Computer Language for Statistical Data Analysis - * Copyright (C) 1998--2018 The R Core Team + * Copyright (C) 1998--2022 The R Core Team * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka * based on code (C) 1979 and later Royal Statistical Society * @@ -30,6 +30,12 @@ #include "nmath.h" #include "dpq.h" +/**----------- DEBUGGING ------------- + * + * make CFLAGS='-DDEBUG_qbeta ...' + *MM (w/ Debug, w/o Optimization): if not manually adding -dDEBUG_qbeta , use + (cd `R-devel-qbeta-dbg RHOME`/src/nmath ; make qbeta.o; cd ../..; make R) +*/ #ifdef DEBUG_qbeta # define R_ifDEBUG_printf(...) REprintf(__VA_ARGS__) #else @@ -63,6 +69,7 @@ double qbeta(double alpha, double p, double q, int lower_tail, int log_p) double qbet[2];// = { qbeta(), 1 - qbeta() } qbeta_raw(alpha, p, q, lower_tail, log_p, + // swap_01 , log_q_cut , n_N MLOGICAL_NA, USE_LOG_X_CUTOFF, n_NEWTON_FREE, qbet); return qbet[0]; } @@ -75,7 +82,7 @@ static const double // Too extreme: inaccuracy in pbeta(); e.g for qbeta(0.95, 1e-9, 20): // -> in pbeta() --> bgrat(..... b*z == 0 underflow, hence inaccurate pbeta() /* DBL_very_MIN = 0x0.0000001p-1022, // = 2^-1050 = 2^(-1022 - 28) */ - /* DBL_log_v_MIN = -1050. * M_LN2, // = log(DBL_very_MIN) */ + /* DBL_log_v_MIN = -1050. * M_LN2, // = log(DBL_very_MIN) ~= -727.8045 */ // the most extreme -- not ok, as pbeta() then behaves strangely, // e.g., for qbeta(0.95, 1e-8, 20): /* DBL_very_MIN = 0x0.0000000000001p-1022, // = 2^-1074 = 2^(-1022 -52) */ @@ -141,7 +148,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, return_q_1; } - // check alpha {*before* transformation which may all accuracy}: + // check alpha {*before* transformation which may lose all accuracy}: if((log_p && alpha > 0) || (!log_p && (alpha < 0 || alpha > 1))) { // alpha is outside R_ifDEBUG_printf("qbeta(alpha=%g, %g, %g, .., log_p=%d): %s%s\n", @@ -183,9 +190,19 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, logbeta = lbeta(p, q); swap_tail = (swap_choose) ? (p_ > 0.5) : swap_01; + + R_ifDEBUG_printf( + "qbeta(%g, %g, %g, lower_t=%d, log_p=%d, swap_01=%d, log_q_cut=%g, n_N=%d):%s\n" + " swap_tail=%s :", + alpha, p,q, lower_tail, log_p, swap_01, log_q_cut, n_N, + (log_p && (p_ == 0. || p_ == 1.)) ? (p_==0.?" p_=0":" p_=1") : "", + (swap_tail ? "TRUE": "F")); + + int n_maybe_swaps = 0; +maybe_swap: // change tail; default (swap_01 = NA): afterwards 0 < a <= 1/2 - if(swap_tail) { /* change tail, swap p <-> q :*/ - a = R_DT_CIv(alpha); // = 1 - p_ < 1/2 + if(swap_tail) { /* change tail, swap copies of {p,q}: p <-> q :*/ + a = R_DT_CIv(alpha); // = 1 - p_ , is < 1/2 if(swap_choose) /* la := log(a), but without numerical cancellation: */ la = R_DT_Clog(alpha); pp = q; qq = p; @@ -195,6 +212,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, la = R_DT_log(alpha); pp = p; qq = q; } + n_maybe_swaps++; /* calculate the initial approximation */ @@ -209,37 +227,37 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, */ double acu = fmax2(acu_min, pow(10., -13. - 2.5/(pp * pp) - 0.5/(a * a))); // try to catch "extreme left tail" early - double tx, u0 = (la + log(pp) + logbeta) / pp; // = log(x_0) + double tx, + u0 = (la + log(pp) + logbeta) / pp, // = log(x_0) + rp = pp*(1.-qq)/(pp+1.); static const double log_eps_c = M_LN2 * (1. - DBL_MANT_DIG);// = log(DBL_EPSILON) = -36.04.. - r = pp*(1.-qq)/(pp+1.); t = 0.2; // FIXME: Factor 0.2 is a bit arbitrary; '1' is clearly much too much. - + Rboolean u0_maybe = (M_LN2 * DBL_MIN_EXP < u0 && u0 < -0.01); + /* 1. cannot allow exp(u0) = 0 ==> exp(u1) = exp(u0) = 0 + * 2. must: u0 < 0, but too close to 0 <==> x = exp(u0) = 0.99.. */ R_ifDEBUG_printf( - "qbeta(%g, %g, %g, lower_t=%d, log_p=%d):%s\n" - " swap_tail=%d, la=%#8g, u0=%#8g (bnd: %g (%g)) ", - alpha, p,q, lower_tail, log_p, - (log_p && (p_ == 0. || p_ == 1.)) ? (p_==0.?" p_=0":" p_=1") : "", - swap_tail, la, u0, - (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2., - t*log_eps_c - log(fabs(r)) + " n_maybe_swaps=%d, la=%#8g, u0=%#8g -> u0_maybe=%s (bnd: %g (%g)); ", + n_maybe_swaps, la, u0, (u0_maybe ? "TRUE" : "F"), + u0_maybe ? (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2. : -0., + u0_maybe ? t*log_eps_c - log(fabs(rp)) : -0. ); - if(M_LN2 * DBL_MIN_EXP < u0 && // cannot allow exp(u0) = 0 ==> exp(u1) = exp(u0) = 0 - u0 < -0.01 && // (must: u0 < 0, but too close to 0 <==> x = exp(u0) = 0.99..) + double u_n = 1.; // to be log(xinbta) <==> xinbta = exp(u_n). 1 is impossible + if(u0_maybe && // qq <= 2 && // <--- "arbitrary" - // u0 < t*log_eps_c - log(fabs(r)) && + // u0 < t*log_eps_c - log(fabs(rp)) && u0 < (t*log_eps_c - log(fabs(pp*(1.-qq)*(2.-qq)/(2.*(pp+2.)))))/2.) { // TODO: maybe jump here from below, when initial u "fails" ? // L_tail_u: // MM's one-step correction (cheaper than 1 Newton!) - r = r*exp(u0);// = r*x0 - if(r > -1.) { - u = u0 - log1p(r)/pp; - R_ifDEBUG_printf("u1-u0=%9.3g --> choosing u = u1\n", u-u0); + rp = rp*exp(u0);// = rp*x0 + if(rp > -1.) { + u = u0 - log1p(rp)/pp; + R_ifDEBUG_printf("u1-u0=%9.3g --> choosing u = u1 = %g\n", u-u0, u); } else { u = u0; R_ifDEBUG_printf("cannot cheaply improve u0\n"); @@ -252,14 +270,14 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, // y := y_\alpha in AS 64 := Hastings(1955) approximation of qnorm(1 - a) : r = sqrt(-2 * la); y = r - (const1 + const2 * r) / (1. + (const3 + const4 * r) * r); - + R_ifDEBUG_printf("y_a(AS64)=%g\n ", y); if (pp > 1 && qq > 1) { // use Carter(1947), see AS 109, remark '5.' r = (y * y - 3.) / 6.; s = 1. / (pp + pp - 1.); t = 1. / (qq + qq - 1.); h = 2. / (s + t); w = y * sqrt(h + r) / h - (t - s) * (r + 5. / 6. - 2. / (3. * h)); - R_ifDEBUG_printf("p,q > 1 => w=%g", w); + R_ifDEBUG_printf("p,q > 1 => w=%g <= 300 ? ", w); if(w > 300) { // exp(w+w) is huge or overflows t = w+w + log(qq) - log(pp); // = argument of log1pexp(.) u = // log(xinbta) = - log1p(qq/pp * exp(w+w)) = -log(1 + exp(t)) @@ -274,33 +292,44 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, r = qq + qq; /* A slightly more stable version of t := \chi^2_{alpha} of AS 64 * t = 1. / (9. * qq); t = r * R_pow_di(1. - t + y * sqrt(t), 3); */ - t = 1. / (3. * sqrt(qq)); + t = 1. / (3. * sqrt(qq)); // = sqrt(t) of formula above t = r * R_pow_di(1. + t*(-t + y), 3);// = \chi^2_{alpha} of AS 64 - s = 4. * pp + r - 2.;// 4p + 2q - 2 = numerator of new t = (...) / chi^2 - R_ifDEBUG_printf("min(p,q) <= 1: t=%g", t); + s = 4. * pp + r - 2.;// 4p + 2q - 2 = numerator of new t' = s / t = s / chi^2 + R_ifDEBUG_printf("min(p,q) <= 1: t=%g, s=%g", t, s); if (t == 0 || (t < 0. && s >= t)) { // cannot use chisq approx // x0 = 1 - { (1-a)*q*B(p,q) } ^{1/q} {AS 65} // xinbta = 1. - exp((log(1-a)+ log(qq) + logbeta) / qq); - double l1ma;/* := log(1-a), directly from alpha (as 'la' above): - * FIXME: not worth it? log1p(-a) always the same ?? */ - if(swap_tail) - l1ma = R_DT_log(alpha); - else - l1ma = R_DT_Clog(alpha); - R_ifDEBUG_printf(" t <= 0 : log1p(-a)=%.15g, better l1ma=%.15g\n", log1p(-a), l1ma); + double l1ma = /* := log(1-a), directly from alpha (as 'la' above); + though only seen very small improvements */ + swap_tail ? R_DT_log(alpha) : R_DT_Clog(alpha); + + R_ifDEBUG_printf( + " t <= 0 => AS65: log1p(-a)=%7g, better l1ma=%.15g, relD=%g\n", + log1p(-a), l1ma, log1p(-a)/l1ma - 1); + double xx = (l1ma + log(qq) + logbeta) / qq; + R_ifDEBUG_printf(" xx = (l1ma + log(qq) + logbeta) / qq = %.10g; ", xx); if(xx <= 0.) { xinbta = -expm1(xx); u = R_Log1_Exp (xx);// = log(xinbta) = log(1 - exp(...A...)) + R_ifDEBUG_printf(" xx <= 0 useful ==> u, xinbta\n"); + } else { // xx > 0 ==> 1 - e^xx < 0 .. is nonsense - R_ifDEBUG_printf(" xx=%g > 0: xinbta:= 1-e^xx < 0\n", xx); - xinbta = 0; u = ML_NEGINF; /// FIXME can do better? + R_ifDEBUG_printf(" xx > 0: xinbta:= 1-e^xx < 0 'nonsense'; "); + // Try MM's one-step correction (or else u0) + double r_ = rp*exp(u0); + if(r_ > -1.) { + u = u0 - log1p(r_)/pp; R_ifDEBUG_printf("u:= u0 - log1p(r_)/pp = %g\n"); + } else { + u = u0; R_ifDEBUG_printf("u:= u0 = %g\n"); + } + xinbta = exp(u); } } else { t = s / t; - R_ifDEBUG_printf(" t > 0 or s < t < 0: new t = %g ( > 1 ?)\n", t); + R_ifDEBUG_printf(" t > 0 or s < t < 0: new t := s/t = %g ( > 1 ?)\n", t); if (t <= 1.) { // cannot use chisq, either - u = (la + log(pp) + logbeta) / pp; + u = u0; xinbta = exp(u); } else { // (1+x0)/(1-x0) = t, solved for x0 : xinbta = 1. - 2. / (t + 1.); @@ -310,25 +339,29 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, } // Problem: If initial u is completely wrong, we make a wrong decision here - if(swap_choose && + if(swap_choose && // vvvv/ why -exp(*)? u on log-x scale! Swapping can be very good, but needs smart if(..) (( swap_tail && u >= -exp( log_q_cut)) || // ==> "swap back" (!swap_tail && u >= -exp(4*log_q_cut) && pp / qq < 1000.) // ==> "swap now" )) { - // "revert swap" -- and use_log_x + if(swap_tail) + R_ifDEBUG_printf("\"swap back\" as u = %g >= -exp(log_q_cut);", u); + else + R_ifDEBUG_printf("\"swap now\" as u = %g >= -exp(4*log_q_cut) && pp/qq = %g < 1000;", + u, pp/qq); + // reverse swap (and typically use_log_x) swap_tail = !swap_tail; - R_ifDEBUG_printf(" u = %g (e^u = xinbta = %.16g) ==> ", u, xinbta); + if(swap_tail) { // "swap now" (much less easily) - a = R_DT_CIv(alpha); // needed ? + a = R_DT_CIv(alpha); // needed ? la = R_DT_Clog(alpha); pp = q; qq = p; } - else { // swap back : + else { // "swap back" : a = p_; la = R_DT_log(alpha); pp = p; qq = q; } - R_ifDEBUG_printf("\"%s\"; la = %g\n", - (swap_tail ? "swap now" : "swap back"), la); + R_ifDEBUG_printf(" ==> la = %g\n", la); // we could redo computations above, but this should be stable u = R_Log1_Exp(u); xinbta = exp(u); @@ -338,7 +371,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, 2) The correction step can go outside (u_n > 0 ==> e^u > 1 is illegal) e.g., for qbeta(0.2066, 0.143891, 0.05) */ - } else R_ifDEBUG_printf("\n"); + } if(!use_log_x) use_log_x = (u < log_q_cut);// <==> xinbta = e^u < exp(log_q_cut) @@ -351,7 +384,6 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, ((bad_init && !bad_u) ? ", ** bad_init **" : ""), (use_log_x ? ", on u = LOG(x) SCALE" : "")); - double u_n = 1.; // -Wall tx = xinbta; // keeping "original initial x" (for now) if(bad_u || u < log_q_cut) { @@ -362,7 +394,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, w = pbeta_raw(DBL_very_MIN, pp, qq, TRUE, log_p); if(w > (log_p ? la : a)) { R_ifDEBUG_printf( - " quantile is left of %g; \"convergence\"\n", DBL_very_MIN); + " quantile is left of %g: boundary \"convergence\"\n", DBL_very_MIN); if(log_p || fabs(w - a) < fabs(0 - a)) { // DBL_very_MIN is better than 0 tx = DBL_very_MIN; u_n = DBL_log_v_MIN;// = log(DBL_very_MIN) @@ -373,8 +405,9 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, use_log_x = log_p; add_N_step = FALSE; goto L_return; } else { - R_ifDEBUG_printf(" pbeta(%g, *) = %g <= %g (= %s) --> continuing\n", - DBL_log_v_MIN, w, (log_p ? la : a), (log_p ? "la" : "a")); + R_ifDEBUG_printf(" pbeta(%g, %g, %g, T, log) = %g <= %g (= %s) --> continuing\n", + DBL_very_MIN, pp,qq, w, + (log_p ? la : a), (log_p ? "la" : "a")); if(u < DBL_log_v_MIN) { u = DBL_log_v_MIN;// = log(DBL_very_MIN) xinbta = DBL_very_MIN; @@ -429,14 +462,23 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, */ w = (y == ML_NEGINF) // y = -Inf well possible: we are on log scale! ? 0. : (y - la) * exp(y - u + logbeta + r * u + t * R_Log1_Exp(u)); - if(!R_FINITE(w)) - break; + R_ifDEBUG_printf("N(i=%2d): u=%#20.16g, lnpb(e^u)=%#15.11g, w=%#12.7g", + i_pb, u, y, w); + if(!R_FINITE(w)) { + // what we should do is ==> go back, get better starting value + R_ifDEBUG_printf(" -- not finite --> %s\n", + (n_maybe_swaps <= 1) ? "goto maybe_swap" : "give up"); + if(n_maybe_swaps <= 1) + goto maybe_swap; + /* else was 'break;' ... + but rather give up returning NaN directly as in "normal scale" Newton */ + ML_WARNING(ME_DOMAIN, ""); + qb[0] = qb[1] = ML_NAN; return; + } if (i_pb >= n_N && w * wprev <= 0.) prev = fmax2(fabs(adj),fpu); - R_ifDEBUG_printf( - "N(i=%2d): u=%#20.16g, pb(e^u)=%#15.9g, w=%#15.9g, %s prev=%g,", - i_pb, u, y, w, - (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev); + R_ifDEBUG_printf(", %s prev=%g,", + (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev); g = 1; for (i_inn=0; i_inn < 1000; i_inn++) { adj = g * w; @@ -459,7 +501,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, // (cancellation in (u_n -u) => may differ from adj: double D = fmin2(fabs(adj), fabs(u_n - u)); /* R_ifDEBUG_printf(" delta(u)=%g\n", u_n - u); */ - R_ifDEBUG_printf(" it{in}=%d, delta(u)=%9.3g, D/|.|=%.3g\n", + R_ifDEBUG_printf(" it{in}=%d, d.(u)=%6.3g, D/|.|=%.3g\n", i_inn, u_n - u, D/fabs(u_n + u)); if (D <= 4e-16 * fabs(u_n + u)) goto L_converged; @@ -473,16 +515,6 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, for (i_pb=0; i_pb < 1000; i_pb++) { y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p); // delta{y} : d_y = y - (log_p ? la : a); -#ifdef IEEE_754 - if(!R_FINITE(y) && !(log_p && y == ML_NEGINF))// y = -Inf is ok if(log_p) -#else - if (errno) -#endif - { // ML_WARN_return_NAN : - ML_WARNING(ME_DOMAIN, ""); - qb[0] = qb[1] = ML_NAN; return; - } - /* w := Newton step size (F(.) - a) / F'(.) or, * -- log: (lF - la) / (F' / F) = exp(lF) * (lF - la) / F' @@ -490,10 +522,25 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, w = log_p ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta)) : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta)); + if(!R_FINITE(w)) { + // what we should do is ==> go back, get better starting value + R_ifDEBUG_printf( + "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g -- is not finite --> %s\n", + i_pb, xinbta, y, w, (n_maybe_swaps <= 2) ? "goto maybe_swap" : "give up"); + if(n_maybe_swaps <= 2) { + if(!log_p && n_maybe_swaps == 2) use_log_x = TRUE; // try now + if(!log_p || n_maybe_swaps <= 1) + goto maybe_swap ; + } + /* else was 'break;' ... + but rather give up returning NaN directly as in "normal scale" Newton */ + ML_WARNING(ME_DOMAIN, ""); + qb[0] = qb[1] = ML_NAN; return; + } if (i_pb >= n_N && w * wprev <= 0.) prev = fmax2(fabs(adj),fpu); R_ifDEBUG_printf( - "N(i=%2d): x0=%#17.15g, pb(x0)=%#15.9g, w=%#15.9g, %s prev=%g,", + "N(i=%2d): x0=%#19.15g, pb(x0)=%#15.11g, w=%#12.7g, %s prev=%g,", i_pb, xinbta, y, w, (i_pb >= n_N && w * wprev <= 0.) ? "new" : "old", prev); g = 1; @@ -530,7 +577,7 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, ML_WARNING(ME_PRECISION, "qbeta"); L_converged: - log_ = log_p || use_log_x; // only for printing + log_ = log_p || use_log_x; R_ifDEBUG_printf(" %s: Final delta(y) = %g%s\n", warned ? "_NO_ convergence" : "converged", y - (log_ ? la : a), (log_ ? " (log_)" : "")); @@ -553,7 +600,8 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, (log_ ? ", log_" : ""), fabs(y - (log_ ? la : a))); } L_return: - if(give_log_q) { // ==> use_log_x , too + // use if (use_log_x) u_n else tx {and u_n is on log scale} + if(give_log_q) { // {currently not used from R's qbeta()} ==> use_log_x , too if(!use_log_x) // (see if claim above is true) MATHLIB_WARNING( "qbeta() L_return, u_n=%g; give_log_q=TRUE but use_log_x=FALSE -- please report!", @@ -569,16 +617,22 @@ qbeta_raw(double alpha, double p, double q, int lower_tail, int log_p, if(add_N_step) { /* add one last Newton step on original x scale, e.g., for qbeta(2^-98, 0.125, 2^-96) */ - xinbta = exp(u_n); + if(u_n != 1.) // u_n has been computed above + xinbta = exp(u_n); y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE, log_p); w = log_p ? (y - la) * exp(y + logbeta + r * log(xinbta) + t * log1p(-xinbta)) : (y - a) * exp( logbeta + r * log(xinbta) + t * log1p(-xinbta)); - tx = xinbta - w; - R_ifDEBUG_printf(" Final Newton correction(non-log scale):\n" - // \n xinbta=%.16g - " xinbta=%.16g, y=%g, w=-Delta(x)=%g. \n=> new x=%.16g\n", - xinbta, y, w, tx); + if(R_FINITE(w)) { + tx = xinbta - w; + R_ifDEBUG_printf(" Final Newton correction(non-log scale):\n" + " xinbta=%.16g, y=%g, w=-Delta(x)=%g. \n=> new x=%.16g\n", + xinbta, y, w, tx); + } else { // Newton step w cannot be used + R_ifDEBUG_printf(" Final Newton correction(non-log scale), canNOT use step 'w':\n" + " xinbta=%.16g, y=%g, w=-Delta(x)=%g\n=> keep x=xinbta\n", xinbta, y, w); + tx = xinbta; + } } else { if(swap_tail) { qb[0] = -expm1(u_n); qb[1] = exp (u_n); diff --git a/src/qbinom.c b/src/qbinom.c index cbc1742..ddda3e4 100644 --- a/src/qbinom.c +++ b/src/qbinom.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 1999-2021 The R Core Team + * Copyright (C) 2003-2021 The R Foundation * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2009 The R Core Team - * Copyright (C) 2003-2009 The R Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -33,41 +33,24 @@ #include "nmath.h" #include "dpq.h" -static double -do_search(double y, double *z, double p, double n, double pr, double incr) -{ - if(*z >= p) { - /* search to the left */ -#ifdef DEBUG_qbinom - REprintf("\tnew z=%7g >= p = %7g --> search to left (y--) ..\n", z,p); -#endif - for(;;) { - double newz; - if(y == 0 || - (newz = pbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p) - return y; - y = fmax2(0, y - incr); - *z = newz; - } - } - else { /* search to the right */ #ifdef DEBUG_qbinom - REprintf("\tnew z=%7g < p = %7g --> search to right (y++) ..\n", z,p); +# define R_DBG_printf(...) REprintf(__VA_ARGS__) +#else +# define R_DBG_printf(...) #endif - for(;;) { - y = fmin2(y + incr, n); - if(y == n || - (*z = pbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p) - return y; - } - } -} +#define _thisDIST_ binom +#define _dist_PARS_DECL_ double n, double pr +#define _dist_PARS_ n, pr +#define _dist_MAX_y n +// === Binomial Y <= n + +#include "qDiscrete_search.h" +// ------------------> do_search() and all called by q_DISCRETE_*() below + double qbinom(double p, double n, double pr, int lower_tail, int log_p) { - double q, mu, sigma, gamma, z, y; - #ifdef IEEE_754 if (ISNAN(p) || ISNAN(n) || ISNAN(pr)) return p + n + pr; @@ -78,57 +61,26 @@ double qbinom(double p, double n, double pr, int lower_tail, int log_p) if(!R_FINITE(p) && !log_p) ML_WARN_return_NAN; - if(n != floor(n + 0.5)) ML_WARN_return_NAN; + n = R_forceint(n); + if (pr < 0 || pr > 1 || n < 0) ML_WARN_return_NAN; R_Q_P01_boundaries(p, 0, n); if (pr == 0. || n == 0) return 0.; + if (pr == 1.) return n; /* covers the full range of the distribution */ - q = 1 - pr; - if(q == 0.) return n; /* covers the full range of the distribution */ - mu = n * pr; - sigma = sqrt(n * pr * q); - gamma = (q - pr) / sigma; - -#ifdef DEBUG_qbinom - REprintf("qbinom(p=%7g, n=%g, pr=%7g, l.t.=%d, log=%d): sigm=%g, gam=%g\n", - p,n,pr, lower_tail, log_p, sigma, gamma); -#endif - /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- - * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ - if(!lower_tail || log_p) { - p = R_DT_qIv(p); /* need check again (cancellation!): */ - if (p == 0.) return 0.; - if (p == 1.) return n; - } - /* temporary hack --- FIXME --- */ - if (p + 1.01*DBL_EPSILON >= 1.) return n; - - /* y := approx.value (Cornish-Fisher expansion) : */ - z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); - y = floor(mu + sigma * (z + gamma * (z*z - 1) / 6) + 0.5); - - if(y > n) /* way off */ y = n; - -#ifdef DEBUG_qbinom - REprintf(" new (p,1-p)=(%7g,%7g), z=qnorm(..)=%7g, y=%5g\n", p, 1-p, z, y); -#endif - z = pbinom(y, n, pr, /*lower_tail*/TRUE, /*log_p*/FALSE); + // (NB: unavoidable cancellation for pr ~= 1) + double + q = 1 - pr, + mu = n * pr, + sigma = sqrt(n * pr * q), + gamma = (q - pr) / sigma; - /* fuzz to ensure left continuity: */ - p *= 1 - 64*DBL_EPSILON; + R_DBG_printf("qbinom(p=%.12g, n=%.15g, pr=%.7g, l.t.=%d, log=%d): sigma=%g, gamma=%g;\n", + p, n,pr, lower_tail, log_p, sigma, gamma); - if(n < 1e5) return do_search(y, &z, p, n, pr, 1); - /* Otherwise be a bit cleverer in the search */ - { - double incr = floor(n * 0.001), oldincr; - do { - oldincr = incr; - y = do_search(y, &z, p, n, pr, incr); - incr = fmax2(1, floor(incr/100)); - } while(oldincr > 1 && incr > n*1e-15); - return y; - } + q_DISCRETE_01_CHECKS(); + q_DISCRETE_BODY(); } diff --git a/src/qf.c b/src/qf.c index c99be89..f22124e 100644 --- a/src/qf.c +++ b/src/qf.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions - * Copyright (C) 1998 Ross Ihaka * Copyright (C) 2000-2015 The R Core Team * Copyright (C) 2005 The R Foundation + * Copyright (C) 1998 Ross Ihaka * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -43,10 +43,10 @@ double qf(double p, double df1, double df2, int lower_tail, int log_p) if (df1 <= df2 && df2 > 4e5) { if(!R_FINITE(df1)) /* df1 == df2 == Inf : */ return 1.; - /* else */ + /* else value for df2 == Inf : */ return qchisq(p, df1, lower_tail, log_p) / df1; } - if (df1 > 4e5) { /* and so df2 < df1 */ + else if (df1 > 4e5) { /* and so df2 < df1 -- return value for df1 == Inf */ return df2 / qchisq(p, df2, !lower_tail, log_p); } diff --git a/src/qnbinom.c b/src/qnbinom.c index b313ce5..2ba5554 100644 --- a/src/qnbinom.c +++ b/src/qnbinom.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2021 The R Core Team + * Copyright (C) 2005-2021 The R Foundation * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2016 The R Core Team - * Copyright (C) 2005-2016 The R Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -26,7 +26,8 @@ * * DESCRIPTION * - * The quantile function of the negative binomial distribution. + * The quantile function of the negative binomial distribution, + * for the (size, prob) parametrizations * * NOTES * @@ -44,31 +45,23 @@ #include "nmath.h" #include "dpq.h" -static double -do_search(double y, double *z, double p, double n, double pr, double incr) -{ - if(*z >= p) { /* search to the left */ - for(;;) { - if(y == 0 || - (*z = pnbinom(y - incr, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) < p) - return y; - y = fmax2(0, y - incr); - } - } - else { /* search to the right */ - for(;;) { - y = y + incr; - if((*z = pnbinom(y, n, pr, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p) - return y; - } - } -} +/*-------- DEBUGGING ------------- make CFLAGS='-DDEBUG_qnbinom ...' + */ +#ifdef DEBUG_qnbinom +# define R_DBG_printf(...) REprintf(__VA_ARGS__) +#else +# define R_DBG_printf(...) +#endif +#define _thisDIST_ nbinom +#define _dist_PARS_DECL_ double size, double prob +#define _dist_PARS_ size, prob + +#include "qDiscrete_search.h" +// ------------------> do_search() and all called by q_DISCRETE_*() below double qnbinom(double p, double size, double prob, int lower_tail, int log_p) { - double P, Q, mu, sigma, gamma, z, y; - #ifdef IEEE_754 if (ISNAN(p) || ISNAN(size) || ISNAN(prob)) return p + size + prob; @@ -78,56 +71,22 @@ double qnbinom(double p, double size, double prob, int lower_tail, int log_p) prob == size/(size+mu) */ if (prob == 0 && size == 0) return 0; - if (prob <= 0 || prob > 1 || size < 0) ML_WARN_return_NAN; - if (prob == 1 || size == 0) return 0; R_Q_P01_boundaries(p, 0, ML_POSINF); - Q = 1.0 / prob; - P = (1.0 - prob) * Q; - mu = size * P; - sigma = sqrt(size * P * Q); - gamma = (Q + P)/sigma; - - /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- - * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ - if(!lower_tail || log_p) { - p = R_DT_qIv(p); /* need check again (cancellation!): */ - if (p == R_DT_0) return 0; - if (p == R_DT_1) return ML_POSINF; - } - /* temporary hack --- FIXME --- */ - if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF; - - /* y := approx.value (Cornish-Fisher expansion) : */ - z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); - y = R_forceint(mu + sigma * (z + gamma * (z*z - 1) / 6)); - - z = pnbinom(y, size, prob, /*lower_tail*/TRUE, /*log_p*/FALSE); + double + Q = 1.0 / prob, + P = (1.0 - prob) * Q, // = (1 - prob) / prob = Q - 1 + mu = size * P, + sigma = sqrt(size * P * Q), + gamma = (Q + P)/sigma; - /* fuzz to ensure left continuity: */ - p *= 1 - 64*DBL_EPSILON; + R_DBG_printf("qnbinom(p=%.12g, size=%.15g, prob=%g, l.t.=%d, log=%d):" + " mu=%g, sigma=%g, gamma=%g;\n", + p, size, prob, lower_tail, log_p, mu, sigma, gamma); - /* If the C-F value is not too large a simple search is OK */ - if(y < 1e5) return do_search(y, &z, p, size, prob, 1); - /* Otherwise be a bit cleverer in the search */ - { - double incr = floor(y * 0.001), oldincr; - do { - oldincr = incr; - y = do_search(y, &z, p, size, prob, incr); - incr = fmax2(1, floor(incr/100)); - } while(oldincr > 1 && incr > y*1e-15); - return y; - } -} - -double qnbinom_mu(double p, double size, double mu, int lower_tail, int log_p) -{ - if (size == ML_POSINF) // limit case: Poisson - return(qpois(p, mu, lower_tail, log_p)); -/* FIXME! Implement properly!! (not losing accuracy for very large size (prob ~= 1)*/ - return qnbinom(p, size, /* prob = */ size/(size+mu), lower_tail, log_p); + q_DISCRETE_01_CHECKS(); + q_DISCRETE_BODY(); } diff --git a/src/qnbinom_mu.c b/src/qnbinom_mu.c new file mode 100644 index 0000000..d1b8ca5 --- /dev/null +++ b/src/qnbinom_mu.c @@ -0,0 +1,83 @@ +/* + * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2021 The R Core Team + * Copyright (C) 2005-2021 The R Foundation + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, a copy is available at + * https://www.R-project.org/Licenses/ + * + * SYNOPSIS + * + * #include + * double qnbinom_mu(double p, double size, double mu, + * int lower_tail, int log_p) + * + * DESCRIPTION + * + * The quantile function of the negative binomial distribution, + * for the (size, mu) parametrizations + * + * METHOD + * + * Uses the Cornish-Fisher Expansion to include a skewness + * correction to a normal approximation. This gives an + * initial value which never seems to be off by more than + * 1 or 2. A search is then conducted of values close to + * this initial start point. + */ + +#include "nmath.h" +#include "dpq.h" + +#ifdef DEBUG_qnbinom +# define R_DBG_printf(...) REprintf(__VA_ARGS__) +#else +# define R_DBG_printf(...) +#endif + +#define _thisDIST_ nbinom_mu +#define _dist_PARS_DECL_ double size, double mu +#define _dist_PARS_ size, mu + +#include "qDiscrete_search.h" +// ------------------> do_search() and all called by q_DISCRETE_*() below + +double qnbinom_mu(double p, double size, double mu, int lower_tail, int log_p) +{ + if (size == ML_POSINF) // limit case: Poisson + return(qpois(p, mu, lower_tail, log_p)); + +#ifdef IEEE_754 + if (ISNAN(p) || ISNAN(size) || ISNAN(mu)) + return p + size + mu; +#endif + + if (mu == 0 || size == 0) return 0; + if (mu < 0 || size < 0) ML_WARN_return_NAN; + + R_Q_P01_boundaries(p, 0, ML_POSINF); + + double + Q = 1 + mu/size, // (size+mu)/size = 1 / prob + P = mu/size, // = (1 - prob) * Q = (1 - prob) / prob = Q - 1 + sigma = sqrt(size * P * Q), + gamma = (Q + P)/sigma; + + R_DBG_printf("qnbinom_mu(p=%.12g, size=%.15g, mu=%g, l.t.=%d, log=%d):" + " mu=%g, sigma=%g, gamma=%g;\n", + p, size, mu, lower_tail, log_p, mu, sigma, gamma); + + q_DISCRETE_01_CHECKS(); + q_DISCRETE_BODY(); +} diff --git a/src/qnchisq.c b/src/qnchisq.c index a7e2937..6893bfa 100644 --- a/src/qnchisq.c +++ b/src/qnchisq.c @@ -1,8 +1,8 @@ /* * R : A Computer Language for Statistical Data Analysis - * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka - * Copyright (C) 2000-2008 The R Core Team + * Copyright (C) 2000--2020 The R Core Team * Copyright (C) 2004 The R Foundation + * Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -46,8 +46,9 @@ double qnchisq(double p, double df, double ncp, int lower_tail, int log_p) R_Q_P01_boundaries(p, 0, ML_POSINF); - pp = R_D_qIv(p); - if(pp > 1 - DBL_EPSILON) return lower_tail ? ML_POSINF : 0.0; + pp = R_D_qIv(p); // exp(p) iff log_p + if(pp > 1 - DBL_EPSILON) + return lower_tail ? ML_POSINF : 0.0; // early under/over flow iff log_p (FIXME) /* Invert pnchisq(.) : * 1. finding an upper and lower bound */ @@ -59,7 +60,7 @@ double qnchisq(double p, double df, double ncp, int lower_tail, int log_p) c = (df + 3*ncp)/(df + 2*ncp); ff = (df + 2 * ncp)/(c*c); ux = b + c * qchisq(p, ff, lower_tail, log_p); - if(ux < 0) ux = 1; + if(ux <= 0.) ux = 1; ux0 = ux; } diff --git a/src/qnorm.c b/src/qnorm.c index 347f1b1..1e08280 100644 --- a/src/qnorm.c +++ b/src/qnorm.c @@ -1,7 +1,7 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2000--2020 The R Core Team * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000--2005 The R Core Team * based on AS 111 (C) 1977 Royal Statistical Society * and on AS 241 (C) 1988 Royal Statistical Society * @@ -81,8 +81,8 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p) (original fortran code used PARAMETER(..) for the coefficients and provided hash codes for checking them...) */ - if (fabs(q) <= .425) {/* 0.075 <= p <= 0.925 */ - r = .180625 - q * q; + if (fabs(q) <= .425) {/* |p~ - 0.5| <= .425 <==> 0.075 <= p~ <= 0.925 */ + r = .180625 - q * q; // = .425^2 - q^2 >= 0 val = q * (((((((r * 2509.0809287301226727 + 33430.575583588128105) * r + 67265.770927008700853) * r + @@ -94,22 +94,18 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p) 21213.794301586595867) * r + 5394.1960214247511077) * r + 687.1870074920579083) * r + 42.313330701600911252) * r + 1.); } - else { /* closer than 0.075 from {0,1} boundary */ - - /* r = min(p, 1-p) < 0.075 */ - if (q > 0) - r = R_DT_CIv(p);/* 1-p */ - else - r = p_;/* = R_DT_Iv(p) ^= p */ - - r = sqrt(- ((log_p && - ((lower_tail && q <= 0) || (!lower_tail && q > 0))) ? - p : /* else */ log(r))); - /* r = sqrt(-log(r)) <==> min(p, 1-p) = exp( - r^2 ) */ + else { /* closer than 0.075 from {0,1} boundary : + * r := log(p~); p~ = min(p, 1-p) < 0.075 : */ + if(log_p && ((lower_tail && q <= 0) || (!lower_tail && q > 0))) { + r = p; + } else { + r = log( (q > 0) ? R_DT_CIv(p) /* 1-p */ : p_ /* = R_DT_Iv(p) ^= p */); + } + // r = sqrt( - log(min(p,1-p)) ) <==> min(p, 1-p) = exp( - r^2 ) : + r = sqrt(-r); #ifdef DEBUG_qnorm REprintf("\t close to 0 or 1: r = %7g\n", r); #endif - if (r <= 5.) { /* <==> min(p,1-p) >= exp(-25) ~= 1.3888e-11 */ r += -1.6; val = (((((((r * 7.7454501427834140764e-4 + @@ -125,7 +121,12 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p) r + 1.6763848301838038494) * r + 2.05319162663775882187) * r + 1.); } - else { /* very close to 0 or 1 */ + else if(r >= 816) { // p is *extremly* close to 0 or 1 - only possibly when log_p =TRUE + // Using the asymptotical formula -- is *not* optimal but uniformly better than branch below + val = r * M_SQRT2; + } + else { // p is very close to 0 or 1: r > 5 <==> min(p,1-p) < exp(-25) = 1.3888..e-11 + // Wichura, p.478: minimax rational approx R_3(t) is for 5 <= t <= 27 (t :== r) r += -5.; val = (((((((r * 2.01033439929228813265e-7 + 2.71155556874348757815e-5) * r + @@ -147,6 +148,3 @@ double qnorm5(double p, double mu, double sigma, int lower_tail, int log_p) } return mu + sigma * val; } - - - diff --git a/src/qpois.c b/src/qpois.c index f911b54..e6f87d3 100644 --- a/src/qpois.c +++ b/src/qpois.c @@ -1,7 +1,7 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 1999-2021 The R Core Team * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2016 The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -33,31 +33,22 @@ #include "nmath.h" #include "dpq.h" -static double -do_search(double y, double *z, double p, double lambda, double incr) -{ - if(*z >= p) { - /* search to the left */ - for(;;) { - if(y == 0 || - (*z = ppois(y - incr, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) < p) - return y; - y = fmax2(0, y - incr); - } - } - else { /* search to the right */ +#ifdef DEBUG_qpois +# define R_DBG_printf(...) REprintf(__VA_ARGS__) +#else +# define R_DBG_printf(...) +#endif - for(;;) { - y = y + incr; - if((*z = ppois(y, lambda, /*l._t.*/TRUE, /*log_p*/FALSE)) >= p) - return y; - } - } -} + +#define _thisDIST_ pois +#define _dist_PARS_DECL_ double lambda +#define _dist_PARS_ lambda + +#include "qDiscrete_search.h" +// ------------------> do_search() and all called by q_DISCRETE_*() below double qpois(double p, double lambda, int lower_tail, int log_p) { - double mu, sigma, gamma, z, y; #ifdef IEEE_754 if (ISNAN(p) || ISNAN(lambda)) return p + lambda; @@ -70,40 +61,16 @@ double qpois(double p, double lambda, int lower_tail, int log_p) if(p == R_DT_0) return 0; if(p == R_DT_1) return ML_POSINF; - mu = lambda; - sigma = sqrt(lambda); - /* gamma = sigma; PR#8058 should be kurtosis which is mu^-0.5 */ - gamma = 1.0/sigma; - - /* Note : "same" code in qpois.c, qbinom.c, qnbinom.c -- - * FIXME: This is far from optimal [cancellation for p ~= 1, etc]: */ - if(!lower_tail || log_p) { - p = R_DT_qIv(p); /* need check again (cancellation!): */ - if (p == 0.) return 0; - if (p == 1.) return ML_POSINF; - } - /* temporary hack --- FIXME --- */ - if (p + 1.01*DBL_EPSILON >= 1.) return ML_POSINF; - - /* y := approx.value (Cornish-Fisher expansion) : */ - z = qnorm(p, 0., 1., /*lower_tail*/TRUE, /*log_p*/FALSE); - y = nearbyint(mu + sigma * (z + gamma * (z*z - 1) / 6)); - - z = ppois(y, lambda, /*lower_tail*/TRUE, /*log_p*/FALSE); + double + mu = lambda, + sigma = sqrt(lambda), + // had gamma = sigma; PR#8058 should be kurtosis which is mu^-0.5 = 1/sigma + gamma = 1.0/sigma; - /* fuzz to ensure left continuity; 1 - 1e-7 may lose too much : */ - p *= 1 - 64*DBL_EPSILON; + R_DBG_printf("qpois(p=%.12g, lambda=%.15g, l.t.=%d, log=%d):" + " mu=%g, sigma=%g, gamma=%g;\n", + p, lambda, lower_tail, log_p, mu, sigma, gamma); - /* If the mean is not too large a simple search is OK */ - if(lambda < 1e5) return do_search(y, &z, p, lambda, 1); - /* Otherwise be a bit cleverer in the search */ - { - double incr = floor(y * 0.001), oldincr; - do { - oldincr = incr; - y = do_search(y, &z, p, lambda, incr); - incr = fmax2(1, floor(incr/100)); - } while(oldincr > 1 && incr > lambda*1e-15); - return y; - } + // never "needed" here (FIXME?): q_DISCRETE_01_CHECKS(); + q_DISCRETE_BODY(); } diff --git a/src/qt.c b/src/qt.c index 915e724..b74eb79 100644 --- a/src/qt.c +++ b/src/qt.c @@ -1,8 +1,8 @@ /* * Mathlib : A C Library of Special Functions + * Copyright (C) 2000-2022 The R Core Team + * Copyright (C) 2003-2022 The R Foundation * Copyright (C) 1998 Ross Ihaka - * Copyright (C) 2000-2013 The R Core Team - * Copyright (C) 2003-2013 The R Foundation * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -143,8 +143,9 @@ double qt(double p, double ndf, int lower_tail, int log_p) b = 48 / (a * a), c = ((20700 * a / b - 98) * a - 16) * a + 96.36, d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * M_PI_2) * ndf; - - Rboolean P_ok1 = P > DBL_MIN || !log_p, P_ok = P_ok1; + Rboolean + P_ok1 = P > DBL_MIN || !log_p, + P_ok = P_ok1; // when true (after check below), use "normal scale": log_p=FALSE if(P_ok1) { y = pow(d * P, 2.0 / ndf); P_ok = (y >= DBL_EPSILON); @@ -190,16 +191,26 @@ double qt(double p, double ndf, int lower_tail, int log_p) * Probably also improvable when lower_tail = FALSE */ if(P_ok1) { + double M = fabs(sqrt(DBL_MAX/2.) - ndf); int it=0; while(it++ < 10 && (y = dt(q, ndf, FALSE)) > 0 && R_FINITE(x = (pt(q, ndf, FALSE, FALSE) - P/2) / y) && - fabs(x) > 1e-14*fabs(q)) + fabs(x) > 1e-14*fabs(q)) { /* Newton (=Taylor 1 term): * q += x; * Taylor 2-term : */ - q += x * (1. + x * q * (ndf + 1) / (2 * (q * q + ndf))); + double F = (fabs(q) < M) ? + q * (ndf + 1) / (2 * (q * q + ndf)) : + (ndf + 1) / (2 * (q + ndf/q)), + del_q = x * (1. + x * F); + if(R_FINITE(del_q) && R_FINITE(q + del_q)) + q += del_q; + else if(R_FINITE(x) && R_FINITE(q + x)) + q += x; + else // FIXME?? if q+x = +/-Inf is *better* than q should use it + break; // cannot improve q with a Newton/Taylor step + } } } - if(neg) q = -q; - return q; + return neg ? -q : q; } diff --git a/src/rmultinom.c b/src/rmultinom.c index f14b195..8358e9c 100644 --- a/src/rmultinom.c +++ b/src/rmultinom.c @@ -75,7 +75,7 @@ void rmultinom(int n, double* prob, int K, int* rN) rN[k] = 0; } if(fabs((double)(p_tot - 1.)) > 1e-7) - MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"), + MATHLIB_ERROR(_("rbinom: probability sum should be 1, but is %g"), (double) p_tot); if (n == 0) return; if (K == 1 && p_tot == 0.) return;/* trivial border case: do as rbinom */ diff --git a/src/stirlerr.c b/src/stirlerr.c index 1beb3dc..c9b3cc6 100644 --- a/src/stirlerr.c +++ b/src/stirlerr.c @@ -4,7 +4,7 @@ * October 23, 2000. * * Merge in to R: - * Copyright (C) 2000, The R Core Team + * Copyright (C) 2000-2021, The R Core Team * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -41,8 +41,9 @@ * = log Gamma(n+1) - (n + 1/2) * log(n) + n - log(2*pi)/2 * * see also lgammacor() in ./lgammacor.c which computes almost the same! + * + * NB: stirlerr(n/2) is called from dt() *and* gamma(n/2) when n is integer and n/2 <= 50 */ - double attribute_hidden stirlerr(double n) { @@ -53,7 +54,7 @@ double attribute_hidden stirlerr(double n) #define S4 0.0008417508417508417508417508/* 1/1188 */ /* - error for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0. + exact values for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0. */ const static double sferr_halves[31] = { 0.0, /* n=0 - wrong, place holder only */