diff --git a/CHANGELOG.md b/CHANGELOG.md index 787c642ad9..c21bb91258 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,6 +19,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), ### Fixed - Fixed the validation of classic (25-character) license keys depending on the time zone: east of UTC, a key that expired the day before the build was released was still accepted, and west of UTC, the console message printed an expiry date one day too early. [#1728](https://github.com/handsontable/hyperformula/pull/1728) +- Fixed the `ERF` and `ERFC` functions losing precision for arguments close to zero and for large arguments, for example `ERF(1E-10)` and `ERFC(5)`; the results now agree with Excel. [#1795](https://github.com/handsontable/hyperformula/pull/1795) - Fixed the `AVERAGEIF` function returning a division-by-zero error when the calculated average was `0`. [#1733](https://github.com/handsontable/hyperformula/pull/1733) - Fixed the localized names of `VSTACK` and `HSTACK` in 14 language packs to match Microsoft Excel. [#1748](https://github.com/handsontable/hyperformula/pull/1748) - Fixed the MAXPOOL and MEDIANPOOL functions throwing an uncaught `TypeError` instead of returning the `#VALUE!` error when the range dimensions are not a whole multiple of the window size and the stride. [#1718](https://github.com/handsontable/hyperformula/pull/1718) diff --git a/src/interpreter/plugin/3rdparty/jstat/jstat.ts b/src/interpreter/plugin/3rdparty/jstat/jstat.ts index 94e10c8dd5..caff10708f 100644 --- a/src/interpreter/plugin/3rdparty/jstat/jstat.ts +++ b/src/interpreter/plugin/3rdparty/jstat/jstat.ts @@ -21,46 +21,67 @@ THE SOFTWARE. */ -export function erf(x: number): number { - const cof = [-1.3026537197817094, 6.4196979235649026e-1, 1.9476473204185836e-2, - -9.561514786808631e-3, -9.46595344482036e-4, 3.66839497852761e-4, - 4.2523324806907e-5, -2.0278578112534e-5, -1.624290004647e-6, - 1.303655835580e-6, 1.5626441722e-8, -8.5238095915e-8, - 6.529054439e-9, 5.059343495e-9, -9.91364156e-10, - -2.27365122e-10, 9.6467911e-11, 2.394038e-12, - -6.886027e-12, 8.94487e-13, 3.13092e-13, - -1.12708e-13, 3.81e-16, 7.106e-15, - -1.523e-15, -9.4e-17, 1.21e-16, - -2.8e-17] - let j = cof.length - 1 - let isneg = false - let d = 0 - let dd = 0 - let t, ty, tmp, res - - if (x === 0) { - return 0 - } - if (x < 0) { - x = -x - isneg = true - } +const SQRT_PI = Math.sqrt(Math.PI) - t = 2 / (2 + x) - ty = 4 * t - 2 +/** + * erf(x) for 0 <= x < 2 from its series exp(-x^2) * 2 / sqrt(pi) * sum(2^n * x^(2n + 1) / (2n + 1)!!). + * All terms are positive, so there is no cancellation, and the relative error stays near 1e-16 down to the + * smallest numbers (the earlier 1 - exp(...) form lost every digit below about 1e-16). + */ +function erfSeries(x: number): number { + const xx = x * x + let term = x + let sum = x + for (let n = 0; n < 500; n++) { + term *= 2 * xx / (2 * n + 3) + sum += term + if (term < sum * 1e-17) { + break + } + } + return 2 / SQRT_PI * Math.exp(-xx) * sum +} - for (; j > 0; j--) { - tmp = d - d = ty * d - dd + cof[j] - dd = tmp +/** + * erfc(x) for x >= 1 from the continued fraction exp(-x^2) / sqrt(pi) / (x + (1/2) / (x + 1 / (x + (3/2) / (x + ...)))), + * evaluated from the tail. Accurate to about 1e-16 relative, also where erfc is tiny and 1 - erf(x) would be 0. + */ +function erfcFraction(x: number): number { + let tail = 0 + for (let k = 300; k >= 1; k--) { + tail = (k / 2) / (x + tail) } + return Math.exp(-x * x) / SQRT_PI / (x + tail) +} - res = t * Math.exp(-x * x + 0.5 * (cof[0] + ty * d) - dd) - return isneg ? res - 1 : 1 - res +export function erf(x: number): number { + if (x === 0 || Number.isNaN(x)) { + return x + } + const abs = Math.abs(x) + let result + if (abs < 2) { + result = erfSeries(abs) + } else if (abs < 6.5) { + result = 1 - erfcFraction(abs) + } else { + result = 1 + } + return x < 0 ? -result : result } export function erfc(x: number): number { - return 1 - erf(x) + if (Number.isNaN(x)) { + return x + } + if (x < 0) { + return 2 - erfc(-x) + } + if (x < 1) { + return 1 - erfSeries(x) + } + // erfc(27.3) is below the smallest double + return x < 27.3 ? erfcFraction(x) : 0 } function erfcinv(p: number): number {