pxfm/err/
erfcxf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 9/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::f_fmla;
30use crate::exponents::core_expdf;
31use crate::polyeval::{f_estrin_polyeval8, f_polyeval6};
32
33#[inline]
34fn core_erfcx(x: f32) -> f64 {
35    // x here is already always > 1
36    let dx = x as f64;
37    if x < 8. {
38        // Rational approximant generated by Wolfram Mathematica:
39        // <<FunctionApproximations`
40        // ClearAll["Global`*"]
41        // f[x_]:=Exp[x^2]Erfc[x]
42        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{1,8},7,7},WorkingPrecision->75,MaxIterations->100]
43        // num=Numerator[approx];
44        // den=Denominator[approx];
45        // coeffs=CoefficientList[num,z];
46        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
47        // coeffs=CoefficientList[den,z];
48        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
49        let p_num = f_estrin_polyeval8(
50            dx,
51            f64::from_bits(0x3ff00000804c8f8f),
52            f64::from_bits(0x3ffb7307ea8fdbeb),
53            f64::from_bits(0x3ff7081ba7bc735c),
54            f64::from_bits(0x3fe767338b33532a),
55            f64::from_bits(0x3fce3c8288507fd6),
56            f64::from_bits(0x3fa7ca2cb4ae697f),
57            f64::from_bits(0x3f72b11b0dfb2348),
58            f64::from_bits(0xbd9f64f0c15c479b),
59        );
60        let p_den = f_estrin_polyeval8(
61            dx,
62            f64::from_bits(0x3ff0000000000000),
63            f64::from_bits(0x4006c071e850132e),
64            f64::from_bits(0x400d30326bc347ee),
65            f64::from_bits(0x40060d8d56bada75),
66            f64::from_bits(0x3ff56643fc4580eb),
67            f64::from_bits(0x3fdb0e194e72a513),
68            f64::from_bits(0x3fb5154759b61be3),
69            f64::from_bits(0x3f8090b063cce524),
70        );
71        return p_num / p_den;
72    }
73    // for large x erfcx(x) ~ 1/sqrt(pi) / x * R(1/x)
74    const ONE_OVER_SQRT_PI: f64 = f64::from_bits(0x3fe20dd750429b6d);
75    let r = 1. / dx;
76    // Rational approximant generated by Wolfram Mathematica:
77    // <<FunctionApproximations`
78    // ClearAll["Global`*"]
79    // f[x_]:=Exp[1/x^2]Erfc[1/x]/x*Sqrt[Pi]
80    // {err0,approx}=MiniMaxApproximation[f[z],{z,{2^-12,1/8},5,5},WorkingPrecision->75,MaxIterations->100]
81    // num=Numerator[approx][[1]];
82    // den=Denominator[approx][[1]];
83    // coeffs=CoefficientList[num,z];
84    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
85    // coeffs=CoefficientList[den,z];
86    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
87    let p_num = f_polyeval6(
88        r,
89        f64::from_bits(0x3ff0000000000002),
90        f64::from_bits(0xbfd09caf2bb541c3),
91        f64::from_bits(0x40132238367ae454),
92        f64::from_bits(0xc0060bc62c3711b1),
93        f64::from_bits(0x40024a90d229158d),
94        f64::from_bits(0xc0013665d8ff3813),
95    );
96    let p_den = f_polyeval6(
97        r,
98        f64::from_bits(0x3ff0000000000000),
99        f64::from_bits(0xbfd09caf2bb5101d),
100        f64::from_bits(0x4015223836772f2c),
101        f64::from_bits(0xc00715911b5f5f5c),
102        f64::from_bits(0x4010b66411ec4e1f),
103        f64::from_bits(0xc00b325c767ed436),
104    );
105    (r * ONE_OVER_SQRT_PI) * (p_num / p_den)
106}
107
108/// Scaled complementary error function (exp(x^2)*erfc(x))
109///
110/// ulp 0.5
111pub fn f_erfcxf(x: f32) -> f32 {
112    let ux = x.to_bits().wrapping_shl(1);
113    if ux >= 0xffu32 << 24 || ux <= 0x6499_999au32 {
114        // |x| == 0, |x| == inf, |x| == NaN, |x| <= 1.19209290e-08
115        if ux <= 0x6499_999au32 {
116            // |x| == 0, |x| <= 1.19209290e-08
117            return 1.;
118        }
119        if x.is_infinite() {
120            return if x.is_sign_positive() {
121                0.
122            } else {
123                f32::INFINITY
124            };
125        }
126        return f32::NAN; // x == NaN
127    }
128    let ax = x.to_bits() & 0x7fff_ffff;
129    if x <= -9.382415 {
130        // x <= -9.382415
131        return f32::INFINITY;
132    }
133    if ax <= 0x34000000u32 {
134        // |x| < ulp(1) we use taylor series at 0
135        // erfcx(x) ~ 1-(2 x)/Sqrt[\[Pi]]+x^2-(4 x^3)/(3 Sqrt[\[Pi]])+x^4/2-(8 x^5)/(15 Sqrt[\[Pi]])+O[x]^6
136        #[cfg(any(
137            all(
138                any(target_arch = "x86", target_arch = "x86_64"),
139                target_feature = "fma"
140            ),
141            all(target_arch = "aarch64", target_feature = "neon")
142        ))]
143        {
144            use crate::common::f_fmlaf;
145            const M_TWO_OVER_SQRT_PI: f32 = f32::from_bits(0xbf906ebb);
146            return f_fmlaf(x, M_TWO_OVER_SQRT_PI, 1.);
147        }
148        #[cfg(not(any(
149            all(
150                any(target_arch = "x86", target_arch = "x86_64"),
151                target_feature = "fma"
152            ),
153            all(target_arch = "aarch64", target_feature = "neon")
154        )))]
155        {
156            use crate::common::f_fmla;
157            const M_TWO_OVER_SQRT_PI: f64 = f64::from_bits(0xbff20dd750429b6d);
158            let dx = x as f64;
159            return f_fmla(dx, M_TWO_OVER_SQRT_PI, 1.) as f32;
160        }
161    }
162
163    if ax <= 0x3f800000u32 {
164        // |x| <= 1
165        let dx = x as f64;
166        // Generated by Wolfram Mathematica:
167        // <<FunctionApproximations`
168        // ClearAll["Global`*"]
169        // f[x_]:=Exp[x^2]Erfc[x]
170        // {err0,approx}=MiniMaxApproximation[f[z],{z,{-1,1},7,7},WorkingPrecision->75,MaxIterations->100]
171        // num=Numerator[approx][[1]];
172        // den=Denominator[approx][[1]];
173        // coeffs=CoefficientList[num,z];
174        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
175        // coeffs=CoefficientList[den,z];
176        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
177        let p_num = f_estrin_polyeval8(
178            dx,
179            f64::from_bits(0x3feffffffffffff8),
180            f64::from_bits(0x3ff26c328bd2dc5f),
181            f64::from_bits(0x3fe6f91b9fa5f58c),
182            f64::from_bits(0x3fd09edf3fcf5ee1),
183            f64::from_bits(0x3faddb3bcedbff91),
184            f64::from_bits(0x3f7e43b5dd4b7587),
185            f64::from_bits(0x3f3baab6b3e61d7b),
186            f64::from_bits(0xbe83e7d629825321),
187        );
188        let p_den = f_estrin_polyeval8(
189            dx,
190            f64::from_bits(0x3ff0000000000000),
191            f64::from_bits(0x40023d04ee0abc28),
192            f64::from_bits(0x400252b377263d61),
193            f64::from_bits(0x3ff510af7f826479),
194            f64::from_bits(0x3fddfc089c4731ed),
195            f64::from_bits(0x3fba79b040e28b0a),
196            f64::from_bits(0x3f8aea2f3579235a),
197            f64::from_bits(0x3f485d2875b4f88c),
198        );
199        return (p_num / p_den) as f32;
200    }
201
202    let erfcx_abs_x = core_erfcx(f32::from_bits(ax));
203    if x < 0. {
204        // exp(x^2)erfc(-x) = 2*exp(x^2) - erfcx(|x|)
205        let dx = x as f64;
206        return f_fmla(2., core_expdf(dx * dx), -erfcx_abs_x) as f32;
207    }
208    erfcx_abs_x as f32
209}
210
211#[cfg(test)]
212mod tests {
213    use super::*;
214    #[test]
215    fn test_erfcx() {
216        assert_eq!(f_erfcxf(5.19209290e-09), 1.0);
217        assert_eq!(f_erfcxf(1.19209290e-08), 1.0);
218        assert_eq!(f_erfcxf(f32::EPSILON), 0.9999999);
219        assert_eq!(f_erfcxf(12.1), 0.046469606);
220        assert_eq!(f_erfcxf(7.1), 0.07869752);
221        assert_eq!(f_erfcxf(1.1), 0.40173045);
222        assert_eq!(f_erfcxf(-0.23), 1.3232007);
223        assert_eq!(f_erfcxf(-1.4325), 15.234794);
224        assert_eq!(f_erfcxf(-10.), f32::INFINITY);
225        assert_eq!(f_erfcxf(f32::INFINITY), 0.);
226        assert_eq!(f_erfcxf(f32::NEG_INFINITY), f32::INFINITY);
227        assert!(f_erfcxf(f32::NAN).is_nan());
228    }
229}
pxfm/err/erfcxf.rs

pxfm/err/
erfcxf.rs
