pxfm/err/
erfcx.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::double_double::DoubleDouble;
30use crate::pow_exec::exp_dd_fast;
31
32#[inline]
33fn core_erfcx(x: f64) -> DoubleDouble {
34    if x <= 8. {
35        // Rational approximant for erfcx generated by Wolfram Mathematica:
36        // <<FunctionApproximations`
37        // ClearAll["Global`*"]
38        // f[x_]:=Exp[x^2]Erfc[x]
39        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{1, 8},11,11},WorkingPrecision->75,MaxIterations->100]
40        // num=Numerator[approx];
41        // den=Denominator[approx];
42        // coeffs=CoefficientList[num,z];
43        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
44        // coeffs=CoefficientList[den,z];
45        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
46        const P: [(u64, u64); 12] = [
47            (0xbc836faeb9a312bb, 0x3ff000000000ee8e),
48            (0x3c91842f891bec6a, 0x4002ca20a78aaf8f),
49            (0x3c7916e8a1c30681, 0x4005e955f70aed5b),
50            (0x3cabad150d828d82, 0x4000646f5807ad07),
51            (0xbc6f482680d66d9c, 0x3ff1449e03ed381c),
52            (0xbc7188796156ae19, 0x3fdaa7e997e3b034),
53            (0xbc5c8af0642761e3, 0x3fbe836282058d4a),
54            (0xbc372829be2d072f, 0x3f99a2b2adc2ec05),
55            (0x3c020cc8b96000ab, 0x3f6e6cc3d120a955),
56            (0x3bdd138e6c136806, 0x3f3743d6735eaf13),
57            (0xbb9fbd22f0675122, 0x3ef1c1d36ebe29a2),
58            (0xb89093cc981c934c, 0xbc43c18bc6385c74),
59        ];
60
61        let x2 = DoubleDouble::from_exact_mult(x, x);
62        let x4 = x2 * x2;
63        let x8 = x4 * x4;
64
65        let e0 = DoubleDouble::mul_f64_add(
66            DoubleDouble::from_bit_pair(P[1]),
67            x,
68            DoubleDouble::from_bit_pair(P[0]),
69        );
70        let e1 = DoubleDouble::mul_f64_add(
71            DoubleDouble::from_bit_pair(P[3]),
72            x,
73            DoubleDouble::from_bit_pair(P[2]),
74        );
75        let e2 = DoubleDouble::mul_f64_add(
76            DoubleDouble::from_bit_pair(P[5]),
77            x,
78            DoubleDouble::from_bit_pair(P[4]),
79        );
80        let e3 = DoubleDouble::mul_f64_add(
81            DoubleDouble::from_bit_pair(P[7]),
82            x,
83            DoubleDouble::from_bit_pair(P[6]),
84        );
85        let e4 = DoubleDouble::mul_f64_add(
86            DoubleDouble::from_bit_pair(P[9]),
87            x,
88            DoubleDouble::from_bit_pair(P[8]),
89        );
90        let e5 = DoubleDouble::mul_f64_add(
91            DoubleDouble::from_bit_pair(P[11]),
92            x,
93            DoubleDouble::from_bit_pair(P[10]),
94        );
95
96        let f0 = DoubleDouble::mul_add(x2, e1, e0);
97        let f1 = DoubleDouble::mul_add(x2, e3, e2);
98        let f2 = DoubleDouble::mul_add(x2, e5, e4);
99
100        let g0 = DoubleDouble::mul_add(x4, f1, f0);
101
102        let p_num = DoubleDouble::mul_add(x8, f2, g0);
103
104        const Q: [(u64, u64); 12] = [
105            (0x0000000000000000, 0x3ff0000000000000),
106            (0xbc95d65be031374e, 0x400bd10c4fb1dbe5),
107            (0x3cb2d8f661db08a0, 0x4016a649ff973199),
108            (0x3ca32cbcfdc0ea93, 0x4016daab399c1ffc),
109            (0xbca2982868536578, 0x400fd61ab892d14c),
110            (0xbca2e29199e17fd9, 0x40001f56c4d495a3),
111            (0x3c412ce623a1790a, 0x3fe852b582135164),
112            (0x3c61152eaf4b0dc5, 0x3fcb760564da7cde),
113            (0xbc1b57ff91d81959, 0x3fa6e146988df835),
114            (0x3c17183d8445f19a, 0x3f7b06599b5e912f),
115            (0xbbd0ada61b85ff98, 0x3f449e39467b73d0),
116            (0xbb658d84fc735e67, 0x3eff794442532b51),
117        ];
118
119        let e0 = DoubleDouble::mul_f64_add_f64(
120            DoubleDouble::from_bit_pair(Q[1]),
121            x,
122            f64::from_bits(0x3ff0000000000000),
123        );
124        let e1 = DoubleDouble::mul_f64_add(
125            DoubleDouble::from_bit_pair(Q[3]),
126            x,
127            DoubleDouble::from_bit_pair(Q[2]),
128        );
129        let e2 = DoubleDouble::mul_f64_add(
130            DoubleDouble::from_bit_pair(Q[5]),
131            x,
132            DoubleDouble::from_bit_pair(Q[4]),
133        );
134        let e3 = DoubleDouble::mul_f64_add(
135            DoubleDouble::from_bit_pair(Q[7]),
136            x,
137            DoubleDouble::from_bit_pair(Q[6]),
138        );
139        let e4 = DoubleDouble::mul_f64_add(
140            DoubleDouble::from_bit_pair(Q[9]),
141            x,
142            DoubleDouble::from_bit_pair(Q[8]),
143        );
144        let e5 = DoubleDouble::mul_f64_add(
145            DoubleDouble::from_bit_pair(Q[11]),
146            x,
147            DoubleDouble::from_bit_pair(Q[10]),
148        );
149
150        let f0 = DoubleDouble::mul_add(x2, e1, e0);
151        let f1 = DoubleDouble::mul_add(x2, e3, e2);
152        let f2 = DoubleDouble::mul_add(x2, e5, e4);
153
154        let g0 = DoubleDouble::mul_add(x4, f1, f0);
155
156        let p_den = DoubleDouble::mul_add(x8, f2, g0);
157        return DoubleDouble::div(p_num, p_den);
158    }
159
160    // for large x erfcx(x) ~ 1/sqrt(pi) / x * R(1/x)
161    const ONE_OVER_SQRT_PI: DoubleDouble =
162        DoubleDouble::from_bit_pair((0x3c61ae3a914fed80, 0x3fe20dd750429b6d));
163    let r = DoubleDouble::from_quick_recip(x);
164    // Rational approximant generated by Wolfram:
165    // <<FunctionApproximations`
166    // ClearAll["Global`*"]
167    // f[x_]:=Exp[1/x^2]Erfc[1/x]/x*Sqrt[Pi]
168    // {err0,approx}=MiniMaxApproximation[f[z],{z,{2^-23,1/8},8,8},WorkingPrecision->75,MaxIterations->100]
169    // num=Numerator[approx][[1]];
170    // den=Denominator[approx][[1]];
171    // coeffs=CoefficientList[num,z];
172    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
173    // coeffs=CoefficientList[den,z];
174    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
175    const P: [(u64, u64); 9] = [
176        (0xbb1d2ee37e46a4cd, 0x3ff0000000000000),
177        (0x3ca2e575a4ce3d30, 0x4001303ab00c8bac),
178        (0xbccf38381e5ee521, 0x4030a97aeed54c9f),
179        (0xbcc3a2842df0dd3d, 0x4036f7733c9fd2f9),
180        (0xbcfeaf46506f16ed, 0x4051c5f382750553),
181        (0x3ccbb9f5e11d176a, 0x404ac0081e0749e0),
182        (0xbcf374f8966ae2a5, 0x4052082526d99a5c),
183        (0x3cbb5530b924f224, 0x402feabbf6571c29),
184        (0xbcbcdd50a3ca4776, 0x40118726e1f2d204),
185    ];
186    const Q: [(u64, u64); 9] = [
187        (0x0000000000000000, 0x3ff0000000000000),
188        (0x3ca2e4613c9e0017, 0x4001303ab00c8bac),
189        (0xbcce5f17cf14e51d, 0x4031297aeed54c9f),
190        (0xbcdf7e0fed176f92, 0x40380a76e7a09bb2),
191        (0x3cfc57b67a2797af, 0x4053bb22e04faf3e),
192        (0xbcd3e63b7410b46b, 0x404ff46317ae9483),
193        (0xbce122c15db2653f, 0x405925ef8a428c36),
194        (0x3ce174ebe3e52c8e, 0x4040f49acfe692e1),
195        (0xbcda0e267ce6e2e6, 0x40351a07878bfbd3),
196    ];
197    let mut p_num = DoubleDouble::mul_add(
198        DoubleDouble::from_bit_pair(P[8]),
199        r,
200        DoubleDouble::from_bit_pair(P[7]),
201    );
202    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[6]));
203    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[5]));
204    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[4]));
205    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[3]));
206    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[2]));
207    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[1]));
208    p_num = DoubleDouble::mul_add(p_num, r, DoubleDouble::from_bit_pair(P[0]));
209
210    let mut p_den = DoubleDouble::mul_add(
211        DoubleDouble::from_bit_pair(Q[8]),
212        r,
213        DoubleDouble::from_bit_pair(Q[7]),
214    );
215    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[6]));
216    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[5]));
217    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[4]));
218    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[3]));
219    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[2]));
220    p_den = DoubleDouble::mul_add(p_den, r, DoubleDouble::from_bit_pair(Q[1]));
221    p_den = DoubleDouble::mul_add_f64(p_den, r, f64::from_bits(0x3ff0000000000000));
222
223    let v0 = DoubleDouble::quick_mult(ONE_OVER_SQRT_PI, r);
224    let v1 = DoubleDouble::div(p_num, p_den);
225    DoubleDouble::quick_mult(v0, v1)
226}
227
228/// Scaled complementary error function (exp(x^2)*erfc(x))
229pub fn f_erfcx(x: f64) -> f64 {
230    let ux = x.to_bits().wrapping_shl(1);
231
232    if ux >= 0x7ffu64 << 53 || ux <= 0x7960000000000000u64 {
233        // x == NaN, x == inf, x == 0, |x| <= f64::EPSILON
234        if x.is_nan() {
235            return f64::NAN;
236        }
237        if x.to_bits().wrapping_shl(1) == 0 {
238            return 1.;
239        }
240        if x.is_infinite() {
241            return if x.is_sign_positive() {
242                0.
243            } else {
244                f64::INFINITY
245            };
246        }
247
248        if ux <= 0x7888f5c28f5c28f6u64 {
249            // |x| <= 2.2204460492503131e-18
250            return 1.;
251        }
252
253        // |x| <= f64::EPSILON
254        use crate::common::f_fmla;
255        const M_TWO_OVER_SQRT_PI: DoubleDouble =
256            DoubleDouble::from_bit_pair((0xbc71ae3a914fed80, 0xbff20dd750429b6d));
257        return f_fmla(
258            M_TWO_OVER_SQRT_PI.lo,
259            x,
260            f_fmla(M_TWO_OVER_SQRT_PI.hi, x, 1.),
261        );
262    }
263
264    if x.to_bits() >= 0xc03aa449ebc84dd6 {
265        // x <= -sqrt(709.783) ~ -26.6417
266        return f64::INFINITY;
267    }
268
269    let ax = x.to_bits() & 0x7fff_ffff_ffff_ffffu64;
270
271    if ax <= 0x3ff0000000000000u64 {
272        // |x| <= 1
273        // Rational approximant generated by Wolfram Mathematica:
274        // <<FunctionApproximations`
275        // ClearAll["Global`*"]
276        // f[x_]:=Exp[x^2]Erfc[x]
277        // {err0,approx}=MiniMaxApproximation[f[z],{z,{-1, 1},10,10},WorkingPrecision->75,MaxIterations->100]
278        // num=Numerator[approx][[1]];
279        // den=Denominator[approx][[1]];
280        // coeffs=CoefficientList[num,z];
281        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
282        // coeffs=CoefficientList[den,z];
283        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
284        const P: [(u64, u64); 11] = [
285            (0xbb488611350b1950, 0x3ff0000000000000),
286            (0xbc86ae482c7f2342, 0x3ff9c5d39e89602f),
287            (0x3c6702d70b807254, 0x3ff5a4c406d6468b),
288            (0x3c7fe41fc43cfed5, 0x3fe708e7f401bd0c),
289            (0x3c73a4a355172c6d, 0x3fd0d9a0c1a7126c),
290            (0x3c5f4c372faa270f, 0x3fb154722e30762e),
291            (0xbc04c0227976379e, 0x3f88ecebb62ce646),
292            (0xbbdc9ea151b9eb33, 0x3f580ea84143877b),
293            (0xbb6dae7001a91491, 0x3f1c3c5f95579b0a),
294            (0x3b6aca5e82c52897, 0x3ecea4db51968d9e),
295            (0x3a41c4edd175d2af, 0x3dbc0dccea7fc8ed),
296        ];
297
298        let x2 = DoubleDouble::from_exact_mult(x, x);
299        let x4 = x2 * x2;
300        let x8 = x4 * x4;
301
302        let q0 = DoubleDouble::mul_f64_add(
303            DoubleDouble::from_bit_pair(P[1]),
304            x,
305            DoubleDouble::from_bit_pair(P[0]),
306        );
307        let q1 = DoubleDouble::mul_f64_add(
308            DoubleDouble::from_bit_pair(P[3]),
309            x,
310            DoubleDouble::from_bit_pair(P[2]),
311        );
312        let q2 = DoubleDouble::mul_f64_add(
313            DoubleDouble::from_bit_pair(P[5]),
314            x,
315            DoubleDouble::from_bit_pair(P[4]),
316        );
317        let q3 = DoubleDouble::mul_f64_add(
318            DoubleDouble::from_bit_pair(P[7]),
319            x,
320            DoubleDouble::from_bit_pair(P[6]),
321        );
322        let q4 = DoubleDouble::mul_f64_add(
323            DoubleDouble::from_bit_pair(P[9]),
324            x,
325            DoubleDouble::from_bit_pair(P[8]),
326        );
327
328        let r0 = DoubleDouble::mul_add(x2, q1, q0);
329        let r1 = DoubleDouble::mul_add(x2, q3, q2);
330
331        let s0 = DoubleDouble::mul_add(x4, r1, r0);
332        let s1 = DoubleDouble::mul_add(x2, DoubleDouble::from_bit_pair(P[10]), q4);
333        let p_num = DoubleDouble::mul_add(x8, s1, s0);
334
335        const Q: [(u64, u64); 11] = [
336            (0x0000000000000000, 0x3ff0000000000000),
337            (0xbc7bae414cad99c8, 0x4005e9d57765fdce),
338            (0x3c8fa553bed15758, 0x400b8c670b3fbcda),
339            (0x3ca6c7ad610f1019, 0x4004f2ca59958153),
340            (0x3c87787f336cc4e6, 0x3ff55c267090315a),
341            (0xbc6ef55d4b2c4150, 0x3fde8b84b64b6f4e),
342            (0x3c570d63c94be3a3, 0x3fbf0d5e36017482),
343            (0x3c1882a745ef572e, 0x3f962f73633506c1),
344            (0xbc0850bb6fc82764, 0x3f65593e0dc46acd),
345            (0xbbb9dc0097d7d776, 0x3f290545603e2f94),
346            (0xbb776e5781e3889d, 0x3edb29c49d18cf89),
347        ];
348
349        let q0 = DoubleDouble::mul_f64_add_f64(
350            DoubleDouble::from_bit_pair(Q[1]),
351            x,
352            f64::from_bits(0x3ff0000000000000),
353        );
354        let q1 = DoubleDouble::mul_f64_add(
355            DoubleDouble::from_bit_pair(Q[3]),
356            x,
357            DoubleDouble::from_bit_pair(Q[2]),
358        );
359        let q2 = DoubleDouble::mul_f64_add(
360            DoubleDouble::from_bit_pair(Q[5]),
361            x,
362            DoubleDouble::from_bit_pair(Q[4]),
363        );
364        let q3 = DoubleDouble::mul_f64_add(
365            DoubleDouble::from_bit_pair(Q[7]),
366            x,
367            DoubleDouble::from_bit_pair(Q[6]),
368        );
369        let q4 = DoubleDouble::mul_f64_add(
370            DoubleDouble::from_bit_pair(Q[9]),
371            x,
372            DoubleDouble::from_bit_pair(Q[8]),
373        );
374
375        let r0 = DoubleDouble::mul_add(x2, q1, q0);
376        let r1 = DoubleDouble::mul_add(x2, q3, q2);
377
378        let s0 = DoubleDouble::mul_add(x4, r1, r0);
379        let s1 = DoubleDouble::mul_add(x2, DoubleDouble::from_bit_pair(Q[10]), q4);
380        let p_den = DoubleDouble::mul_add(x8, s1, s0);
381
382        let v = DoubleDouble::div(p_num, p_den);
383        return v.to_f64();
384    }
385
386    let mut erfcx_abs_x = core_erfcx(f64::from_bits(ax));
387    if x < 0. {
388        // exp(x^2)erfc(-x) = 2*exp(x^2) - erfcx(|x|)
389        erfcx_abs_x = DoubleDouble::from_exact_add(erfcx_abs_x.hi, erfcx_abs_x.lo);
390        let d2x = DoubleDouble::from_exact_mult(x, x);
391        let expd2x = exp_dd_fast(d2x);
392        return DoubleDouble::mul_f64_add(expd2x, 2., -erfcx_abs_x).to_f64();
393    }
394    erfcx_abs_x.to_f64()
395}
396
397#[cfg(test)]
398mod tests {
399    use crate::f_erfcx;
400
401    #[test]
402    fn test_erfcx() {
403        assert_eq!(f_erfcx(2.2204460492503131e-18), 1.0);
404        assert_eq!(f_erfcx(-2.2204460492503131e-18), 1.0);
405        assert_eq!(f_erfcx(-f64::EPSILON), 1.0000000000000002);
406        assert_eq!(f_erfcx(f64::EPSILON), 0.9999999999999998);
407        assert_eq!(f_erfcx(-173.), f64::INFINITY);
408        assert_eq!(f_erfcx(-9.4324165432), 8.718049147018359e38);
409        assert_eq!(f_erfcx(9.4324165432), 0.059483265496416374);
410        assert_eq!(f_erfcx(-1.32432512125), 11.200579112797806);
411        assert_eq!(f_erfcx(1.32432512125), 0.3528722004785406);
412        assert_eq!(f_erfcx(-0.532431235), 2.0560589406595384);
413        assert_eq!(f_erfcx(0.532431235), 0.5994337293294584);
414        assert_eq!(f_erfcx(1e-26), 1.0);
415        assert_eq!(f_erfcx(-0.500000000023073), 1.952360489253639);
416        assert_eq!(f_erfcx(-175.), f64::INFINITY);
417        assert_eq!(f_erfcx(f64::INFINITY), 0.);
418        assert_eq!(f_erfcx(f64::NEG_INFINITY), f64::INFINITY);
419        assert!(f_erfcx(f64::NAN).is_nan());
420    }
421}
pxfm/err/erfcx.rs

pxfm/err/
erfcx.rs
