diff --git a/src/__tests__/__snapshots__/index.test.ts.snap b/src/__tests__/__snapshots__/index.test.ts.snap index 7b0f8eb18..3efb5b356 100644 --- a/src/__tests__/__snapshots__/index.test.ts.snap +++ b/src/__tests__/__snapshots__/index.test.ts.snap @@ -192,6 +192,7 @@ exports[`existence of exported functions 1`] = ` "createFromToArray", "createRandomArray", "createStepArray", + "erfc", "getCombinations", "getCombinationsIterator", "getRescaler", diff --git a/src/utils/__tests__/erfc.test.ts b/src/utils/__tests__/erfc.test.ts new file mode 100644 index 000000000..4222dba7d --- /dev/null +++ b/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)); +}); diff --git a/src/utils/erfc.ts b/src/utils/erfc.ts new file mode 100644 index 000000000..e59029260 --- /dev/null +++ b/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; +} diff --git a/src/utils/index.ts b/src/utils/index.ts index 58b8388cb..9935ba5c6 100644 --- a/src/utils/index.ts +++ b/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';