feat: implement complementary error function (erfc) and add tests

src/__tests__/__snapshots__/index.test.ts.snap

@@ -192,6 +192,7 @@ exports[`existence of exported functions 1`] = `
   "createFromToArray",
   "createRandomArray",
   "createStepArray",
+  "erfc",
   "getCombinations",
   "getCombinationsIterator",
   "getRescaler",

src/utils/__tests__/erfc.test.ts

@@ -0,0 +1,19 @@
+import { expect, test } from 'vitest';
+
+import { erfc } from '../erfc.ts';
+
+test('should match the exact special values', () => {
+  expect(erfc(0)).toBeCloseTo(1, 12);
+  expect(erfc(Number.POSITIVE_INFINITY)).toBeCloseTo(0, 12);
+  expect(erfc(Number.NEGATIVE_INFINITY)).toBeCloseTo(2, 12);
+});
+
+test('should respect the symmetry relation erfc(-x) = 2 - erfc(x)', () => {
+  expect(erfc(-1)).toBeCloseTo(2 - erfc(1), 10);
+  expect(erfc(-2)).toBeCloseTo(2 - erfc(2), 10);
+});
+
+test('should be monotone decreasing on the real line', () => {
+  expect(erfc(0.5)).toBeGreaterThan(erfc(1));
+  expect(erfc(1)).toBeGreaterThan(erfc(2));
+});

src/utils/erfc.ts

@@ -0,0 +1,49 @@
+/**
+ * Evaluates the complementary error function erfc(x).
+ *
+ * This implementation uses the Abramowitz-Stegun approximation:
+ *
+ * erfc(x) ≈ (a1 t + a2 t^2 + ... + a5 t^5) * exp(-x^2)
+ *
+ * with t = 1 / (1 + p * |x|), and an exact symmetry relation for negative values.
+ *
+ * The approximation is accurate for typical numeric work and is suitable for use
+ * in a Newton iteration for the inverse complementary error function.
+ *
+ * @param x - input value.
+ * @returns the complementary error function value for x.
+ */
+export function erfc(x: number): number {
+  if (Number.isNaN(x)) {
+    return Number.NaN;
+  }
+
+  if (x === 0) {
+    return 1;
+  }
+
+  if (x === Number.POSITIVE_INFINITY) {
+    return 0;
+  }
+
+  if (x === Number.NEGATIVE_INFINITY) {
+    return 2;
+  }
+
+  const sign = x < 0 ? -1 : 1;
+  const ax = Math.abs(x);
+
+  const p = 0.3275911;
+
+  const a1 = 0.254829592;
+  const a2 = -0.284496736;
+  const a3 = 1.421413741;
+  const a4 = -1.453152027;
+  const a5 = 1.061405429;
+
+  const t = 1 / (1 + p * ax);
+  const polynomial = ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t;
+  const result = polynomial * Math.exp(-ax * ax);
+
+  return sign >= 0 ? result : 2 - result;
+}

src/utils/index.ts

@@ -2,6 +2,7 @@ export * from './createArray.ts';
 export * from './createFromToArray.ts';
 export * from './createRandomArray.ts';
 export * from './createStepArray.ts';
+export * from './erfc.ts';
 export { clearFFTCache, setFFTCacheMaxSize } from './fftCache.ts';
 export * from './getCombinations.ts';
 export * from './getCombinationsIterator.ts';