Fix the accuracy of ERF and ERFC near zero and for large arguments

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)

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 {