pxfm/err/inverff.rs
1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/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::logs::simple_fast_log;
30use crate::polyeval::{
31 f_estrin_polyeval7, f_estrin_polyeval8, f_estrin_polyeval9, f_polyeval3, f_polyeval5,
32 f_polyeval10, f_polyeval11,
33};
34
35#[inline]
36pub(crate) fn erfinv_core(z: f64, ax: u32, sign: f32) -> f32 {
37 if ax <= 0x3c1ba5e3u32 {
38 // 0.0095
39 // for small |x| using taylor series first 3 terms
40 let z2 = z * z;
41 // Generated by SageMath:
42 // from mpmath import mp, erf
43 //
44 // mp.prec = 100
45 //
46 // def inverf_series(n_terms):
47 // from mpmath import taylor
48 // series_erf = taylor(mp.erfinv, 0, n_terms)
49 // return series_erf
50 //
51 // ser = inverf_series(10)
52 // for i in range(1, len(ser), 2):
53 // k = ser[i]
54 // print("f64::from_bits(" + double_to_hex(RealField(100)(k)) + "),")
55 let p = f_polyeval3(
56 z2,
57 f64::from_bits(0x3fec5bf891b4ef6b),
58 f64::from_bits(0x3fcdb29fb2fee5e4),
59 f64::from_bits(0x3fc053c2c0ab91c5),
60 ) * z;
61 return f32::copysign(p as f32, sign);
62 } else if ax <= 0x3d75c28fu32 {
63 // 0.06
64 // for |x| < 0.06 using taylor series first 5 terms
65 let z2 = z * z;
66 // Generated by SageMath:
67 // from mpmath import mp, erf
68 //
69 // mp.prec = 100
70 //
71 // def inverf_series(n_terms):
72 // from mpmath import taylor
73 // series_erf = taylor(mp.erfinv, 0, n_terms)
74 // return series_erf
75 //
76 // ser = inverf_series(10)
77 // for i in range(1, len(ser), 2):
78 // k = ser[i]
79 // print("f64::from_bits(" + double_to_hex(RealField(100)(k)) + "),")
80 let p = f_polyeval5(
81 z2,
82 f64::from_bits(0x3fec5bf891b4ef6b),
83 f64::from_bits(0x3fcdb29fb2fee5e4),
84 f64::from_bits(0x3fc053c2c0ab91c5),
85 f64::from_bits(0x3fb62847c47dda48),
86 f64::from_bits(0x3fb0a13189c6ef7a),
87 ) * z;
88 return f32::copysign(p as f32, sign);
89 }
90
91 if ax <= 0x3f400000u32 {
92 // |x| <= 0.75
93 let z2 = z * z;
94
95 // First step rational approximant is generated, but it's ill-conditioned, thus
96 // we're using taylor expansion to create Newton form at the point.
97 //
98 // <<FunctionApproximations`
99 // ClearAll["Global`*"]
100 // f[x_]:=InverseErf[Sqrt[x]]/Sqrt[x]
101 // {err0, approx}=MiniMaxApproximation[f[z],{z,{0.06,0.75},8,7},WorkingPrecision->70]
102 // num=Numerator[approx][[1]];
103 // den=Denominator[approx][[1]];
104 // poly=num;
105 // coeffs=CoefficientList[poly,z];
106 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
107 let r = z2 - 0.5625;
108 // x0=SetPrecision[0.5625,75];
109 // NumberForm[Series[num[x],{x,x0,50}], ExponentFunction->(Null&)]
110 // coeffs=Table[SeriesCoefficient[num[x],{x,x0,k}],{k,0,8}];
111 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
112 let p_num = f_estrin_polyeval9(
113 r,
114 f64::from_bits(0x3fa329348a73d9d4),
115 f64::from_bits(0xbfd2cb089b644580),
116 f64::from_bits(0x3fed229149f732d6),
117 f64::from_bits(0xbff6a233d2028bff),
118 f64::from_bits(0x3ff268adbfbb6023),
119 f64::from_bits(0xbfddac401c7d70f4),
120 f64::from_bits(0x3fb3b1bd759d5046),
121 f64::from_bits(0xbf67aeb45bad547e),
122 f64::from_bits(0xbf01ccc7434d381b),
123 );
124 // x0=SetPrecision[0.5625,75];
125 // NumberForm[Series[den[x],{x,x0,50}], ExponentFunction->(Null&)]
126 // coeffs=Table[SeriesCoefficient[den[x],{x,x0,k}],{k,0,7}];
127 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
128 let p_den = f_estrin_polyeval8(
129 r,
130 f64::from_bits(0x3fa1aac2ee4b1413),
131 f64::from_bits(0xbfd279342e281c99),
132 f64::from_bits(0x3feef89a353c6d1b),
133 f64::from_bits(0xbffa8f1b7cd6d0a7),
134 f64::from_bits(0x3ff89ce6289819a1),
135 f64::from_bits(0xbfe7db5282a4a2e1),
136 f64::from_bits(0x3fc543f9a928db4a),
137 f64::from_bits(0xbf888fd2990e88db),
138 );
139 let k = (p_num / p_den) * z;
140 f32::copysign(k as f32, sign)
141 } else if ax <= 0x3f580000u32 {
142 // |x| <= 0.84375
143 let z2 = z * z;
144
145 // First step rational approximant is generated, but it's ill-conditioned, thus
146 // we're using taylor expansion to create Newton form at the point.
147 //
148 // <<FunctionApproximations`
149 // ClearAll["Global`*"]
150 // f[x_]:=InverseErf[Sqrt[x]]/Sqrt[x]
151 // {err0, approx}=MiniMaxApproximation[f[z],{z,{0.75,0.84375},6,6},WorkingPrecision->70]
152 // num=Numerator[approx][[1]];
153 // den=Denominator[approx][[1]];
154 // poly=num;
155 // coeffs=CoefficientList[poly,z];
156 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
157 let r = z2 - 0.84375;
158 // x0=SetPrecision[0.84375,75];
159 // NumberForm[Series[num[x],{x,x0,50}], ExponentFunction->(Null&)]
160 // coeffs=Table[SeriesCoefficient[num[x],{x,x0,k}],{k,0,9}];
161 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
162 let p_num = f_polyeval10(
163 r,
164 f64::from_bits(0x3f116d07e62cbb74),
165 f64::from_bits(0xbf5c38d390052412),
166 f64::from_bits(0x3f92d6f96f84efe3),
167 f64::from_bits(0xbfbac9189cae446b),
168 f64::from_bits(0x3fd5dd124fb25677),
169 f64::from_bits(0xbfe49845d46b80ab),
170 f64::from_bits(0x3fe556c4913f60f8),
171 f64::from_bits(0xbfd59e527704e33b),
172 f64::from_bits(0x3fb07614a5e6c9f1),
173 f64::from_bits(0xbf60ce54a2d8a789),
174 );
175 // x0=SetPrecision[0.84375,75];
176 // NumberForm[Series[den[x],{x,x0,50}], ExponentFunction->(Null&)]
177 // coeffs=Table[SeriesCoefficient[den[x],{x,x0,k}],{k,0,9}];
178 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
179 let p_den = f_polyeval10(
180 r,
181 f64::from_bits(0x3f09fbdd1c987d1e),
182 f64::from_bits(0xbf5602ad17d419f4),
183 f64::from_bits(0x3f8efe31ea5bc71d),
184 f64::from_bits(0xbfb77e5f1bd26730),
185 f64::from_bits(0x3fd4c3f03e4f5478),
186 f64::from_bits(0xbfe5aa87dfc5e757),
187 f64::from_bits(0x3fe9c6406f9abc0b),
188 f64::from_bits(0xbfdff2f008b4db05),
189 f64::from_bits(0x3fc1123be5319800),
190 f64::from_bits(0xbf83be49c2d5cb9e),
191 );
192 let k = (p_num / p_den) * z;
193 f32::copysign(k as f32, sign)
194 } else if ax <= 0x3f700000u32 {
195 // |x| <= 0.9375
196 // First step rational approximant is generated, but it's ill-conditioned, thus
197 // we're using taylor expansion to create Newton form at the point.
198 //
199 // <<FunctionApproximations`
200 // ClearAll["Global`*"]
201 // f[x_]:=InverseErf[Sqrt[x]]/Sqrt[x]
202 // {err0, approx}=MiniMaxApproximation[f[z],{z,{0.84375,0.9375},10,9},WorkingPrecision->70]
203 // num=Numerator[approx][[1]];
204 // den=Denominator[approx][[1]];
205 // coeffs=CoefficientList[poly,z];
206 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
207 let x2 = z * z;
208 let r = x2 - 0.87890625;
209 // x0=SetPrecision[0.87890625,75];
210 // NumberForm[Series[num[x],{x,x0,50}], ExponentFunction->(Null&)]
211 // coeffs=Table[SeriesCoefficient[num[x],{x,x0,k}],{k,0,9}];
212 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
213 let p_num = f_polyeval11(
214 r,
215 f64::from_bits(0x3ec70f1cbf8a758b),
216 f64::from_bits(0xbf1c9dff87b698d0),
217 f64::from_bits(0x3f5dfe7be00cc21c),
218 f64::from_bits(0xbf913fd09c5a3682),
219 f64::from_bits(0x3fb7ab0095693976),
220 f64::from_bits(0xbfd3b3ca6a3c9919),
221 f64::from_bits(0x3fe3533be6d1d8c8),
222 f64::from_bits(0xbfe48208ef308ac7),
223 f64::from_bits(0x3fd361a82dab69d1),
224 f64::from_bits(0xbfa2401965a98195),
225 f64::from_bits(0xbf54ba4d14ca54e3),
226 );
227 // x0=SetPrecision[0.87890625,75];
228 // NumberForm[Series[den[x],{x,x0,50}], ExponentFunction->(Null&)]
229 // coeffs=Table[SeriesCoefficient[den[x],{x,x0,k}],{k,0,9}];
230 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
231 let p_den = f_polyeval10(
232 r,
233 f64::from_bits(0x3ec0699f391e2327),
234 f64::from_bits(0xbf151ec184941078),
235 f64::from_bits(0x3f5717bb379a3c6e),
236 f64::from_bits(0xbf8beed3755c3484),
237 f64::from_bits(0x3fb46148b4a431ef),
238 f64::from_bits(0xbfd25690b7bc76fa),
239 f64::from_bits(0x3fe3f1b2f4ee0d9d),
240 f64::from_bits(0xbfe888a7a4511975),
241 f64::from_bits(0x3fdd84db18f2a240),
242 f64::from_bits(0xbfb844807521be56),
243 );
244 let f = z * (p_num / p_den);
245 f32::copysign(f as f32, sign)
246 } else {
247 // Rational approximation generated by Wolfram Mathematica:
248 // for inverf(x) = sqrt(-log(1-x))*R(1/sqrt(-log(1-x)))
249 //
250 // <<FunctionApproximations`
251 // ClearAll["Global`*"]
252 // f[x_]:=InverseErf[Exp[-1/(x^2)]*(-1+Exp[1/(x^2)])]/(Sqrt[-Log[1-(Exp[-1/(x^2)]*(-1+Exp[1/(x^2)]))]] )
253 // {err0, approx,err1}=MiniMaxApproximation[f[z],{z,{0.2,0.9999999},7,6},WorkingPrecision->90]
254 // num=Numerator[approx];
255 // den=Denominator[approx];
256 // poly=num;
257 // coeffs=CoefficientList[poly,z];
258 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
259 // poly=den;
260 // coeffs=CoefficientList[poly,z];
261 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
262 let zeta = -simple_fast_log(1. - z);
263 let zeta_sqrt = zeta.sqrt();
264 let rcp_zeta = (1. / zeta) * zeta_sqrt;
265 let p_num = f_estrin_polyeval8(
266 rcp_zeta,
267 f64::from_bits(0x3ff00072876c578e),
268 f64::from_bits(0x40314e00c10282da),
269 f64::from_bits(0x404f4a1412af03f6),
270 f64::from_bits(0x404c895cc0d9b1b3),
271 f64::from_bits(0x404545794620bfaf),
272 f64::from_bits(0x403264d21ea21354),
273 f64::from_bits(0x3fc5a5141dd19237),
274 f64::from_bits(0xbf8c2e49707c21ec),
275 );
276 let p_den = f_estrin_polyeval7(
277 rcp_zeta,
278 f64::from_bits(0x3ff0000000000000),
279 f64::from_bits(0x403151312c313d77),
280 f64::from_bits(0x405032345fa3d0cd),
281 f64::from_bits(0x4053e0a81d4c5f09),
282 f64::from_bits(0x4054fa20c5e0731c),
283 f64::from_bits(0x404620d7f94d4804),
284 f64::from_bits(0x4035d7400867b81f),
285 );
286 let r = zeta_sqrt * (p_num / p_den);
287 f32::copysign(r as f32, sign)
288 }
289}
290
291/// Inverse error function
292///
293/// Max ulp 0.5
294pub fn f_erfinvf(x: f32) -> f32 {
295 let ax = x.to_bits() & 0x7fff_ffff;
296
297 if ax >= 0x3f800000u32 || ax <= 0x3400_0000u32 {
298 // |x| >= 1, |x| == 0, |x| <= f32::EPSILON
299 if ax == 0 {
300 // |x| == 0
301 return 0.;
302 }
303
304 if ax <= 0x3400_0000u32 {
305 // |x| <= f32::EPSILON
306 // inverf(x) ~ Sqrt[Pi]x/2+O[x]^3
307 const SQRT_PI_OVER_2: f64 = f64::from_bits(0x3fec5bf891b4ef6b);
308 return (x as f64 * SQRT_PI_OVER_2) as f32;
309 }
310
311 if ax == 0x3f800000u32 {
312 // |x| == 1
313 return if x.is_sign_negative() {
314 f32::NEG_INFINITY
315 } else {
316 f32::INFINITY
317 };
318 }
319
320 // |x| > 1
321 return f32::NAN; // |x| == NaN, |x| == Inf, |x| > 1
322 }
323
324 let z = f32::from_bits(ax) as f64;
325 erfinv_core(z, ax, x)
326}
327
328#[cfg(test)]
329mod tests {
330 use super::*;
331
332 #[test]
333 fn f_test_inv_erff() {
334 assert!(f_erfinvf(f32::NEG_INFINITY).is_nan());
335 assert!(f_erfinvf(f32::INFINITY).is_nan());
336 assert!(f_erfinvf(-1.1).is_nan());
337 assert!(f_erfinvf(1.1).is_nan());
338 assert_eq!(f_erfinvf(f32::EPSILON), 1.05646485e-7);
339 assert_eq!(f_erfinvf(-1.), f32::NEG_INFINITY);
340 assert_eq!(f_erfinvf(1.), f32::INFINITY);
341 assert_eq!(f_erfinvf(0.002), 0.0017724558);
342 assert_eq!(f_erfinvf(-0.002), -0.0017724558);
343 assert_eq!(f_erfinvf(0.02), 0.017726395);
344 assert_eq!(f_erfinvf(-0.02), -0.017726395);
345 assert_eq!(f_erfinvf(0.05), 0.044340387);
346 assert_eq!(f_erfinvf(-0.05), -0.044340387);
347 assert_eq!(f_erfinvf(0.5), 0.47693628);
348 assert_eq!(f_erfinvf(-0.5), -0.47693628);
349 assert_eq!(f_erfinvf(0.76), 0.8308411);
350 assert_eq!(f_erfinvf(-0.76), -0.8308411);
351 assert_eq!(f_erfinvf(0.92), 1.2379221);
352 assert_eq!(f_erfinvf(-0.92), -1.2379221);
353 assert_eq!(f_erfinvf(0.97), 1.5344859);
354 assert_eq!(f_erfinvf(-0.97), -1.5344859);
355 assert_eq!(f_erfinvf(0.99), 1.8213866);
356 assert_eq!(f_erfinvf(-0.99), -1.8213866);
357 assert_eq!(f_erfinvf(0.7560265), 0.82385886);
358 }
359}