News for Package mathr

Version 0.1.5 (2026-10-10)

• 0.1.4 was tagged but never sent to CRAN. Its changes, listed below, reach CRAN with this release.

• The check that 0.1.4 put in front of the functions that work a value at a time missed some arguments: p in Qf, Tol in VMmin, Plot in mlr, add in ellipse, GradTol, StepTol, FnTol and ItnLimit in Optim0, and n in Rnorm, and so in Rmvn, which hands it to Rnorm. Given a vector, each still stopped with one of R's own errors. Every argument that drives an if or a while is now checked, and n in Rmvn besides: 142 argument positions across 70 functions, among them n in trapez1 and N in romb, which gained an if with the fixes below. The help pages of the nine functions name the arguments added. Since Plot, add, Double in Dec2Bin and DecFlag in LOG are TRUE or FALSE, the error, and their help pages, now ask for "a single value" rather than "a single number".

• At a size of one, several functions failed. runifLC() at its default n = 1, SymInv, LDLT and evJacobi on a 1 x 1 matrix, and romb with N = 1 stopped with an error; trapez0 on a single point returned NA; trapez1 with n = 1 returned three times the trapezoidal value, and NashSVD on a one-column matrix returned 0 as its singular value. In most, a 2:n or 1:(n-1) that should be empty at n = 1 counted down instead (in trapez1 it added both end points a second time); NashSVD has no pair of columns to rotate then, and so never set the singular value. Each now handles the case, as Grad and Deriv2 already did, and the romb page no longer asks for N above 2.

• At n = 0, Rnorm returned FALSE, and Rgamma, Rgamma0 and Rbeta returned one deviate for most shapes, where rnorm(0) gives numeric(0). They all give numeric(0) now.

• SQRT(Inf) never returned: it handed Inf to NORM, whose halving loop cannot end on Inf, and NORM looped forever on a negative x as well. SQRT(Inf) is now Inf, with the warning EXP and LOG give; NORM stops with a message on a non-finite x and gives a negative fraction for a negative x.

• evJacobi stopped with "object 's' not found" on some symmetric matrices, diag(c(2, 1)) among them: a test on the rotation sine ran on a pair that needed no rotation, before any sine had been set. The test now runs only when a rotation is made.

• VMmin never returned when its search direction right after resetting the Hessian estimate was not downhill, as at a starting point where the gradient is zero: a misspelt SmallXCount kept its loop alive. Optim0 stopped with "object 'RelStep' not found" whenever it ended on its first iteration; RelStep and DelF are NA then.

• MachEps3 and MachEps4 with MaxIter = 0 stopped with "object 'CurEps' not found"; they return their first trial value, 0.5 and -0.5, with iteration 0. Dec2Bin(0) took the logarithm of 0 and returned NaN bits; it returns all zeros, the pattern of +0, and for -0 the same with the sign bit set. Run.test returned a p-value of 2 when no nonzero value was left; it returns NA.

• Given NA, EXP, LOG and SQRT stopped with R's "missing value where TRUE/FALSE needed", and SQRT did the same for NaN. They now return NA for NA, as exp, log and sqrt do, and SQRT returns NaN for NaN with the warning EXP and LOG already give.

• mlr no longer calls anyNA, which first appeared in R 3.1.0, above the R (>= 3.0.0) that DESCRIPTION declares.

• Source comments that described methods the code does not use are corrected. erfccheb, betaiapprox and gammpapprox forward to erfc, betai and gammp or gammq; erfc goes through gammq, not a Cephes rational approximation; and what the EXP guards keep from overflowing is g*P(y) and Q(y), not y = g*g.

• The Ref line on the betacf, betai, erf, erfc, gammln, gammp, gammq, gcf, gser, inverfc, invbetai and invgammp help pages lacked a full stop and ran into the sentence after it.

• GAMMA and LGAMMA hold the Lanczos coefficients the GNU Scientific Library uses, which only the gammln page said; their own pages now say it too. gammpapprox and betaiapprox are no longer titled as quadrature approximations, which in this package they never were.

• The 0.1.4 entries below are corrected where they were wrong. The EXP error figures were relative errors in units of .Machine$double.eps, not ulp, and the shortcut was taken on every call whose result is finite and nonzero, not on every call. The old EXP also returned Inf at log(.Machine$double.xmax) itself. Its NaN began near 2e78 rather than 1e78, came back instead of 0 for very large negative arguments as well, and was nearly but not quite total past 1e82. The overflow warning below 0.9 * .Machine$double.xmax is new in 0.1.4, not kept from before, and the 0.1.4 argument check did not reach every argument it said it did. Smaller corrections fix the list of bare titles and, under 0.1.0, which functions carry an attribution.

Version 0.1.4 (2026-09-20)

• EXP was not running its own algorithm. The range reduction left xexp = x/ln2 unrounded, so the reduced argument g never grew past about 1e-13, the abs(g) < sqrt(eps) shortcut was taken on every call whose result is finite and nonzero, and the rational approximation under it was dead code there. What came back was 2^(x/ln2) and nothing else, up to 223 eps off in relative terms (about 444 ulp) at the ends of the range, and Inf at log(.Machine$double.xmax) itself, where x/ln2 rounds to 1024. xexp is now rounded to the nearest integer, which is what puts g in the interval abs(g) <= ln2/2 the coefficients were fitted on. The shortcut is gone with it: once g can be nonzero, taking exp(g) as 1 throws away eight digits, and the rational form is exact at g = 0 anyway. The final scaling is done in two halves, because 2^xexp on its own overflows to Inf, and underflows to 0, at exponents where the result is still a finite nonzero number. Worst relative error over the whole range is now about 1.5 eps rather than 223, about 2 ulp.

• EXP returned NaN instead of Inf for very large arguments, and NaN instead of 0 for very large negative ones, from about 2e78 in magnitude and for nearly every argument past about 1e82. The overflow guard sat at 0.9 * .Machine$double.xmax, which is nowhere near where exp overflows, so such an argument reached a range reduction that could not reduce it: g grew with x until the numerator and the denominator of the rational form both overflowed, and their quotient was NaN. The guard is now at log(.Machine$double.xmax) going up and at -1075*ln2 going down, the latter being where exp falls under half the least subnormal. Overflow now warns and returns Inf for every argument past the upper limit, where it used to warn only past 0.9 * .Machine$double.xmax; underflow returns 0 quietly. Both limits agree with exp to the last double.

• EXP(-Inf) returned -Inf. The branch for it had been the +Inf branch with the sign flipped, so it handed the argument back instead of the limit. It returns 0, which is what exp does.

• Given a vector, the functions here that work a value at a time failed with "the condition has length > 1", the message R raises for an if whose condition is not a single value. It names neither the function nor the reason. 130 argument positions across 64 functions, EXP, LOG, SQRT and the distribution functions among them, are now checked up front, and the error names the function and the argument and points at sapply. Arguments that do accept a vector and return the right answer for one, such as mu in Qnorm or x in Dnorm, are untouched. Each of the 64 help pages now names the arguments its own function constrains, so the ones that are scalar throughout and the ones that merely have a scalar control argument, such as standardize in mlr, can be told apart without reading the source.

• The Chol and LDLT examples lost their last line in the rendered help page. A percent sign starts a comment in Rd, so L %*% t(L) was cut back to L. Both are now escaped and the round trip back to M shows.

• Help page titles: invbetai was titled as the inverse of gammp, NORM, DENORM and ellipRange carried their bare names, and seven pairs of pages shared one title (conv and conv0, deconv and deconv0, gammp and gammq, trapez0 and trapez1, Run.test and run.p, Pbinom and Ppois, Qnorm and Qlnorm). Each title now says what its page is about.

• LOG, EXP and SQRT say that current CPUs provide them at the instruction level, so the R versions are there to be read. DESCRIPTION says the same of the package as a whole: it is best studied with its source open.

• Help pages brought back in line with the code. erf and erfc call gammp and gammq, not erfccheb; gammp and gammq call gser or gcf; betai calls betacf and gammln; inverfc goes through the normal quantile; GetP returns solve(H, g); the attribute set by MachEps3 and MachEps4 is "iteration"; the mlr columns are R-Student, COV-Ratio and DFFITS; NORM returns c(f, n); NashSVD takes any real matrix; GAMMA returns Inf at zero as well as at the negative integers; the runifLC recurrence adds k.

• Distribution parameters are named for what they are: mean and standard deviation for the normal and log-normal, location and scale for t, degrees of freedom for chi-square and F, rate for the gamma density, distribution and quantile functions, scale for Rgamma and Rgamma0, mean lambda for Poisson and rate alpha for Rexp.

• Spelling and wording fixed across the help pages, among them Decomposition, Meticulous, Numerical, commutative, array, lengths, chi-square and distribution, together with argument names and function names left over from copied pages.

Version 0.1.3 (2026-08-31)

• Every help page now carries a runnable example. Where a base R equivalent exists the example shows both, so the agreement can be seen at a glance.

• Rbeta was wrong for every case except Beta(1, 1). It built the two gamma deviates with the shape and the scale arguments the other way round, so Rbeta(n, 2, 5) had mean 0.65 where the distribution has 2/7 = 0.286.

• Rgamma divided by bet in its alph == 1 branch where the other two branches multiply, so that one shape came out with mean 1/bet instead of bet.

• Rgamma0 never scaled the accepted deviate by alph, so every shape returned the same Exp(1) times bet: mean 1 rather than 3 for Rgamma0(n, 3, 1). It also now requires alph >= 1, which is the range its rejection step is valid over. Use Rgamma below that.

• simps13 had the Simpson weights 4 and 2 on the wrong points, which left it converging to the wrong value: 21.099 against a true 21.337 even at n = 60. It is now accurate to 1e-12 there.

• Both Simpson rules failed at their smallest legal n, 2 and 3, on a seq() with the wrong sign of step.

• The simps38 example asked for n = 10, which is not a multiple of 3, so the help page showed NULL.

• The deconv0 example divided by a spectrum containing a zero and printed Inf and NaN.

• The erfccheb title no longer claims a Chebyshev approximation, which it stopped using in 0.1.0.

Version 0.1.2 (2026-08-31)

• ChkFx left show.error.message holding a list instead of restoring it. options() returns the previous settings as a list, and that list was being passed back as the value. The restore now passes the list whole and runs from on.exit, so early returns are covered too.

• mlr(Plot = TRUE) no longer opens a device of its own and no longer leaves par() changed. It draws on the current device and restores the previous settings on exit.

• Diagnostics in mlr and deconv go through message() rather than cat(), so they can be suppressed and do not write to stdout.

• The ellipse examples no longer call dev.new().

• Packaging for CRAN: Authors@R replaces the free text author fields, the obsolete LazyLoad field is gone, and the reference manual is no longer shipped under inst/doc, since CRAN builds its own and the copy was 86 percent of the tarball.

Version 0.1.1 (2026-08-31)

• Qpois and Qbinom now follow the quantile convention of base R, the smallest k whose cumulative probability reaches p. They previously returned the largest k whose cumulative probability did not exceed p, one below stats::qpois and stats::qbinom for every p. The two now agree exactly over a wide parameter grid; see tests/.

• Qbinom(p, n, pe) returns n for p = 1, and both functions compare against a tolerance sized to the accuracy of Ppois and Pbinom, about 1e-12. Feeding a cumulative probability straight back in as p can still land on the other side of a tie, because those two are built from the incomplete gamma and beta rather than the saddle point algorithms base R uses.

Version 0.1.0 (2026-08-31)

• First release.

• Calculus, distribution functions, random variate generation, linear algebra, differential equations and optimization, sized for one semester.

• The special functions follow freely licensed reference implementations, each attributed in its own help page, and the distribution and quantile functions are built on them. gammln uses the Lanczos g=7 coefficients of the GNU Scientific Library. gser, gcf, gammp and gammq follow the Cephes series and continued fraction as redistributed in ALGLIB. betacf and betai follow the continued fraction and symmetry transform of the GNU Scientific Library. inverfc uses the Cephes normal quantile ndtri.

• erf and erfc come from the incomplete gamma function through erf(x) = P(1/2, x^2) and erfc(x) = Q(1/2, x^2). That needs no coefficient table of its own and stays accurate far into the tail: erfc holds to 1e-13 relative down to erfc(26) = 5.7e-296, where a rational fit in x alone has long since lost its significant digits.

• invgammp and invbetai are bracketed Newton iterations on the forward functions, following the strategy of gsl_cdf_gamma_Pinv and gsl_cdf_beta_Pinv. Above the median the residual is formed in the upper tail, which avoids the loss of about eight digits that a lower-tail residual suffers for p near 1.

• gammpapprox, betaiapprox and erfccheb are thin wrappers over gammp, betai and erfc.

• Every distribution function is checked against the base R equivalent over a wide parameter grid. See tests/.