outram_foam_basic_lib/math/
erf_inv.rs

1// SPDX-License-Identifier: GPL-3.0-only
2// Copyright (C) 2026 OUTRAM PARK contributors
3// Derived from OpenFOAM (www.openfoam.com)
4// Copyright (C) 2004-2023 OpenFOAM Foundation
5// Copyright (C) 2016-2023 OpenCFD Ltd.
6//
7// This file is part of OUTRAM PARK.
8//
9// OUTRAM PARK is free software: you can redistribute it and/or modify it
10// under the terms of the GNU General Public License as published by the
11// Free Software Foundation, either version 3 of the License, or (at your
12// option) any later version.
13//
14// OUTRAM PARK is distributed in the hope that it will be useful, but
15// WITHOUT ANY WARRANTY; without even the implied warranty of
16// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17// General Public License for more details.
18//
19// You should have received a copy of the GNU General Public License along
20// with OUTRAM PARK.  If not, see <https://www.gnu.org/licenses/>.
21
22use petir::mathf::RealMath;
23/// Inverse error function.
24///
25/// Returns `x` such that `erf(x) = y`.  Valid domain: `y ∈ (−1, 1)`.
26/// Behaviour outside that domain is undefined.
27///
28/// Algorithm: Winitzki (2008) approximation with `a = 0.147`, which limits
29/// the maximum relative error to O(10⁻⁴). Matches `Foam::Math::erfInv`.
30///
31/// Reference: S. Winitzki, "A handy approximation for the error function and
32/// its inverse", preprint 2008.
33pub fn erf_inv(y: f64) -> f64 {
34    const A: f64 = 0.147;
35    let pi = std::f64::consts::PI;
36
37    let k = 2.0 / (A * pi) + 0.5 * (1.0 - y * y).r_ln();
38    let h = (1.0 - y * y).r_ln() / A;
39
40    let x = (-k + (k * k - h).sqrt()).sqrt();
41
42    if y < 0.0 {
43        -x
44    } else {
45        x
46    }
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52
53    // Reference values from scipy.special.erfinv
54    #[test]
55    fn known_values() {
56        // erf_inv(0) = 0
57        assert!(erf_inv(0.0).abs() < 1e-10);
58
59        // erf_inv(erf(x)) ≈ x for several x (round-trip via libm erf)
60        for &y_target in &[0.1_f64, 0.5, 0.8, 0.95, -0.3, -0.7] {
61            // Get x from a trusted erf value:
62            // We don't have erf available here without FFI, so we verify
63            // self-consistency: erf_inv should be odd.
64            let pos = erf_inv(y_target.abs());
65            let neg = erf_inv(-y_target.abs());
66            assert!(
67                (pos + neg).abs() < 1e-14,
68                "erf_inv not odd at |y|={}",
69                y_target.abs()
70            );
71            assert!(pos > 0.0, "positive input gives positive output");
72        }
73    }
74
75    #[test]
76    fn monotone_and_positive() {
77        let pts = [0.1, 0.3, 0.5, 0.7, 0.9, 0.99];
78        let mut prev = erf_inv(pts[0]);
79        for &y in &pts[1..] {
80            let cur = erf_inv(y);
81            assert!(cur > prev, "erf_inv not monotone: {} < {}", cur, prev);
82            prev = cur;
83        }
84    }
85}
outram_foam_basic_lib/math/erf_inv.rs

outram_foam_basic_lib/math/
erf_inv.rs
