Delaporte

While there are breaking changes below, the majority of them would not have been encountered in normal practice. Most users should see no change to results, other than increased speed and fewer failures.

Fixed

• **Breaking:** Internal function log1p now carries a cubic Taylor term. The prior quadratic form left a truncation error of order beta^3 / 3. For large with very small this reached the sixth significant figure of ordinary probabilities. Results for beta > 1e-4 are bitwise-unchanged.
• **Breaking:** Non-finite parameters now return NaN with a warning across d/p/q/rdelap instead of clamp-laundered 0s and 1s.
• **Breaking:** ddelap, pdelap, qdelap, and rdelap with singleton parameters now first attempt to build their PMF and CDF tables using a three-term recurrence relationship derived from the Delaporte distribution's probability generating function. This takes O(K) arithmetic operations instead of O(K^2) log_gamma/exp evaluations. Masses whose leading term underflows---roughly + (1+) > 745---are handled by internal rescaling, allowing accuracy to be retained for very large . This is done natively in linear- or log-space as per the call. This may shift results relative to the old "always sum" algorithm at the level of accumulated rounding: one or two ulps for short tables, growing with table length to roughly the 13th significant digit for very long ones. Both parameters extreme enough to risk overflow in the recurrence fall back to the summation algorithm and vector-valued parameters continue to use the per-element path.
• **Breaking:** ddelap and pdelap now treat arguments below the support the same way base R does. ddelap returns 0 (-Inf with log = TRUE) and pdelap returns a lower-tail CDF of 0 and a survival of 1. Previously these inputs returned NaN with a "NaNs produced" warning.
• **Breaking:** ddelap's internal summand formula subtracted from its running sum before a large canceling pair, lgamma(alpha+i) - lgamma(alpha) (exactly 0 at i = 0), had resolved. When is very small, lgamma(alpha) is itself large, and , when smaller than that intermediate value's rounding unit, was silently discarded before the cancellation completed. This was invisible with log = FALSE, since exp() of the resulting error and exp(0) = 1 are both indistinguishable in linear-space double precision at the affected magnitudes. Parenthesizing the canceling pair so it resolves first fixes both the log-space and linear-space summand formulas. Ordinary (non-degenerate) parameters are bitwise-unchanged.
• **Breaking:** ddelap with log = TRUE now accumulates the PMF in log space using a streaming log-sum-exp instead of taking the log of the linear-space result. This is technically breaking as deep-tail log-probabilities which previously returned -Inf because the linear PMF underflowed below the smallest representable double are now accurate to near full precision. For example, ddelap(2000, 1, 1, 1, log = TRUE) returns approximately -1385.99 instead of -Inf. Hoisting the variables for both log and normal space routines makes them slightly faster than before.
• **Breaking:** pdelap with lower.tail = FALSE now computes the upper tail by direct summation of the PMF when doing so is needed to preserve accuracy, instead of 1 - CDF. Survival probabilities near or below machine epsilon previously returned only floating-point rounding noise. For example, pdelap(60, 1, 1, 1, lower.tail = FALSE) returned roughly 3.3e-16 when the true value is roughly 1.18e-18. The algorithms are now accurate to near full double precision. This is technically breaking, but users should see no change outside of the deep tail which was likely not of interest to most. This summation is used only once the survival probability is small enough that the ordinary complement 1 - CDF would lose meaningful precision: below sqrt(EPS) ~ 1.49e-8. Since the summation's remainder-bound safeguard gives it an O(beta^2) worst-case cost, it is reserved for where it is genuinely needed.
• **Breaking:** pdelap with log.p = TRUE previously computed the linear-space CDF/survival first and took log() of the result, so any value that underflowed to exactly 0 in the deep tail silently returned -Inf even though the true log-probability is finite. log.p = TRUE is now computed directly in log space throughout, matching the log-space guarantee ddelap(..., log = TRUE) already provides, and is accurate to near full double precision even where the linear-space value underflows. pdelap derives log-probabilities in the saturated bands---log-CDF values within about sqrt(.Machine$double.eps) of 0, and their log-survival mirror---from the opposite, relatively accurate side via log1p(-exp(.)) rather than reporting the accumulation's own rounding noise. On the log scale those values now agree with an independent dnbinom* dpois convolution oracle to about 1e-13 relative. In linear space, the original algorithm was already accurate to an ulp, so linear-scale results are unaffected, as is everything outside the bands. All three computational routes---recurrence table, degenerate-ratio elemental fill, and vector-parameter recycling---apply the same repair. A consequence is that qdelap/pdelap round trips with log.p = TRUE are now exact in both tails wherever the input is representable.
• **Breaking:** qdelap for probabilities within accumulated rounding of 1 combined with highly overdispersed parameters previously ground its lookup-table build indefinitely, as the cumulative sum plateaus below such targets no matter how far it extends. These targets are now resolved from the accurate upper-tail survival table, returning the true finite quantile instantly. For example, qdelap(1 - 1e-12, 5, 5, 7000) now returns 7630. Ordinary probabilities are unaffected.
• **Breaking:** qdelap with lower.tail = FALSE and/or log.p = TRUE now searches in the native space of the supplied probabilities instead of transforming them to linear lower-tail space first. The old 1 - p flip collapsed survival probabilities below machine epsilon to exactly 1, returning Inf where qpois resolves a finite quantile, and the old exp() of log-probabilities collapsed deep log-space targets the same way; near-1 targets were inflated by an ulp, stepping round-tripped quantiles up by one. Upper-tail targets are now searched against a survival table built as pdelap, and log-space targets against log-probability tables accumulated by streaming log-sum-exp. The default (lower.tail = TRUE, log.p = FALSE) path is bitwise-unchanged. Interior search targets carry one-sided tolerances sized to their measured noise sources, in the spirit of base R's 64 * EPS fuzz in qpois.c. The vector-parameter path retains base R's fuzz verbatim.