pnchisq_ss() gets new optional argument ssr = ss(....)for flexibility.
New gam1() calls mpfr-ized direct code or the double-only C code gam1d()
Add Abramowitz & Stegun, p.945, pbeta() approximations 26.5.20 and 26.5.21 as pbetaAS_eq20() and pbetaAS_eq21(). For “didactical” reasons, also provide the simple moment matching normal approximation pbetaNorm2(). For now, also low-level versions .pbeta.eq20(), ..., .pbetaSeq21().
BUG FIXES
ldexp({}, 1) or pow({}, 2) no longer segfault.
pbetaRv1() for verbose = TRUE and/ormpfr-number arguments.
dnbinomR() and dnbinom.mu(): fix vectorization logic, using ifelse().
pbetaRv1() for verbose = TRUE and/ormpfr-number arguments.
dnbinomR() and dnbinom.mu(): fix vectorization logic, using ifelse().
Version 0.6-0 (2025-07-08, svn r342)
NEW FEATURES
added functions dltgammaInc() and lgammaP11(), both from TOMS Algorithm 1006 by Rémy Abergel and Lionel Moisan (2020).
C-level logcf() now uses long double, i.e., R'sLDOUBLE, to possibly improve for smallish eps. logcf() also got an optional maxit argument. For the pure-R version, renamed the simplistic and slow one from logcfR() to logcfR_vec() and the smarter faster one from logcfR.() to logcfR(). They have given identical results in test cases, and logcfR() should really be the one used practically.
bd0C() has new version = "R_2025_0510" where asbd0() has just been updated to the latest version also in R-devel svn revisions 88190,1,3 and 88208. New bd0_p1l1(), also as version in dpois_raw();bd0(), and the now four bd0_<meth>() functions work with x=Inf and properly vectorize.
The log1pmx vignette now shows that for dpois(x, 48) at least, using the two parts ofebd0() gives clearly the most accurate values (used in R since ca. 2022). On the other hand, new evaluations at the end of ‘tests/bd0-tst.R’ indicate that increasing the (only inDPQ) optional cutoff from delta = 0.1 to aboutdelta = 0.4 gets much more accurate bd0(x, np) values, notably when |x - np| is smallish for relatively large x and np whereebd0() is clearly less accurate.
BUG FIXES
rchkPROTECT fix in ‘bd0.c’.
Strictly not comparing Rboolean withNA_LOGICAL.
log1pmx(Inf) is correct (= -Inf) now.
dpois_raw(1, {}) now works (returning numeric(0)).
MISC
‘tests/pow-tst.R’ also work on "bizarre" platforms whereosVersion is NULL.
‘tests/stirlerr-tst.R’ uses many more if(doExtras) ... to run faster without “extras”.
Version 0.5-9 (2024-08-23, svn r333)
NEW FEATURES
dbinom_raw() gets a new more accurateversion = "R4.4", using the new pow1p() function, from R-devel, i.e., R 4.4.0 and newer.
stirlerr_simpl() gets argument version, a string specifying the exact function to be used, currently with options ("R3", "lgamma1p", "MM2", "n0"). The first being the direct formula hardwired in R up to versions 4.3.z, "lgamma1p" is to uselgamma1p(n) instead of lgamma(n+1), whereas"MM2" is a version using lgamma(n) and slightly more accurate for n \in [1, 5], approximately. stirlerr() correspondingly gets argument direct.ver to be passed to stirlerr_simpl(., version=*) in the “direct” formula case, when n < cutoffs[1].
stirlerr() gets argument order allowing currently up to order 20 terms (constants S0..S19) in Stirling's series (and not using cutoffs or direct formula).
New stirlerrC() interfacing to C code which should “parallelize” R's own C level stirlerr(). More documentation and tests; experiments added in ‘../Misc/stirlerr-trms.R’; nice explorations to ‘../tests/stirlerr-tst.R’.
New functions pow(x,y) and pow_di(x,y) computingx^y as in R's C API ‘Rmathlib’.
dgamma.R(x, *) gets new optional argument dpois_r_args.
New function expm1x() providing exp(x) - 1 - xnumerically accurately also for small |x|.
New experimental Ixpq() implementing Gil et al.(2023)'s direct incomplete Beta function, i.e., pbeta() version.TODO: make it work with mpfr() !
New function pntGST23_T6() and pntGST23_1()implementing Gil et al.(2023)'s first asymptotic formula (Theorem 6), and based on Temme's Maple code for their direct computation (1). pnt*(): CDF of the non-central t-distribution.
New function pntVW13() implementing Viktor Witkovský(2013)'s algorithm, from his published Matlab code.
Data set pt_Witkovsky_Tab1 containing (the free software subset of) Viktor Witkovský(2013)'s Table 1 ofpt(x, df=nu, ncp=delta) examples; updated with octave results for both the original 2013 matlab code and the 2022 corrected one.
lgamma1p_series(x, k), the Taylor series approximations for lgamma1p(x) := \log \Gamma(1+x) are now available up to k = 15, also for mpfr() numbers x.
New gam1d() and gamln1() providing more TOMS 708 auxiliary routines.
BUG FIXES
fix format fixes for when 'long != long long'.
dgamma.R(x, ..) now does vectorize in x, as documented.
minor adjustment of pntR()'s default use.pnorm.
‘src/bd0.c’ needed one more PROTECT().
Misc
new ‘tests/dt-ex.R’ exploring dt() accuracy issues, mostly historical, but newly for df=\nu \ll 1.
Version 0.5-8 (2023-11-30, svn r278)
BUG FIXES
ebd0() and C level ebd0C() now work with long vectors and hence return data.frame(yh = *, yl = *) instead of 2-row matrix.
format fixes: "%lld" and (long long) notably forR_xlen_t.
Deprecated phyperBinMolenaar(); it has been synonym to phyperBinMolenaar.1().
Fix LaTeX warnings/errors by removing amazingly many\cr in ‘man/*.Rd’ files.
Version 0.5-5 (2023-06-23, svn r264)
BUG FIXES
qbinomR(p, *), qnbinomR(p, *), and qpoisR(p, *): Fix R_Q_P01_boundaries(p, ..) for trivial boundaries, e.g.,p \in \{0,1\} for log.p=FALSE.
Misc
A new flang Fortran compiler does not know outdatedderfc(), but needs current standard erfc().
Version 0.5-4 (2023-04-12, svn r261)
NEW FEATURES
New R function qnormCappr().
New R function rexpm1() from TOMS 708.
Provide also qnormR(*, version = "1.0.x") and "1.0_noN", for history's sake.
New R function rlog1() from TOMS 708.
BUG FIXES
qnormUappr() & qnormUappr6() were _wrongly_negative for lp in [-.693147, 0) (corresponding to p > 1/2).
Misc
Tweaks to the ‘qnorm-asymp’ vignette, related to the JSS submission.
‘tests/pqnorm_extreme.R’ plots more about thepnormAsymp() relative errors.
Apple clang 14.0.3 needs looser tolerance for logcf{R}()and dbinom_raw() comparisons.
Minimally document the internal.D_*() and.DT_*() utilities, as we provide them to experts and e.g., package DPQmpfr.
Version 0.5-3 (2022-12-01, svn r240)
NEW FEATURES
New R functions qntR() and qtU(),Vectorize()d versions of qntR1() and qtU1(), respectively.
New qtR() <- Vectorize(qntR1, *); pure Rimplementation of R's Mathlib C-level qt(), but additionally allowing tweaks, used for fixing R's bug PR#18630. Added (optional) log-scale Newton steps via logNewton = log.p, needing more tests.
New qtNappr() – from the remark in R's ‘nmath/qt.c’ about very large df approximation, now with all 4 large-df terms from Abramowitz & Stegun's (26.7.5).
New gammaVer() to exemplify R's (partly historical) versions of gamma() implementations.
New qnormUappr6(), providing the ‘6 coefficients’-approximation of A. & S. to qnorm().
qnormR() gets new version = "2022-08-04" which uses MM's newly derived asymptotic approximations to qnorm().
New qnormAsymp() providing new asymptotic approximations to qnorm().
Added vignette “Asymptotic Tail Formulas for Gaussian Quantiles”.
New chebychevPoly() (and auxiliaries) for evaluation of Chebyshev polynomials.
Version 0.5-2 (2022-06-08, svn r212)
NEW FEATURES
qbetaAppr() (and the qbetaAppr.{1,3}()auxiliaries) now get a log.p argument.
all(?) qbeta.Appr*() functions now accept alower.tail argument (with default TRUE).
New R function dpsifn(), as interface to R's C APIdpsifn(), the workhorse of all R's psigamma()functions.
dpois_raw(x, lambda, *) gets new cutoff small.x__lambdato not use bd0() or ebd0() when x is much smaller than lambda.
Simple functions, mostly for didactical and comparison purposes,dpois_simpl(), dpois_simpl0() andstirlerr_simpl(), all of which use the “simple” direct formula which leads to numerical cancellation typically.
New R function bpser() interfacing to the C function of the same in R's Mathlib pbeta() or bratio() in file ‘nmath/toms708.c’.
BUG FIXES
Updated URL of Abramowitz & Stegun.
Fixed \ escape in one ‘man/*.Rd’
Misc
Our ‘tests/*.R’ no longer rely on the Matrix‘test-tools.R’ collection but rather use our own.
In ‘tests/chisq-nonc-ex.R’ (and ‘DESCRIPTION’), we've replaced akima by interp as the latter has a FOSS licence.
Version 0.5-1 (2021-12-10, svn r197)
NEW FEATURES
bd0(x, np, delta, ..) now usesif (|x-np| <= delta * (x+np)) (less-equal instead of strictly less), and hence setting delta = 0 is now allowed, using the series expansion only for x == np, useful e.g., for the case of highly accurate mpfr-numbers. Similarly, in log1pmx(x, .., eps2, ..), eps2 = 0 is now allowed.
BUG FIXES
consistency between ebd0() and ebd0C(): multiplication * e coming last ("bug" did not show ..).
okLongDouble() should no longer fail on M1 mac.
ebd0C() now checking |yl| < 5.5 before returning; was horrendously wrong on Windows compiled with -mnative as long as it included (the default) -mfma (FMA:= Fused-Multiply-Add).
fix rchk issue: "need" PROTECT(.) while calling allocating lgamma1p().
Version 0.5-0 (2021-09-10, svn r183)
NEW FEATURES
new ebd0C() interface to C version; helps to fix lapsus in pure R version ebd0(): ebd0(x, M) for large M now checks for overflow inM/x, and notably the case f * 2^-(e+10) =: fg == Inf.
In bd0() work around underflow of (x-np)/(x+np).
new R functions frexp() and ldexp() for getting and setting base-2 representations of numbers, and new R function modf() to split number into integer and fractional part, all interfacing the C99 (math lib) standard functions of the same name.
logcfR.(x) has been “vectorized” in x even though it's iterative with different number of iterations for eachx[i], and is hence considerably faster when x is an"mpfr" vector. TODO: Consider renaming the two logcfR versions; at least neither is deprecated!
BUG FIXES
fixed embarrassing typo (two i's) bug in logcfR()
stirlerr(n) now also works when is.integer(n) andn is large enough for n * n to overflow (to NA).
fix typo/thinko in dnbinomR()
ditto in logcfR(*, trace=TRUE) iteration report.
dpois_raw(x, *) now works up to maximal x, now preventing overflow in previous 2*pi*x computation.
dnbinomR() and dnbinom.mu() fix for x > size when log=TRUE, notably for x >> size.
Version 0.4-4 (2021-05-22, svn r175)
NEW FEATURES
new bd0_*() versions of bd0(), based mainly on log1pmx().
REFACTORING
To be more modular, our ‘test-tools.R’ no longersource() those of Matrix.
BUG FIXES
‘src/bd0.c’ gave warning with some compiler settings with-Wself-assign.
‘test-tools.R’ readRDS_() thinko fixed.
Version 0.4-3 (2021-05-05, svn r166)
NEW FEATURES
Provide R functionslog1mexpC(), log1pexpC(), log1pmxC(), and lgamma1pC(), all interfacing to R's C API (‘Rmath.h’), aka ‘Rmathlib’.
New pnormAsymp() for asymptotic (typically upper tail, log scale) approximation of pnorm(). With Rmpfr, we can see how accurate these are in detail.
New dnbinomR() finding better code for R itself,dnbinom.mu() and dbinom_raw(); also new utility functions bd0() and stirlerr(), both vectorized, and also working with "mpfr"-numbers, such that dpois_raw()now does so, too. Additionally, an “extended” version of bd0() calledebd0() in pure R, where a C version was proposed by Morten Welinder in R's bugzilla, PR#15628. Experimentally, also provide p1l1() and its Taylor series approximations p1l1ser() which could be employed forbd0() / ebd0() instead of the current algorithms.
Several C level utilities to be .Call()ed from R, from R's mathlib, ldexp() and frexp() even from C math.
Pure R implementations (somewhat experimental) of correspondingR Mathlib C functions: qbinomR(), qnbinomR(), and qpoisR()each with several tuning parameters for the underlying algorithm, notably the root-finding parts.
newly, logcf() now based on C code, perfectly vectorizes; the pure R version, now called logcfR(x, *)currently still runs the iterations simultaneously for all 'x[i]' and hence convergence and rescaling happen by “group decision”, efficiently but undesirable for strict comparisons. logcfR(x, *) and log1pmx(x, *) now also work for "mpfr"-numbers x, and log1pmx() gets optionaleps2 = 1e-2 and minL1 = -0.791 arguments, the defaults of which may change, as I think at least the latter to not be perfect.
Now lb_chiAsymp(nu, order) works up to order 8.
Provide the first parts of a new vignette (‘../vignettes/log1pmx-etc.Rnw’) on log1pmx(),bd0(), and stirlerr(), which should contain part of Loader(2000)'s report and new findings of improved bd0() andstirlerr() computations.
Version 0.4-2 (2020-11-07, svn r151)
NEW FEATURES
New functions pnormL*() and pnormU*() for (mathematically proven) lower and upper bound to pnorm(), notably also for investigation with log.p=TRUE.
qnormR(), implementing current R's qnorm() in pure R, notably with trace and version options.
Version 0.4-1 (2020-06-17, svn r145)
TESTING
Reorganize tests; notably to become less platform dependent.
Version 0.4-0 (2020-06-15, svn r137)
NEW FEATURES
Many new phyper*() functions and helpers for them, such as Bernoulli numbers Bern() and asymptoticlgammaAsymp().
Notably phyperR2() which is a pure R version of R's own (C code based) phyper().
Version 0.3-5 (2019-10-18, svn r131)
NEW FEATURES
pnbetaAS310() gained a useAS226 option.
New okLongDouble() function, notably for detecting that with a valgrinded version of R, long double C arithmetic seems to silently switch to (simple) double precision.
BUG FIXES
long double printing from C now should happen correctly on all platforms including 32- and 64-bit Windows.
Version 0.3-4 (2019-10-16, svn r125)
NEW FEATURES
more efficient (internally vectorized) dntJKBf(). Consequently, dntJKBf1() is deprecated now.
pntR() (and pntR1()) get new optionuse.pnorm (the default of which had been hard coded previously).
BUG FIXES
fix thinko in any_mpfr() and all_mpfr().
pnchisqRC()'s C code gave severe valgrind warnings; fixed printing of long double etc; also added special MinGW deal in Windows.
ppoisD() behaves differently in a valgrinded version; for now, reproducible only when using valgrind on non-valgrinded installed package.
Version 0.3-3 (2019-09-24, svn r110)
NEW FEATURES
Renamed (and generalized / merged) many functions to have less "." in names.
New pnbetaAS310() function interfacing to my_corrected_ C version of 'ASA 310' (2007).
New algdiv() function interfacing to the 'TOMS 708' version of our logQab().
New pl2curves() which generalizes (somewhat) previous function p.dnchiB().
Made newton() more flexible with new xMin andxMax arguments which are notably useful for q*()(quantile function) computations. Correspondingly replaced previous qchisq2() andqchisqA() with new qchisqN() (‘N’ewton).
new pnchisqRC(), a version of R's C level non-central chi-squared, with additional options.
“new” logspace.add() and logspace.sub().
“new” pnchisqT93() (plus its two auxiliaries), implementing Temme(1993)'s approximations.
“new” pnchisqBolKuz() and qnchisqBolKuz()implementing Bol'shev and Kuznetzov (1963)'s approximations.
“new” pnchi1sq() and pnchi2sq() with “exact” formulas for the special cases df=1 anddf=3.
simplified formula in dtWV().
BUG FIXES
qnchisqPearson(pp, df=DF, ncp=100) andqnchisqSankaran_d(*) no longer return NaN for very large DF = 1e200.
pnchisq() now also has default verbose = 0 as all other such functions with verbose (or trace[.lev]optional argument.
Version 0.3-0 [2018-08-28]
NEW FEATURES
Move many of my up to 15 years old DPQ computation utilities into a package, to become public, “bloggable”, etc.