pxfm/err/
inverfc.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::double_double::DoubleDouble;
31use crate::logs::{fast_log_d_to_dd, fast_log_dd};
32use crate::polyeval::{f_polyeval4, f_polyeval5};
33
34#[cold]
35fn inverf_0p06_to_0p75(x: DoubleDouble) -> DoubleDouble {
36    // First step rational approximant is generated, but it's ill-conditioned, thus
37    // we're using taylor expansion to create Newton form at the point.
38    // Generated in Wolfram Mathematica:
39    // <<FunctionApproximations`
40    // ClearAll["Global`*"]
41    // f[x_]:=InverseErf[x]/x
42    // g[x_] =f[Sqrt[x]];
43    // {err0,approx}=MiniMaxApproximation[g[z],{z,{0.06,0.75},9,9},WorkingPrecision->75, MaxIterations->100]
44    // num=Numerator[approx][[1]];
45    // den=Denominator[approx][[1]];
46    // poly=den;
47    // coeffs=CoefficientList[poly,z];
48    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
49    // x0=SetPrecision[0.5625,75];
50    // NumberForm[Series[num[x],{x,x0,50}], ExponentFunction->(Null&)]
51    // coeffs=Table[SeriesCoefficient[num[x],{x,x0,k}],{k,0,9}];
52    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]];
53    const P: [(u64, u64); 10] = [
54        (0xbc3e06eda42202a0, 0x3f93c2fc5d00e0c8),
55        (0xbc6eb374406b33b4, 0xbfc76fcfd022e3ff),
56        (0xbc857822d7ffd282, 0x3fe6f8443546010a),
57        (0x3c68269c66dfb28a, 0xbff80996754ceb79),
58        (0x3c543dce8990a9f9, 0x3ffcf778d5ef0504),
59        (0xbc72fc55f73765f6, 0xbff433be821423d0),
60        (0xbc66d05fb37c8592, 0x3fdf15f19e9d8da4),
61        (0x3c56dfb85e83a2c5, 0xbfb770b6827e0829),
62        (0x3bff1472ecdfa403, 0x3f7a98a2980282bb),
63        (0x3baffb33d69d6276, 0xbf142a246fd2c07c),
64    ];
65    let x2 = DoubleDouble::quick_mult(x, x);
66    let vz = DoubleDouble::full_add_f64(x2, -0.5625);
67
68    let vx2 = vz * vz;
69    let vx4 = vx2 * vx2;
70    let vx8 = vx4 * vx4;
71
72    let p0 = DoubleDouble::mul_add(
73        vz,
74        DoubleDouble::from_bit_pair(P[1]),
75        DoubleDouble::from_bit_pair(P[0]),
76    );
77    let p1 = DoubleDouble::mul_add(
78        vz,
79        DoubleDouble::from_bit_pair(P[3]),
80        DoubleDouble::from_bit_pair(P[2]),
81    );
82    let p2 = DoubleDouble::mul_add(
83        vz,
84        DoubleDouble::from_bit_pair(P[5]),
85        DoubleDouble::from_bit_pair(P[4]),
86    );
87    let p3 = DoubleDouble::mul_add(
88        vz,
89        DoubleDouble::from_bit_pair(P[7]),
90        DoubleDouble::from_bit_pair(P[6]),
91    );
92    let p4 = DoubleDouble::mul_add(
93        vz,
94        DoubleDouble::from_bit_pair(P[9]),
95        DoubleDouble::from_bit_pair(P[8]),
96    );
97
98    let q0 = DoubleDouble::mul_add(vx2, p1, p0);
99    let q1 = DoubleDouble::mul_add(vx2, p3, p2);
100
101    let r0 = DoubleDouble::mul_add(vx4, q1, q0);
102    let num = DoubleDouble::mul_add(vx8, p4, r0);
103    // Generated in Wolfram Mathematica:
104    // <<FunctionApproximations`
105    // ClearAll["Global`*"]
106    // f[x_]:=InverseErf[x]/x
107    // g[x_] =f[Sqrt[x]];
108    // {err0,approx}=MiniMaxApproximation[g[z],{z,{0.06,0.75},9,9},WorkingPrecision->75, MaxIterations->100]
109    // num=Numerator[approx][[1]];
110    // den=Denominator[approx][[1]];
111    // coeffs=CoefficientList[poly,z];
112    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
113    // x0=SetPrecision[0.5625,75];
114    // NumberForm[Series[den[x],{x,x0,50}], ExponentFunction->(Null&)]
115    // coeffs=Table[SeriesCoefficient[den[x],{x,x0,k}],{k,0,9}];
116    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]];
117    const Q: [(u64, u64); 10] = [
118        (0xbc36337f24e57cb9, 0x3f92388d5d757e3a),
119        (0xbc63dfae43d60e0b, 0xbfc6ca7da581358c),
120        (0xbc77656389bd0e62, 0x3fe7c82ce417b4e0),
121        (0xbc93679667bef2f0, 0xbffad58651fd1a51),
122        (0x3ca2c6cb9eb17fb4, 0x4001bdb67e93a242),
123        (0xbc9b58961ba253bc, 0xbffbdaeff6fbb81c),
124        (0x3c7861f549c6aa61, 0x3fe91b12cf47da3a),
125        (0xbc696dfd665b2f5e, 0xbfc7c5d0ffb7f1da),
126        (0x3c1552b0ec0ba7b3, 0x3f939ada247f7609),
127        (0xbbcaa226fb7b30a8, 0xbf41be65038ccfe6),
128    ];
129
130    let p0 = DoubleDouble::mul_add(
131        vz,
132        DoubleDouble::from_bit_pair(Q[1]),
133        DoubleDouble::from_bit_pair(Q[0]),
134    );
135    let p1 = DoubleDouble::mul_add(
136        vz,
137        DoubleDouble::from_bit_pair(Q[3]),
138        DoubleDouble::from_bit_pair(Q[2]),
139    );
140    let p2 = DoubleDouble::mul_add(
141        vz,
142        DoubleDouble::from_bit_pair(Q[5]),
143        DoubleDouble::from_bit_pair(Q[4]),
144    );
145    let p3 = DoubleDouble::mul_add(
146        vz,
147        DoubleDouble::from_bit_pair(Q[7]),
148        DoubleDouble::from_bit_pair(Q[6]),
149    );
150    let p4 = DoubleDouble::mul_add(
151        vz,
152        DoubleDouble::from_bit_pair(Q[9]),
153        DoubleDouble::from_bit_pair(Q[8]),
154    );
155
156    let q0 = DoubleDouble::mul_add(vx2, p1, p0);
157    let q1 = DoubleDouble::mul_add(vx2, p3, p2);
158
159    let r0 = DoubleDouble::mul_add(vx4, q1, q0);
160    let den = DoubleDouble::mul_add(vx8, p4, r0);
161
162    let r = DoubleDouble::div(num, den);
163    DoubleDouble::quick_mult(r, x)
164}
165
166#[inline]
167fn inverf_asympt_small(z: DoubleDouble, zeta_sqrt: DoubleDouble) -> DoubleDouble {
168    // Generated in Wolfram Mathematica:
169    // <<FunctionApproximations`
170    // ClearAll["Global`*"]
171    // f[x_]:=InverseErf[Exp[-1/(x^2)]*(-1+Exp[1/(x^2)])]/(Sqrt[-Log[1-(Exp[-1/(x^2)]*(-1+Exp[1/(x^2)]))]] )
172    // {err0, approx,err1}=MiniMaxApproximation[f[z],{z,{0.2,0.9999999},10,10},WorkingPrecision->90]
173    // num=Numerator[approx];
174    // den=Denominator[approx];
175    // poly=num;
176    // coeffs=CoefficientList[poly,z];
177    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
178    const P: [(u64, u64); 11] = [
179        (0x3c936555853a8b2c, 0x3ff0001df06a2515),
180        (0x3cea488e802db3c3, 0x404406ba373221da),
181        (0xbce27d42419754e3, 0x407b0442e38a9597),
182        (0xbd224a407624cbdf, 0x409c9277e31ef446),
183        (0x3d4f16ce65d6fea0, 0x40aec3ec005b1d8a),
184        (0x3d105bc37bc61b58, 0x40b46be8f860f4d9),
185        (0x3d5ca133dcdecaa0, 0x40b3826e6a32dad7),
186        (0x3d1d52013ba8aa38, 0x40aae93a603cf3ea),
187        (0xbd07a75306df0fc3, 0x4098ab8357dc2e51),
188        (0x3d1bb6770bb7a27e, 0x407ebead00879010),
189        (0xbbfcbff4a9737936, 0x3f8936117ccbff83),
190    ];
191
192    let z2 = DoubleDouble::quick_mult(z, z);
193    let z4 = DoubleDouble::quick_mult(z2, z2);
194    let z8 = DoubleDouble::quick_mult(z4, z4);
195
196    let q0 = DoubleDouble::mul_add(
197        DoubleDouble::from_bit_pair(P[1]),
198        z,
199        DoubleDouble::from_bit_pair(P[0]),
200    );
201    let q1 = DoubleDouble::mul_add(
202        DoubleDouble::from_bit_pair(P[3]),
203        z,
204        DoubleDouble::from_bit_pair(P[2]),
205    );
206    let q2 = DoubleDouble::mul_add(
207        DoubleDouble::from_bit_pair(P[5]),
208        z,
209        DoubleDouble::from_bit_pair(P[4]),
210    );
211    let q3 = DoubleDouble::mul_add(
212        DoubleDouble::from_bit_pair(P[7]),
213        z,
214        DoubleDouble::from_bit_pair(P[6]),
215    );
216    let q4 = DoubleDouble::mul_add(
217        DoubleDouble::from_bit_pair(P[9]),
218        z,
219        DoubleDouble::from_bit_pair(P[8]),
220    );
221
222    let r0 = DoubleDouble::mul_add(z2, q1, q0);
223    let r1 = DoubleDouble::mul_add(z2, q3, q2);
224
225    let s0 = DoubleDouble::mul_add(z4, r1, r0);
226    let s1 = DoubleDouble::mul_add(z2, DoubleDouble::from_bit_pair(P[10]), q4);
227    let num = DoubleDouble::mul_add(z8, s1, s0);
228
229    // See numerator generation above:
230    // poly=den;
231    // coeffs=CoefficientList[poly,z];
232    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
233    const Q: [(u64, u64); 11] = [
234        (0x0000000000000000, 0x3ff0000000000000),
235        (0xbc75b1109d4a3262, 0x40440782efaab17f),
236        (0x3d1f7775b207d84f, 0x407b2da74b0d39f2),
237        (0xbd3291fdbab49501, 0x409dac8d9e7c90b2),
238        (0xbd58d8fdd27707a9, 0x40b178dfeffa3192),
239        (0xbd57fc74ad705ce0, 0x40bad19b686f219f),
240        (0x3d4075510031f2cd, 0x40be70a598208cea),
241        (0xbd5442e109152efb, 0x40b9683ef36ae330),
242        (0x3d5398192933962e, 0x40b04b7c4c3ca8ee),
243        (0x3d2d04d03598e303, 0x409bd0080799fbf1),
244        (0x3d2a988eb552ef44, 0x40815a46f12bafe3),
245    ];
246
247    let q0 = DoubleDouble::mul_add_f64(
248        DoubleDouble::from_bit_pair(Q[1]),
249        z,
250        f64::from_bits(0x3ff0000000000000),
251    );
252    let q1 = DoubleDouble::mul_add(
253        DoubleDouble::from_bit_pair(Q[3]),
254        z,
255        DoubleDouble::from_bit_pair(Q[2]),
256    );
257    let q2 = DoubleDouble::mul_add(
258        DoubleDouble::from_bit_pair(Q[5]),
259        z,
260        DoubleDouble::from_bit_pair(Q[4]),
261    );
262    let q3 = DoubleDouble::mul_add(
263        DoubleDouble::from_bit_pair(Q[7]),
264        z,
265        DoubleDouble::from_bit_pair(Q[6]),
266    );
267    let q4 = DoubleDouble::mul_add(
268        DoubleDouble::from_bit_pair(Q[9]),
269        z,
270        DoubleDouble::from_bit_pair(Q[8]),
271    );
272
273    let r0 = DoubleDouble::mul_add(z2, q1, q0);
274    let r1 = DoubleDouble::mul_add(z2, q3, q2);
275
276    let s0 = DoubleDouble::mul_add(z4, r1, r0);
277    let s1 = DoubleDouble::mul_add(z2, DoubleDouble::from_bit_pair(Q[10]), q4);
278    let den = DoubleDouble::mul_add(z8, s1, s0);
279    let r = DoubleDouble::div(num, den);
280    DoubleDouble::quick_mult(r, zeta_sqrt)
281}
282
283// branch for |x| > 0.9999 for extreme tail
284#[cold]
285fn inverf_asympt_long(z: DoubleDouble, zeta_sqrt: DoubleDouble) -> DoubleDouble {
286    // First step rational approximant is generated, but it's ill-conditioned, thus
287    // we're using taylor expansion to create Newton form at the point.
288    // Generated in Wolfram Mathematica:
289    // <<FunctionApproximations`
290    // ClearAll["Global`*"]
291    // f[x_]:=InverseErf[Exp[-1/(x^2)]*(-1+Exp[1/(x^2)])]/(Sqrt[-Log[1-(Exp[-1/(x^2)]*(-1+Exp[1/(x^2)]))]] )
292    // {err0, approx}=MiniMaxApproximation[f[z],{z,{0.2,0.9999999},13,13},WorkingPrecision->90]
293    // num=Numerator[approx][[1]];
294    // den=Denominator[approx][[1]];
295    // poly=num;
296    // coeffs=CoefficientList[poly,z];
297    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
298    const P: [(u64, u64); 14] = [
299        (0x3c97612f9b24a614, 0x3ff0000ba84cc7a5),
300        (0xbcee8fe2da463412, 0x40515246546f5d88),
301        (0x3d2fa4a2b891b526, 0x40956b6837159b11),
302        (0x3d5d673ffad4f817, 0x40c5a1aa3be58652),
303        (0x3d8867a1e5506f88, 0x40e65ebb1e1e7c75),
304        (0xbd9bbc0764ed8f5b, 0x40fd2064a652e5c2),
305        (0xbda78e569c0d237f, 0x410a385c627c461c),
306        (0xbdab3123ebc465d7, 0x4110f05ca2b65fe5),
307        (0x3d960def35955192, 0x4110bb079af2fe08),
308        (0xbd97904816054836, 0x410911c24610c11c),
309        (0xbd937745e9192593, 0x40fc603244adca35),
310        (0xbd65fbc476d63050, 0x40e6399103188c21),
311        (0xbd61016ef381cce6, 0x40c6482b44995b89),
312        (0x3c326105c49e5a1a, 0xbfab44bd8b4e3138),
313    ];
314
315    let z2 = z * z;
316    let z4 = z2 * z2;
317    let z8 = z4 * z4;
318
319    let g0 = DoubleDouble::mul_add(
320        z,
321        DoubleDouble::from_bit_pair(P[1]),
322        DoubleDouble::from_bit_pair(P[0]),
323    );
324    let g1 = DoubleDouble::mul_add(
325        z,
326        DoubleDouble::from_bit_pair(P[3]),
327        DoubleDouble::from_bit_pair(P[2]),
328    );
329    let g2 = DoubleDouble::mul_add(
330        z,
331        DoubleDouble::from_bit_pair(P[5]),
332        DoubleDouble::from_bit_pair(P[4]),
333    );
334    let g3 = DoubleDouble::mul_add(
335        z,
336        DoubleDouble::from_bit_pair(P[7]),
337        DoubleDouble::from_bit_pair(P[6]),
338    );
339    let g4 = DoubleDouble::mul_add(
340        z,
341        DoubleDouble::from_bit_pair(P[9]),
342        DoubleDouble::from_bit_pair(P[8]),
343    );
344    let g5 = DoubleDouble::mul_add(
345        z,
346        DoubleDouble::from_bit_pair(P[11]),
347        DoubleDouble::from_bit_pair(P[10]),
348    );
349    let g6 = DoubleDouble::mul_add(
350        z,
351        DoubleDouble::from_bit_pair(P[13]),
352        DoubleDouble::from_bit_pair(P[12]),
353    );
354
355    let h0 = DoubleDouble::mul_add(z2, g1, g0);
356    let h1 = DoubleDouble::mul_add(z2, g3, g2);
357    let h2 = DoubleDouble::mul_add(z2, g5, g4);
358
359    let q0 = DoubleDouble::mul_add(z4, h1, h0);
360    let q1 = DoubleDouble::mul_add(z4, g6, h2);
361
362    let num = DoubleDouble::mul_add(z8, q1, q0);
363
364    // See numerator generation above:
365    // poly=den;
366    // coeffs=CoefficientList[poly,z];
367    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
368    const Q: [(u64, u64); 14] = [
369        (0x0000000000000000, 0x3ff0000000000000),
370        (0xbcfc7b886ee61417, 0x405152838f711f3c),
371        (0xbd33f933c14e831a, 0x409576cb78cab36e),
372        (0x3d33fb09e2c4898a, 0x40c5e8a2c7602ced),
373        (0x3d7be430c664bf7e, 0x40e766fdc8c7638c),
374        (0x3dac662e74cdfc0e, 0x4100276b5f47b5f1),
375        (0x3da67d06e82a8495, 0x410f843887f8a24a),
376        (0x3dbbf2e22fc2550a, 0x4116d04271703e08),
377        (0xbdb2fb3aed100853, 0x4119aff4ed32b74b),
378        (0x3dba75e7b7171c3c, 0x4116b5eb8bf386bd),
379        (0x3dab2d8b8c1937eb, 0x410f71c38e84cb34),
380        (0xbda4e2e8a50b7370, 0x4100ca04b0f36b94),
381        (0xbd86ed6df34fdaf9, 0x40e9151ded4cf4b7),
382        (0x3d6938ea702c0328, 0x40c923ee1ab270c4),
383    ];
384
385    let g0 = DoubleDouble::mul_add(
386        z,
387        DoubleDouble::from_bit_pair(Q[1]),
388        DoubleDouble::from_bit_pair(Q[0]),
389    );
390    let g1 = DoubleDouble::mul_add(
391        z,
392        DoubleDouble::from_bit_pair(Q[3]),
393        DoubleDouble::from_bit_pair(Q[2]),
394    );
395    let g2 = DoubleDouble::mul_add(
396        z,
397        DoubleDouble::from_bit_pair(Q[5]),
398        DoubleDouble::from_bit_pair(Q[4]),
399    );
400    let g3 = DoubleDouble::mul_add(
401        z,
402        DoubleDouble::from_bit_pair(Q[7]),
403        DoubleDouble::from_bit_pair(Q[6]),
404    );
405    let g4 = DoubleDouble::mul_add(
406        z,
407        DoubleDouble::from_bit_pair(Q[9]),
408        DoubleDouble::from_bit_pair(Q[8]),
409    );
410    let g5 = DoubleDouble::mul_add(
411        z,
412        DoubleDouble::from_bit_pair(Q[11]),
413        DoubleDouble::from_bit_pair(Q[10]),
414    );
415    let g6 = DoubleDouble::mul_add(
416        z,
417        DoubleDouble::from_bit_pair(Q[13]),
418        DoubleDouble::from_bit_pair(Q[12]),
419    );
420
421    let h0 = DoubleDouble::mul_add(z2, g1, g0);
422    let h1 = DoubleDouble::mul_add(z2, g3, g2);
423    let h2 = DoubleDouble::mul_add(z2, g5, g4);
424
425    let q0 = DoubleDouble::mul_add(z4, h1, h0);
426    let q1 = DoubleDouble::mul_add(z4, g6, h2);
427
428    let den = DoubleDouble::mul_add(z8, q1, q0);
429    let r = DoubleDouble::div(num, den);
430
431    DoubleDouble::quick_mult(r, zeta_sqrt)
432}
433
434#[inline]
435fn erf_core(x: DoubleDouble) -> DoubleDouble {
436    // x is always positive, here, should be filtered out before the call
437
438    if x.hi <= 0.0095 {
439        // 0.0095
440        // for small |x| using taylor series first 3 terms
441        // Generated by SageMath:
442        // from mpmath import mp, erf
443        //
444        // mp.prec = 100
445        //
446        // def inverf_series(n_terms):
447        //     from mpmath import taylor
448        //     series_erf = taylor(mp.erfinv, 0, n_terms)
449        //     return series_erf
450        //
451        // ser = inverf_series(10)
452        // for i in range(1, len(ser), 2):
453        //     k = ser[i]
454        //     print("f64::from_bits(" + double_to_hex(RealField(100)(k)) + "),")
455        let z2 = DoubleDouble::quick_mult(x, x);
456        let p = f_fmla(
457            z2.hi,
458            f64::from_bits(0x3fb62847c47dda48),
459            f64::from_bits(0x3fc053c2c0ab91c5),
460        );
461        let mut r = DoubleDouble::mul_f64_add(
462            z2,
463            p,
464            DoubleDouble::from_bit_pair((0xbc33ea2ef8dde075, 0x3fcdb29fb2fee5e4)),
465        );
466        r = DoubleDouble::mul_add(
467            z2,
468            r,
469            DoubleDouble::from_bit_pair((0xbc8618f13eb7ca89, 0x3fec5bf891b4ef6b)),
470        );
471        // (rh + rl) * z = rh * z + rl*z
472        let v = DoubleDouble::quick_mult(r, x);
473        return v;
474    } else if x.hi <= 0.06 {
475        // 0.06
476        // for |x| < 0.06 using taylor series first 5 terms
477        // Generated by SageMath:
478        // from mpmath import mp, erf
479        //
480        // mp.prec = 100
481        //
482        // def inverf_series(n_terms):
483        //     from mpmath import taylor
484        //     series_erf = taylor(mp.erfinv, 0, n_terms)
485        //     return series_erf
486        //
487        // ser = inverf_series(10)
488        // for i in range(1, len(ser), 2):
489        //     k = ser[i]
490        //     print("f64::from_bits(" + double_to_hex(RealField(100)(k)) + "),")
491        let z2 = DoubleDouble::quick_mult(x, x);
492        let p = f_polyeval4(
493            z2.hi,
494            f64::from_bits(0x3fb62847c47dda48),
495            f64::from_bits(0x3fb0a13189c6ef7a),
496            f64::from_bits(0x3faa7c85c89bb08b),
497            f64::from_bits(0x3fa5eeb1d488e312),
498        );
499        let mut r = DoubleDouble::mul_f64_add(
500            z2,
501            p,
502            DoubleDouble::from_bit_pair((0x3c2cec68daff0d80, 0x3fc053c2c0ab91c5)),
503        );
504        r = DoubleDouble::mul_add(
505            z2,
506            r,
507            DoubleDouble::from_bit_pair((0xbc33ea2ef8dde075, 0x3fcdb29fb2fee5e4)),
508        );
509        r = DoubleDouble::mul_add(
510            z2,
511            r,
512            DoubleDouble::from_bit_pair((0xbc8618f13eb7ca89, 0x3fec5bf891b4ef6b)),
513        );
514        // (rh + rl) * z = rh * z + rl*z
515        let v = DoubleDouble::quick_mult(r, x);
516        return v;
517    }
518
519    if x.hi <= 0.75 {
520        // |x| < 0.75
521
522        // First step rational approximant is generated, but it's ill-conditioned, thus
523        // we're using taylor expansion to create Newton form at the point.
524        // Generated in Wolfram Mathematica:
525        // <<FunctionApproximations`
526        // ClearAll["Global`*"]
527        // f[x_]:=InverseErf[x]/x
528        // g[x_] =f[Sqrt[x]];
529        // {err0,approx}=MiniMaxApproximation[g[z],{z,{0.06,0.75},9,9},WorkingPrecision->75, MaxIterations->100]
530        // num=Numerator[approx][[1]];
531        // den=Denominator[approx][[1]];
532        // poly=den;
533        // coeffs=CoefficientList[poly,z];
534        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
535        // x0=SetPrecision[0.5625,75];
536        // NumberForm[Series[num[x],{x,x0,50}], ExponentFunction->(Null&)]
537        // coeffs=Table[SeriesCoefficient[num[x],{x,x0,k}],{k,0,9}];
538        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]];
539        const P: [(u64, u64); 5] = [
540            (0xbc3e06eda42202a0, 0x3f93c2fc5d00e0c8),
541            (0xbc6eb374406b33b4, 0xbfc76fcfd022e3ff),
542            (0xbc857822d7ffd282, 0x3fe6f8443546010a),
543            (0x3c68269c66dfb28a, 0xbff80996754ceb79),
544            (0x3c543dce8990a9f9, 0x3ffcf778d5ef0504),
545        ];
546        let x2 = DoubleDouble::quick_mult(x, x);
547        let vz = DoubleDouble::full_add_f64(x2, -0.5625);
548        let ps_num = f_polyeval5(
549            vz.hi,
550            f64::from_bits(0xbff433be821423d0),
551            f64::from_bits(0x3fdf15f19e9d8da4),
552            f64::from_bits(0xbfb770b6827e0829),
553            f64::from_bits(0x3f7a98a2980282bb),
554            f64::from_bits(0xbf142a246fd2c07c),
555        );
556        let mut num = DoubleDouble::mul_f64_add(vz, ps_num, DoubleDouble::from_bit_pair(P[4]));
557        num = DoubleDouble::mul_add(vz, num, DoubleDouble::from_bit_pair(P[3]));
558        num = DoubleDouble::mul_add(vz, num, DoubleDouble::from_bit_pair(P[2]));
559        num = DoubleDouble::mul_add(vz, num, DoubleDouble::from_bit_pair(P[1]));
560        num = DoubleDouble::mul_add(vz, num, DoubleDouble::from_bit_pair(P[0]));
561
562        // Generated in Wolfram Mathematica:
563        // <<FunctionApproximations`
564        // ClearAll["Global`*"]
565        // f[x_]:=InverseErf[x]/x
566        // g[x_] =f[Sqrt[x]];
567        // {err0,approx}=MiniMaxApproximation[g[z],{z,{0.06,0.75},9,9},WorkingPrecision->75, MaxIterations->100]
568        // num=Numerator[approx][[1]];
569        // den=Denominator[approx][[1]];
570        // coeffs=CoefficientList[poly,z];
571        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
572        // x0=SetPrecision[0.5625,75];
573        // NumberForm[Series[den[x],{x,x0,50}], ExponentFunction->(Null&)]
574        // coeffs=Table[SeriesCoefficient[den[x],{x,x0,k}],{k,0,9}];
575        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]];
576        const Q: [(u64, u64); 5] = [
577            (0xbc36337f24e57cb9, 0x3f92388d5d757e3a),
578            (0xbc63dfae43d60e0b, 0xbfc6ca7da581358c),
579            (0xbc77656389bd0e62, 0x3fe7c82ce417b4e0),
580            (0xbc93679667bef2f0, 0xbffad58651fd1a51),
581            (0x3ca2c6cb9eb17fb4, 0x4001bdb67e93a242),
582        ];
583
584        let ps_den = f_polyeval5(
585            vz.hi,
586            f64::from_bits(0xbffbdaeff6fbb81c),
587            f64::from_bits(0x3fe91b12cf47da3a),
588            f64::from_bits(0xbfc7c5d0ffb7f1da),
589            f64::from_bits(0x3f939ada247f7609),
590            f64::from_bits(0xbf41be65038ccfe6),
591        );
592
593        let mut den = DoubleDouble::mul_f64_add(vz, ps_den, DoubleDouble::from_bit_pair(Q[4]));
594        den = DoubleDouble::mul_add(vz, den, DoubleDouble::from_bit_pair(Q[3]));
595        den = DoubleDouble::mul_add(vz, den, DoubleDouble::from_bit_pair(Q[2]));
596        den = DoubleDouble::mul_add(vz, den, DoubleDouble::from_bit_pair(Q[1]));
597        den = DoubleDouble::mul_add(vz, den, DoubleDouble::from_bit_pair(Q[0]));
598        let r = DoubleDouble::div(num, den);
599        let k = DoubleDouble::quick_mult(r, x);
600        let err = f_fmla(
601            k.hi,
602            f64::from_bits(0x3c70000000000000), // 2^-56
603            f64::from_bits(0x3c40000000000000), // 2^-59
604        );
605        let ub = k.hi + (k.lo + err);
606        let lb = k.hi + (k.lo - err);
607        if ub == lb {
608            return k;
609        }
610        return inverf_0p06_to_0p75(x);
611    }
612
613    let q = DoubleDouble::full_add_f64(-x, 1.0);
614
615    let mut zeta = fast_log_dd(q);
616    zeta = DoubleDouble::from_exact_add(zeta.hi, zeta.lo);
617    zeta = -zeta;
618    let zeta_sqrt = zeta.fast_sqrt();
619    let rz = zeta_sqrt.recip();
620
621    if x.hi < 0.9999 {
622        inverf_asympt_small(rz, zeta_sqrt)
623    } else {
624        inverf_asympt_long(rz, zeta_sqrt)
625    }
626}
627
628#[cold]
629fn inverfc_extra_small(x: f64) -> DoubleDouble {
630    // Reversed order for erfinv with direct identity without subtraction.
631    let q = x;
632
633    let mut zeta = fast_log_d_to_dd(q);
634    zeta = DoubleDouble::from_exact_add(zeta.hi, zeta.lo);
635    zeta = -zeta;
636    let zeta_sqrt = zeta.fast_sqrt();
637    let rz = zeta_sqrt.recip();
638    if x >= 0.0001 {
639        inverf_asympt_small(rz, zeta_sqrt)
640    } else {
641        inverf_asympt_long(rz, zeta_sqrt)
642    }
643}
644
645/// Complementary inverse error function
646pub fn f_erfcinv(x: f64) -> f64 {
647    let ix = x.to_bits();
648
649    if ix >= 0x4000000000000000u64 || ix == 0 {
650        // |x| == NaN, x == inf, |x| == 0, x < 0
651        if ix.wrapping_shl(1) == 0 {
652            return f64::INFINITY;
653        }
654        if ix == 0x4000000000000000u64 {
655            return f64::NEG_INFINITY;
656        }
657        return f64::NAN; // x == NaN, x == Inf, x > 2
658    }
659
660    if x == 1. {
661        return 0.;
662    }
663
664    // we compute erfcinv through identity
665    // erfcinv(x) = -erfinv(1-x)
666
667    static SIGN: [f64; 2] = [1.0, -1.0];
668
669    if x < 0.1 {
670        return inverfc_extra_small(x).to_f64();
671    }
672
673    let dx = if x > 1. {
674        DoubleDouble::from_full_exact_sub(2., x)
675    } else {
676        DoubleDouble::new(0., x)
677    };
678    let sign = SIGN[(x > 1.) as usize];
679
680    let mut dx = DoubleDouble::full_add_f64(-dx, 1.);
681    dx = DoubleDouble::from_exact_add(dx.hi, dx.lo);
682    erf_core(dx).to_f64() * sign
683}
684
685#[cfg(test)]
686mod tests {
687    use super::*;
688
689    #[test]
690    fn test_inverfc() {
691        assert_eq!(f_erfcinv(0.12), 1.0993909519492193);
692        assert_eq!(f_erfcinv(1.0000000000027623e-13), 5.261512368864527);
693        assert_eq!(f_erfcinv(1.0001200000182189), -0.00010634724760131264);
694        assert_eq!(f_erfcinv(0.7001200000182189), 0.2723481758403576);
695        assert_eq!(f_erfcinv(1.5231200000182189), -0.502985998867995);
696        assert_eq!(f_erfcinv(1.99545434324323243), -2.0064739778442213);
697        assert_eq!(f_erfcinv(1.), 0.);
698        assert!(f_erfcinv(2.05).is_nan());
699        assert!(f_erfcinv(-0.01).is_nan());
700        assert!(f_erfcinv(f64::NAN).is_nan());
701        assert!(f_erfcinv(f64::NEG_INFINITY).is_nan());
702        assert!(f_erfcinv(f64::INFINITY).is_nan());
703    }
704}
pxfm/err/inverfc.rs

pxfm/err/
inverfc.rs
