thermite_special/specialized/
pd.rs

1use thermite::{
2    math::{
3        TranscendentalMathWithPolicy,
4        policy::{
5            DenormalBehavior, PrecisionPolicy,
6            policies::{CheckOverflow, ExtraPrecision, WorstPrecision},
7        },
8        specialized::SpecializedTranscendentalMath,
9        specialized::reference::{is_reference, map1, map1x2},
10    },
11    prelude::*,
12};
13
14use crate::RealSpecialMathWithPolicy as _;
15
16use super::*;
17
18impl<V: FloatVectorWithBits<Element = f64>> SpecializedSpecialMath<f64> for V
19where
20    V: TranscendentalMathWithPolicy<Element = f64>,
21    V: SpecializedTranscendentalMath<f64>,
22    // Pins the projection so `V`'s real-special methods (whose impl requires
23    // `Primal = V`) resolve. A type parameter's `Primal` will not normalize
24    // through the blanket impl on its own.
25    V: thermite::math::PrimalProjection<Primal = V>,
26    V: thermite::math::RealMathWithPolicy<Element = f64>,
27{
28    /// Through the strict `FloatVectorWithBits::two_sum`, so it survives `algebraic-scalar`.
29    #[inline(always)]
30    fn exp_two_sum(a: Self, b: Self) -> (Self, Self) {
31        a.two_sum(b)
32    }
33
34    // `EXACT_FMA = true`: a real vector's multiply-add is a single rounding, so the seed's
35    // `x*x` residual is real. See `generic::hermite::seed`.
36    #[inline(always)]
37    fn hermite_function_n<P: Policy, const N: usize>(mut x: Self) -> Self {
38        #[cfg(not(target_arch = "spirv"))]
39        if let Some(new_x) = thermite::math::specialized::FlushDenormals::<P>::flush_denormals([x]) {
40            x = new_x[0];
41        }
42
43        generic::hermite::hermite_function_n::<P, _, _, N, true>(x)
44    }
45
46    #[inline(always)]
47    fn hermite_function<P: Policy>(mut x: Self, n: u32) -> Self {
48        #[cfg(not(target_arch = "spirv"))]
49        if let Some(new_x) = thermite::math::specialized::FlushDenormals::<P>::flush_denormals([x]) {
50            x = new_x[0];
51        }
52
53        generic::hermite::hermite_function::<P, _, _, true>(x, n)
54    }
55
56    #[inline(always)]
57    fn hermite_function_series_n<P: Policy, const N: usize>(self, coeffs: &[Self::Element; N]) -> Self {
58        generic::hermite::hermite_function_series::<P, _, _, N, true>(self, coeffs)
59    }
60
61    #[inline(always)]
62    fn hermite_function_series<P: Policy>(self, coeffs: &[Self::Element]) -> Self {
63        generic::hermite::hermite_function_series_slice::<P, _, _, true>(self, coeffs)
64    }
65
66    #[inline(always)]
67    fn zetac<P: Policy>(self) -> Self {
68        generic::zeta::zeta_impl::<P, _, _, true>(self)
69    }
70
71    #[inline(always)]
72    fn polylog<P: Policy>(self, order: crate::PolylogOrder<f64, i64>) -> Self {
73        generic::polylog::polylog_impl::<P, f64, Self>(self, order)
74    }
75
76    #[inline(always)]
77    fn zeta<P: Policy>(self) -> Self {
78        generic::zeta::zeta_impl::<P, _, _, false>(self)
79    }
80
81    #[inline(always)]
82    fn zeta_with_deriv<P: Policy, const ZETAC: bool>(self) -> (Self, Self) {
83        generic::zeta::zeta_core::<P, _, _, ZETAC, true>(self)
84    }
85
86    #[inline(always)]
87    fn bessel_i<P: Policy, const N: i32>(self) -> Self {
88        bessel_i_dispatch::<P, Self, N, false>(self)
89    }
90
91    #[inline(always)]
92    fn bessel_i_scaled<P: Policy, const N: i32>(self) -> Self {
93        bessel_i_dispatch::<P, Self, N, true>(self)
94    }
95
96    #[inline(always)]
97    fn bessel_k<P: Policy, const N: i32>(self) -> Self {
98        bessel_k_dispatch::<P, Self, N, false>(self)
99    }
100
101    #[inline(always)]
102    fn bessel_k_scaled<P: Policy, const N: i32>(self) -> Self {
103        bessel_k_dispatch::<P, Self, N, true>(self)
104    }
105
106    #[inline(always)]
107    fn bessel_j<P: Policy, const N: i32>(self) -> Self {
108        use crate::tables::bessel::jy::{BESSEL_J0_F64, BESSEL_J1_F64};
109        // `Reference` is contractually bit-identical to libm, lane by lane, and unlike most
110        // of this crate, libm actually has these (the C/POSIX XSI set) at every order, so the
111        // arm exists. `I`/`K` have no libm counterpart and therefore no reference arm.
112        if const { is_reference::<P>() } {
113            let v = if const { N == 0 } {
114                map1(self, libm::j0)
115            } else if const { N.unsigned_abs() == 1 } {
116                map1(self, libm::j1)
117            } else {
118                map1(self, |v| libm::jn(N.abs(), v))
119            };
120            // A sign flip is exact, so reflecting libm's own value keeps the tier's
121            // bit-identity promise rather than trading it for a second algorithm.
122            return if const { bessel_reflect_negates(N) } { -v } else { v };
123        }
124        let v = if const { N == 0 } {
125            generic::bessel::jy::bessel_j0_impl::<P, f64, _, _, _, _>(self, &BESSEL_J0_F64)
126        } else if const { N.unsigned_abs() == 1 } {
127            generic::bessel::jy::bessel_j1_impl::<P, f64, _, _, _, _>(self, &BESSEL_J1_F64)
128        } else {
129            generic::bessel::jy::bessel_jn_pair_impl::<P, f64, _, _, _, _, _, _, _, N>(
130                self,
131                &BESSEL_J0_F64,
132                &BESSEL_J1_F64,
133            )
134            .1
135        };
136        // `J_{-n} = (-1)^n J_n`. Every arm above evaluated at `|N|`.
137        if const { bessel_reflect_negates(N) } { -v } else { v }
138    }
139
140    #[inline(always)]
141    fn bessel_y<P: Policy, const N: i32>(self) -> Self {
142        use crate::tables::bessel::jy::{BESSEL_J0_F64, BESSEL_J1_F64, BESSEL_Y0_F64, BESSEL_Y1_F64};
143        if const { is_reference::<P>() } {
144            let v = if const { N == 0 } {
145                map1(self, libm::y0)
146            } else if const { N.unsigned_abs() == 1 } {
147                map1(self, libm::y1)
148            } else {
149                map1(self, |v| libm::yn(N.abs(), v))
150            };
151            return if const { bessel_reflect_negates(N) } { -v } else { v };
152        }
153        let v = if const { N.unsigned_abs() >= 2 } {
154            // Y is the dominant solution, so upward recurrence is stable and costs exactly
155            // |N|-1 steps. No trip count question at all, unlike J.
156            let y0 = generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, false>(
157                self,
158                &BESSEL_Y0_F64,
159                &BESSEL_J0_F64,
160            );
161            let y1 = generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, true>(
162                self,
163                &BESSEL_Y1_F64,
164                &BESSEL_J1_F64,
165            );
166            generic::bessel::jy::bessel_yn_recur::<f64, _, N>(self, y0, y1).1
167        } else if const { N == 0 } {
168            generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, false>(
169                self,
170                &BESSEL_Y0_F64,
171                &BESSEL_J0_F64,
172            )
173        } else {
174            generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, true>(
175                self,
176                &BESSEL_Y1_F64,
177                &BESSEL_J1_F64,
178            )
179        };
180        // `Y_{-n} = (-1)^n Y_n`, the same reflection `J` gets.
181        if const { bessel_reflect_negates(N) } { -v } else { v }
182    }
183
184    #[inline(always)]
185    fn bessel_i_with_deriv<P: Policy, const N: i32, const SCALED: bool>(self) -> (Self, Self) {
186        bessel_i_deriv_dispatch::<P, Self, N, SCALED>(self)
187    }
188
189    #[inline(always)]
190    fn bessel_k_with_deriv<P: Policy, const N: i32, const SCALED: bool>(self) -> (Self, Self) {
191        bessel_k_deriv_dispatch::<P, Self, N, SCALED>(self)
192    }
193
194    #[inline(always)]
195    fn bessel_iv<P: Policy, const SCALED: bool>(self, order: crate::BesselOrder<Self, Self::Signed>) -> Self {
196        // Half-integer order is elementary: hyperbolic seeds and the same two recurrence
197        // directions the integer kernel uses. See `generic::bessel_half`.
198        let order = order.simplify();
199        if let crate::BesselOrder::HalfInteger(k) = order {
200            return generic::bessel::half::bessel_ik_half::<P, f64, _, SCALED>(
201                Self::from_signed_integer(k) * Self::HALF,
202                self,
203                crate::tables::bessel::BESSEL_I0_F64.far_threshold,
204            )
205            .0;
206        }
207        // `Thirds` and `Real` take the table-free arms in `generic::bessel_ik_nu`.
208        let Some(n) = order.as_integer() else {
209            return generic::bessel::ik_real::bessel_ik_real::<P, f64, _, _, 25, 25, SCALED, true>(
210                order.to_real(),
211                self,
212                &crate::tables::lgamma1p::LGAMMA1P_F64,
213                crate::tables::bessel::BESSEL_I0_F64.far_threshold,
214            )
215            .0;
216        };
217        // `I_{-n} = I_n` for integer `n`, so only the magnitude matters and no sign is owed
218        // afterwards. `J`/`Y` below are the ones that reflect.
219        let nf = Self::from_signed_integer(n).abs();
220        let v = generic::bessel::ik::bessel_iv_impl::<P, f64, _, _, _, _, SCALED>(
221            self,
222            nf,
223            &crate::tables::bessel::BESSEL_I0_F64,
224        );
225        // Orders 0 and 1 have closed forms, and the ratio ladder is measurably worse at them:
226        // it reaches order 1 as `I_0 * r_1`, paying the continued fraction for a value the
227        // table gives directly. Measured 5.49 ULP against 3.04 before this select was added.
228        let i1 =
229            generic::bessel::ik::bessel_i1_impl::<P, _, _, _, _, SCALED>(self, &crate::tables::bessel::BESSEL_I1_F64);
230        nf.cmp_le(Self::ONE).select(
231            nf.cmp_le(Self::ZERO).select(
232                generic::bessel::ik::bessel_i0_impl::<P, _, _, _, _, SCALED>(
233                    self,
234                    &crate::tables::bessel::BESSEL_I0_F64,
235                ),
236                i1,
237            ),
238            v,
239        )
240    }
241
242    #[inline(always)]
243    fn bessel_kv<P: Policy, const SCALED: bool>(self, order: crate::BesselOrder<Self, Self::Signed>) -> Self {
244        // Half-integer order is elementary (see `generic::bessel_half`).
245        let order = order.simplify();
246        if let crate::BesselOrder::HalfInteger(k) = order {
247            return generic::bessel::half::bessel_ik_half::<P, f64, _, SCALED>(
248                Self::from_signed_integer(k) * Self::HALF,
249                self,
250                crate::tables::bessel::BESSEL_I0_F64.far_threshold,
251            )
252            .1;
253        }
254        // `Thirds` and `Real` take the table-free arms in `generic::bessel_ik_nu`.
255        let Some(n) = order.as_integer() else {
256            // `NEED_I = false`: this entry wants only `K`, which is the cheap half. Skipping
257            // `I` skips the continued fraction and the asymptotic series both.
258            return generic::bessel::ik_real::bessel_ik_real::<P, f64, _, _, 25, 25, SCALED, false>(
259                order.to_real(),
260                self,
261                &crate::tables::lgamma1p::LGAMMA1P_F64,
262                crate::tables::bessel::BESSEL_I0_F64.far_threshold,
263            )
264            .1;
265        };
266        // `K_{-n} = K_n`, as with `I`.
267        let nf = Self::from_signed_integer(n).abs();
268        generic::bessel::ik::bessel_kv_impl::<P, f64, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, SCALED>(
269            self,
270            nf,
271            &crate::tables::bessel::BESSEL_K0_F64,
272            &crate::tables::bessel::BESSEL_K1_F64,
273            &crate::tables::bessel::BESSEL_I0_F64,
274            &crate::tables::bessel::BESSEL_I1_F64,
275        )
276    }
277
278    #[inline(always)]
279    fn bessel_jv<P: Policy>(self, order: crate::BesselOrder<Self, Self::Signed>) -> Self {
280        // Half-integer order is elementary (see `generic::bessel_half`). `simplify` has
281        // already turned an even numerator into `Integer`, so anything still tagged
282        // `HalfInteger` here is a genuine half-odd order.
283        let order = order.simplify();
284        if let crate::BesselOrder::HalfInteger(k) = order {
285            return generic::bessel::half::bessel_jy_half::<P, f64, _>(Self::from_signed_integer(k) * Self::HALF, self)
286                .0;
287        }
288        // `Thirds` and `Real` take the table-free arms in `generic::bessel_nu`, which cover
289        // the whole axis at any real order. Thirds are not specialised beyond that, and
290        // deliberately: see the module docs there.
291        let Some(n) = order.as_integer() else {
292            let nu = order.to_real();
293            return generic::bessel::jy_real::bessel_jy_real::<P, f64, _, 29, 25, 25, 17, 1>(
294                nu,
295                self,
296                Self::ZERO,
297                &crate::tables::lgamma1p::LGAMMA1P_F64,
298            )
299            .0;
300        };
301        // The const form routes `Reference` to libm. So must this one, or the tier silently
302        // stops meaning "bit-identical to libm" as soon as the order moves into a register.
303        if const { is_reference::<P>() } {
304            let mut out = self;
305            let mut i = 0;
306            while i < Self::LANES {
307                // Reflected here rather than handed to libm signed, so the tier means the
308                // same thing at negative order as the const form does.
309                let k = n.extractv(i);
310                let r = libm::jn(k.unsigned_abs() as i32, self.extractv(i));
311                out = out.insertv(i, if k < 0 && k % 2 != 0 { -r } else { r });
312                i += 1;
313            }
314            return out;
315        }
316        let (nf, flip) = bessel_reflect_v(Self::from_signed_integer(n));
317        generic::bessel::jy::bessel_jv_impl::<P, f64, _, _, _, _, _, _, _>(
318            self,
319            nf,
320            &crate::tables::bessel::jy::BESSEL_J0_F64,
321            &crate::tables::bessel::jy::BESSEL_J1_F64,
322        )
323        .neg_c(flip)
324    }
325
326    #[inline(always)]
327    fn bessel_yv<P: Policy>(self, order: crate::BesselOrder<Self, Self::Signed>) -> Self {
328        // Half-integer order is elementary (see `generic::bessel_half`).
329        let order = order.simplify();
330        if let crate::BesselOrder::HalfInteger(k) = order {
331            return generic::bessel::half::bessel_jy_half::<P, f64, _>(Self::from_signed_integer(k) * Self::HALF, self)
332                .1;
333        }
334        let Some(n) = order.as_integer() else {
335            let nu = order.to_real();
336            return generic::bessel::jy_real::bessel_jy_real::<P, f64, _, 29, 25, 25, 17, 1>(
337                nu,
338                self,
339                Self::ZERO,
340                &crate::tables::lgamma1p::LGAMMA1P_F64,
341            )
342            .1;
343        };
344        if const { is_reference::<P>() } {
345            let mut out = self;
346            let mut i = 0;
347            while i < Self::LANES {
348                let k = n.extractv(i);
349                let r = libm::yn(k.unsigned_abs() as i32, self.extractv(i));
350                out = out.insertv(i, if k < 0 && k % 2 != 0 { -r } else { r });
351                i += 1;
352            }
353            return out;
354        }
355        let (nf, flip) = bessel_reflect_v(Self::from_signed_integer(n));
356        generic::bessel::jy::bessel_yv_impl::<P, f64, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _>(
357            self,
358            nf,
359            &crate::tables::bessel::jy::BESSEL_Y0_F64,
360            &crate::tables::bessel::jy::BESSEL_Y1_F64,
361            &crate::tables::bessel::jy::BESSEL_J0_F64,
362            &crate::tables::bessel::jy::BESSEL_J1_F64,
363        )
364        .neg_c(flip)
365    }
366
367    impl_sph_bessel_entries!(f64, crate::tables::bessel::BESSEL_I0_F64);
368
369    impl_airy_entries!(
370        f64,
371        29,
372        25,
373        25,
374        17,
375        1,
376        &crate::tables::lgamma1p::LGAMMA1P_F64,
377        &crate::tables::bessel::airy::AIRY_ZERO_F64,
378        crate::tables::bessel::BESSEL_I0_F64
379    );
380
381    #[inline(always)]
382    fn bessel_j_with_deriv<P: Policy, const N: i32>(self) -> (Self, Self) {
383        // Order N-1 comes from the recurrence, which walks through it either way: forward
384        // passes it on the last step, downward keeps the shorter product.
385        let (prev, v) = if const { N == 0 } {
386            // J_{-1} = -J_1, so the identity still holds and the N/x term simply vanishes.
387            (-Self::bessel_j::<P, 1>(self), Self::bessel_j::<P, 0>(self))
388        } else if const { N.unsigned_abs() == 1 } {
389            (Self::bessel_j::<P, 0>(self), Self::bessel_j::<P, 1>(self))
390        } else {
391            generic::bessel::jy::bessel_jn_pair_impl::<P, f64, _, _, _, _, _, _, _, N>(
392                self,
393                &crate::tables::bessel::jy::BESSEL_J0_F64,
394                &crate::tables::bessel::jy::BESSEL_J1_F64,
395            )
396        };
397        let d = if const { N == 0 } {
398            prev
399        } else {
400            prev - v * (Self::splat(N.unsigned_abs() as f64) / self)
401        };
402        // The pair above is at `|N|`. Reflecting a negative order scales the function by a
403        // constant `(-1)^n`, so differentiating both sides carries the identical sign.
404        if const { bessel_reflect_negates(N) } {
405            (-v, -d)
406        } else {
407            (v, d)
408        }
409    }
410
411    #[inline(always)]
412    fn bessel_y_with_deriv<P: Policy, const N: i32>(self) -> (Self, Self) {
413        let (prev, v) = if const { N == 0 } {
414            (-Self::bessel_y::<P, 1>(self), Self::bessel_y::<P, 0>(self))
415        } else if const { N.unsigned_abs() == 1 } {
416            (Self::bessel_y::<P, 0>(self), Self::bessel_y::<P, 1>(self))
417        } else {
418            let y0 = generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, false>(
419                self,
420                &crate::tables::bessel::jy::BESSEL_Y0_F64,
421                &crate::tables::bessel::jy::BESSEL_J0_F64,
422            );
423            let y1 = generic::bessel::jy::bessel_y_impl::<P, f64, _, _, _, _, _, _, _, _, true>(
424                self,
425                &crate::tables::bessel::jy::BESSEL_Y1_F64,
426                &crate::tables::bessel::jy::BESSEL_J1_F64,
427            );
428            generic::bessel::jy::bessel_yn_recur::<f64, _, N>(self, y0, y1)
429        };
430        let d = if const { N == 0 } {
431            prev
432        } else {
433            prev - v * (Self::splat(N.unsigned_abs() as f64) / self)
434        };
435        // The pair above is at `|N|`. Reflecting a negative order scales the function by a
436        // constant `(-1)^n`, so differentiating both sides carries the identical sign.
437        if const { bessel_reflect_negates(N) } {
438            (-v, -d)
439        } else {
440            (v, d)
441        }
442    }
443
444    type ExpIntDetails = Self;
445    const LAGUERRE_PRODUCT_SEED_CAP: i32 = 170;
446
447    #[inline(always)]
448    fn chebyshev_n<P: Policy, const K: usize, const N: usize>(self, coeffs: &[f64; N]) -> Self {
449        // See the trait default: the kernel reads `N = 0` as "runtime length", so the
450        // empty-series rejection belongs to the entry point.
451        const {
452            assert!(N >= 1, "chebyshev_n: N must be at least 1");
453        }
454
455        // Real vectors have copysign and a real nearest endpoint, so the Reinsch form is
456        // available, but the kernel still gates it on the policy asking for `Best` or better.
457        generic::chebyshev::chebyshev_series::<P, _, _, K, N, true>(self, coeffs)
458    }
459
460    #[inline(always)]
461    fn chebyshev<P: Policy, const K: usize>(self, coeffs: &[f64]) -> Self {
462        // Reinsch available here too, on the same terms. See `chebyshev_n`.
463        generic::chebyshev::chebyshev_series::<P, _, _, K, 0, true>(self, coeffs)
464    }
465
466    // TEMP(bessel_j): disabled until orders beyond J_0 exist. See thermite-special/src/lib.rs.
467    //fn bessel_j<P: Policy, const N: i32>(self) -> Self {
468    //    todo!()
469    //}
470
471    #[inline(always)]
472    fn lambert_w<P: Policy>(self) -> (Self, Self) {
473        // Computes both W_0(x) and W_{-1}(x) simultaneously.
474        //
475        // Lambert W_0(x): principal branch, defined for x >= -1/e, returns values >= -1.
476        // Lambert W_{-1}(x): secondary real branch, defined for -1/e <= x < 0, returns values <= -1.
477        // Both satisfy w*e^w = x.
478        //
479        // Uses Halley's method with piecewise initial approximations, interleaving
480        // iterations for both branches to maximize instruction-level parallelism.
481        // f64 needs more iterations than f32 due to 52-bit mantissa.
482        //
483        // Halley's iteration for w*exp(w) = x:
484        //   ew = exp(w), f = w*ew - x, wp1 = w + 1
485        //   Denominator rewritten to avoid an extra division:
486        //     d = 2*wp1^2*ew - (w+2)*f
487        //   w' = w - 2*wp1*f / d
488
489        // For initial guess and first Halley iterations, use fast and loose precision
490        type Approx<P> = WorstPrecision<CheckOverflow<P, false>>;
491
492        let x = self;
493
494        // --- Initial approximation (piecewise) ---
495        //
496        // Branch-point region (x near -1/e): damped Puiseux series.
497        // See ps.rs lambert_w for full derivation.
498
499        let p0 = x.mul_adde(Self::E, Self::ONE); // ex + 1
500        let p = (p0 + p0).sqrt(); // sqrt(2(ex+1))
501
502        // p*(1 + p*(-1/3 + p*11/72))
503        let puiseux_numer = p * p.mul_adde(
504            p.mul_adde(
505                thermite::const_splat!(f64: 11.0 / 72.0),
506                thermite::const_splat!(f64: -1.0 / 3.0),
507            ),
508            Self::ONE,
509        );
510
511        // 1 + K*p_0*p
512        let puiseux_denom = p0.mul_adde(p * thermite::const_splat!(f64: 0.12991546098765432), Self::ONE);
513
514        let puiseux = puiseux_numer / puiseux_denom;
515
516        // W_0 branch: -1 + series, W_{-1} branch: -1 - series
517        let w0_branch = puiseux + Self::NEG_ONE;
518        let wm1_branch = Self::NEG_ONE - puiseux;
519
520        // W_0 middle region: ex/(2+ex), exact at x = -1/e and x = 0.
521        let ex = x * Self::E;
522        let w0_mid = ex / (Self::TWO + ex);
523
524        // Shared ln for asymptotic regions
525        let lnx = x.abs().ln_p::<Approx<P>>();
526
527        // W_0 asymptotic (x > e): L_1 - L_2 + L_2/L_1 where L_1 = ln(x), L_2 = ln(L_1).
528        // The L_2/L_1 correction is 0 at x = e (since L_2 = ln(1) = 0), so it doesn't
529        // overshoot near the transition, but closes the gap at large x.
530        let l2 = lnx.ln_p::<Approx<P>>();
531        let w0_asymptotic = (lnx - l2) + (l2 / lnx);
532
533        // W_{-1} asymptotic (x near 0^-): L_1 - L_2 where L_1 = ln(-x), L_2 = ln(-L_1)
534        // lnx = ln(|x|) = ln(-x) since x < 0; this is negative for small |x|.
535        // -lnx is positive, so (-lnx).ln() = ln(-ln(-x)) = L_2.
536        let wm1_asymptotic = lnx - (-lnx).ln_p::<Approx<P>>();
537
538        // Select initial guesses
539        let near_branch = x.cmp_lt(Self::splat(-0.1));
540        let large = x.cmp_gt(Self::E);
541        let mut w0 = near_branch.select(w0_branch, large.select(w0_asymptotic, w0_mid));
542
543        let near_branch_m1 = x.cmp_lt(Self::splat(-0.25));
544        let mut wm1 = near_branch_m1.select(wm1_branch, wm1_asymptotic);
545
546        // --- Interleaved Halley iterations ---
547        // Use cheap exp for warmup iteration, full-precision exp for the final ones.
548
549        #[inline(always)]
550        fn halley_step<P: Policy, W>(w: W, x: W) -> W
551        where
552            W: FloatVectorWithBits<Element = f64> + SpecializedTranscendentalMath<f64>,
553        {
554            // Use exp(-w) to avoid overflow/underflow in e^w for extreme w.
555            // g = w - x*e^{-w} = f*e^{-w}, d = (w^2+2w+2) + (w+2)*x*e^{-w}
556            // g and d are both single FMAs off enw, independent of each other.
557            let enw = (-w).exp_p::<P>();
558
559            let wp1 = w + W::ONE;
560            let q = wp1.mul_adde(wp1, W::ONE); // (w+1)^2 + 1 = w^2 + 2w + 2
561            let wp2h_x = wp1.mul_adde(x, x); // (w+2)*x - no exp dependency
562            let g = x.nmul_adde(enw, w); // w - x*e^{-w}
563            let d = wp2h_x.mul_adde(enw, q); // (w+2)*x*e^{-w} + (w^2+2w+2)
564            (wp1 + wp1).nmul_adde(g / d, w)
565        }
566
567        #[rustfmt::skip]
568        let num_iters = if const { P::POLICY.precision.ge(PrecisionPolicy::Best) } { 3 } else { 2 };
569
570        // warmup iteration with the looser precision to get close enough
571        // for the main iterations to converge in the target precision
572        w0 = halley_step::<Approx<P>, Self>(w0, x);
573        wm1 = halley_step::<Approx<P>, Self>(wm1, x);
574
575        let mut _iter = 0usize;
576        while _iter < num_iters {
577            _iter += 1;
578            w0 = halley_step::<CheckOverflow<P, false>, Self>(w0, x);
579            wm1 = halley_step::<CheckOverflow<P, false>, Self>(wm1, x);
580        }
581
582        // --- Edge cases ---
583        if const { P::POLICY.precision.ge(PrecisionPolicy::Average) } {
584            let x_is_zero = x.is_zero();
585
586            // At x = -1/e, both W_0 and W_{-1} = -1
587            w0 = x.cmp_eq(Self::FRAC_NEG_1_E).select(Self::NEG_ONE, w0);
588            w0 = w0.nz(x_is_zero); // W_0(0) = 0
589
590            wm1 = x.cmp_eq(Self::FRAC_NEG_1_E).select(Self::NEG_ONE, wm1);
591            wm1 = x_is_zero.select(Self::NEG_INFINITY, wm1); // W_{-1}(0) = -inf
592        }
593
594        if const { matches!(P::POLICY.denormal_behavior, DenormalBehavior::Preserve) } {
595            // for subnormal inputs, W_0(x) ≈ x
596            w0 = x.is_subnormal().select(x, w0);
597        }
598
599        if const { P::POLICY.check_overflow } {
600            let in_domain = x.cmp_ge(Self::FRAC_NEG_1_E);
601
602            // W_0 is undefined for x < -1/e, +inf -> +inf
603            w0 = in_domain.select(w0, Self::NAN);
604            w0 = x.cmp_eq(Self::INFINITY).select(Self::INFINITY, w0);
605
606            // W_{-1} is only defined for -1/e <= x < 0
607            wm1 = in_domain.select(wm1, Self::NAN);
608            wm1 = x.cmp_gt(Self::ZERO).select(Self::NAN, wm1);
609        }
610
611        (w0, wm1)
612    }
613
614    #[inline(always)]
615    #[allow(const_item_mutation)]
616    fn erf<P: Policy>(self) -> Self {
617        if const { is_reference::<P>() } {
618            return map1(self, libm::erf);
619        }
620
621        erf_d_internal::<Self, P, false, false>(self, &mut V::EMPTY)
622    }
623
624    #[inline(always)]
625    #[allow(const_item_mutation)]
626    fn erfc<P: Policy>(self) -> Self {
627        if const { is_reference::<P>() } {
628            return map1(self, libm::erfc);
629        }
630
631        erf_d_internal::<Self, P, true, false>(self, &mut V::EMPTY)
632    }
633
634    #[inline(always)]
635    fn erfcx<P: Policy>(self) -> Self {
636        // No libm counterpart at any tier. `erfcx` is not in the C library, and
637        // `exp(x*x) * erfc(x)` is exactly the overflowing form this replaces.
638        super::generic::erfcx::erfcx_internal::<Self, f64, P>(self)
639    }
640
641    #[inline(always)]
642    fn lgamma<P: Policy>(self) -> Self {
643        if const { is_reference::<P>() } {
644            return map1(self, libm::lgamma);
645        }
646
647        Self::lgamma_r::<P>(self).0
648    }
649
650    #[inline(always)]
651    fn tgamma<P: Policy>(self) -> Self {
652        if const { is_reference::<P>() } {
653            return map1(self, libm::tgamma);
654        }
655
656        let z = self;
657
658        if const { P::POLICY.precision.lt(PrecisionPolicy::Average) } {
659            // We have a good lgamma approximation, so use it for tgamma on lower precisions.
660            let (lgamma, sign) = z.lgamma_r_p::<P>();
661
662            // use min(P + 1, Average) precision here. We want decent precision,
663            // but not more than average.
664            return lgamma.exp_p::<ExtraPrecision<P>>() * sign;
665        }
666
667        // 172 is the largest integer whose factorial is finite in f64.
668        generic::gamma::tgamma_impl::<P, _, _, _>(
669            z,
670            &crate::tables::gamma::LANCZOS_F64,
671            172.0,
672            crate::tables::gamma::LN_MAX_F64,
673        )
674    }
675
676    #[inline(always)]
677    fn trigamma<P: Policy>(self) -> Self {
678        generic::trigamma::trigamma_impl::<P, _, _>(self, &crate::tables::gamma::TRIGAMMA_F64)
679    }
680
681    #[inline(always)]
682    fn polygamma<P: Policy>(self, n: u32) -> Self {
683        generic::polygamma::polygamma_impl::<P, _, _>(self, n)
684    }
685
686    #[inline(always)]
687    fn digamma<P: Policy>(self) -> Self {
688        generic::digamma::digamma_impl::<P, _, _, _, _, _, _>(self, &crate::tables::gamma::DIGAMMA_F64)
689    }
690
691    #[inline(always)]
692    fn beta<P: Policy>(a: Self, b: Self) -> Self {
693        generic::gamma::beta_impl::<P, _, _, _>(a, b, &crate::tables::gamma::LANCZOS_F64)
694    }
695
696    #[inline(always)]
697    fn expint_n<P: Policy, const N: usize>(self) -> Self {
698        generic::expint::expint_double_n::<P, f64, Self, N>(self)
699    }
700
701    #[inline(always)]
702    fn expint_primal_n<P: Policy, const N: usize>(self) -> (Self, Self) {
703        generic::expint::expint_double_primal_n::<P, f64, Self, N>(self)
704    }
705
706    #[inline(always)]
707    fn phi_n<P: Policy, const N: usize>(self) -> Self {
708        // Fixed series length for f64. See the f32 twin for the budget split.
709        let terms = const {
710            let needed =
711                super::generic::phi::phi_series_terms(N, f64::EPSILON * P::POLICY.precision.tolerance() as f64 / 32.0);
712            if needed < P::POLICY.max_iterations {
713                needed
714            } else {
715                P::POLICY.max_iterations
716            }
717        };
718        super::generic::phi::phi_internal_n::<Self, f64, P, N, false>(self, terms)
719    }
720
721    #[inline(always)]
722    fn expint<P: Policy>(self, n: u32) -> Self {
723        generic::expint::expint_double::<P, f64, Self>(self, n)
724    }
725
726    #[inline(always)]
727    fn expint_primal<P: Policy>(self, n: u32) -> (Self, Self) {
728        generic::expint::expint_double_primal::<P, f64, Self>(self, n)
729    }
730
731    #[inline(always)]
732    fn phi<P: Policy>(self, n: u32) -> Self {
733        // The const form's term counts, precomputed per policy. The search runs per call
734        // only past the table. Same budget split as `phi_n`.
735        const EPS_SCALE: f64 = 1.0 / 32.0;
736        let table = const {
737            super::generic::phi::phi_terms_table(
738                f64::EPSILON * P::POLICY.precision.tolerance() as f64 * EPS_SCALE,
739                P::POLICY.max_iterations,
740            )
741        };
742        let terms = match table.get(n as usize) {
743            Some(&t) => t,
744            None => {
745                let needed = super::generic::phi::phi_series_terms(
746                    n as usize,
747                    f64::EPSILON * P::POLICY.precision.tolerance() as f64 * EPS_SCALE,
748                );
749                if needed < P::POLICY.max_iterations {
750                    needed
751                } else {
752                    P::POLICY.max_iterations
753                }
754            }
755        };
756        super::generic::phi::phi_internal::<Self, f64, P, false>(self, n, terms)
757    }
758}
759
760impl<V: FloatVectorWithBits<Element = f64>> SpecializedRealSpecialMath<f64> for V
761where
762    V: TranscendentalMathWithPolicy<Element = f64>,
763    V: SpecializedTranscendentalMath<f64>,
764    // Pins the projection: a type parameter's `Primal` will not normalize through
765    // the blanket impl on its own, and the table signatures need `Primal = Self`.
766    V: thermite::math::PrimalProjection<Primal = V>,
767{
768    #[inline(always)]
769    fn fresnel<P: Policy>(self) -> (Self, Self) {
770        use crate::tables::fresnel as t;
771        generic::fresnel::fresnel_with::<P, _, _, _, _, _, _>(
772            self,
773            t::X0_F64,
774            t::MAP_F64,
775            t::CUTOFF_F64,
776            &t::CHEB_C_F64,
777            &t::CHEB_S_F64,
778            &t::AUX_P_F64,
779            &t::AUX_Q_F64,
780        )
781    }
782
783    #[inline(always)]
784    fn sici<P: Policy>(self) -> (Self, Self) {
785        use crate::tables::sici as t;
786        generic::sici::sici_with::<P, _, _, _, _, _, _>(
787            self,
788            t::X0_F64,
789            t::MAP_F64,
790            t::CUTOFF_F64,
791            &t::CHEB_SI_F64,
792            &t::CHEB_CIN_F64,
793            &t::AUX_P_F64,
794            &t::AUX_Q_F64,
795        )
796    }
797
798    // --- Spherical harmonics: the compile-time-table fast paths ---
799    //
800    // A concrete `f32`/`f64` element has a `ShConsts` table, which the generic
801    // defaults cannot assume. Both overrides are guarded by `L <= MAX_SH_DEGREE`,
802    // the extent of the stamped ladder, and fall back to the generic body above it.
803    // A statically-false `if const` arm is dropped before monomorphization, so the
804    // out-of-range table is never built.
805
806    #[inline(always)]
807    fn spherical_harmonics<P: Policy, const L: usize, const N: usize, const CS: bool>(
808        x: Self,
809        y: Self,
810        z: Self,
811        out: &mut [Self; N],
812    ) {
813        // Fully unrolled, constants folded into the instruction stream: no table is
814        // materialized at all, so there is nothing to hoist out of a loop. Above
815        // MAX_SH_DEGREE the kernel routes itself to the general path.
816        sh_impl::<P, f64, Self, L, N, CS>(x, y, z, out);
817    }
818
819    #[inline(always)]
820    fn spherical_harmonics_table<P: Policy, const L: usize, const N: usize, const CS: bool>(
821        table: &mut ShTable<Self, N>,
822    ) {
823        if const { L <= MAX_SH_DEGREE } {
824            // Every coefficient is already a compile-time constant of the right
825            // phase, so building the runtime table is a splat per entry, with none of
826            // the sqrt/divide work the generic default does.
827            let src = &<f64 as ShConsts<L, N, CS>>::TABLE;
828
829            let mut i = 0;
830            while i < N {
831                table.qmm[i] = Self::splat(src.qmm[i]);
832                table.em[i] = Self::splat(src.em[i]);
833                table.a[i] = Self::splat(src.a[i]);
834                table.nb[i] = Self::splat(src.nb[i]);
835                table.f[i] = Self::splat(src.f[i]);
836                table.mf[i] = Self::splat(src.mf[i]);
837                i += 1;
838            }
839        } else {
840            sh_table_impl::<Self, L, N, CS>(table);
841        }
842    }
843
844    #[inline(always)]
845    fn bessel_i_ratio<P: Policy>(self, nu: Self) -> Self {
846        generic::bessel::ratio::bessel_i_ratio_impl::<P, f64, Self>(self, nu)
847    }
848
849    #[inline(always)]
850    fn inv_bessel_i_ratio<P: Policy>(self, nu: Self) -> Self {
851        generic::bessel::ratio::inv_bessel_i_ratio_impl::<P, f64, Self>(self, nu)
852    }
853
854    #[inline(always)]
855    fn bessel_i_ratio_1m<P: Policy>(self, nu: Self) -> Self {
856        generic::bessel::ratio::bessel_i_ratio_1m_impl::<P, f64, Self>(self, nu)
857    }
858
859    #[inline(always)]
860    fn inv_bessel_i_ratio_1m<P: Policy>(self, nu: Self) -> Self {
861        generic::bessel::ratio::inv_bessel_i_ratio_1m_impl::<P, f64, Self>(self, nu)
862    }
863
864    #[inline(always)]
865    fn erfinv<P: Policy>(self) -> Self {
866        // Branchless erfinv: a cheap Winitzki seed refined with Halley iterations
867        // against the (accurate) erfc, which is far friendlier to SIMD than the
868        // many-branch piecewise-rational approach.
869        //
870        // We solve erfc(x) = q for x >= 0, where q = 1 - |y|. The Newton/Halley
871        // residual erf(x) - |y| is evaluated as q - erfc(x): in the tail both
872        // terms are tiny, so their difference keeps full relative precision (the
873        // direct form erf(x) - |y| would cancel two ~1 values down to noise).
874        // Halley is cubic, so the ~1% Winitzki seed reaches full f64 in 2 steps.
875        const ALPHA: f64 = 0.147;
876        const RCP_PI_ALPHA_2: f64 = 4.330746750799873; // 2 / (pi * ALPHA)
877        const RCP_ALPHA: f64 = 1.0 / ALPHA;
878        const SQRT_PI_2: f64 = 0.8862269254527580136490837416706; // sqrt(pi) / 2 = 1 / erf'(0)
879
880        let y = self.flush_denormals_p::<P>();
881        let a = y.abs();
882        let q = Self::ONE - a; // 1 - |y|
883        let omsq = q * (Self::ONE + a); // 1 - y^2, computed without cancellation near |y| = 1
884
885        // Winitzki seed (magnitude): sqrt(sqrt(t1^2 - ln(1-y^2)/alpha) - t1)
886        let lnv = omsq.ln_p::<P>(); // ln(1 - y^2) <= 0
887        let t1 = lnv.mul_adde(Self::HALF, Self::splat(RCP_PI_ALPHA_2));
888        let mut x = (t1.mul_adde(t1, lnv * Self::splat(-RCP_ALPHA)).sqrt() - t1).sqrt();
889
890        // Halley refinement: x -= u / (1 + x*u), u = (erf(x) - |y|) / erf'(x)
891        //   erf(x) - |y| = q - erfc(x),   1/erf'(x) = (sqrt(pi)/2) * exp(x^2)
892        let steps = if const { P::POLICY.precision.ge(PrecisionPolicy::Best) } {
893            2
894        } else {
895            1
896        };
897        let mut i = 0;
898        while i < steps {
899            // erfc(x) already computes exp(-x^2); reuse it so exp(x^2) is just a reciprocal.
900            let mut exp_neg = Self::EMPTY;
901            let erfc = erf_d_internal::<Self, P, true, true>(x, &mut exp_neg);
902            let u = (q - erfc) * Self::splat(SQRT_PI_2) / exp_neg;
903            x -= u / x.mul_adde(u, Self::ONE);
904            i += 1;
905        }
906
907        let mut res = x.copysign(y);
908
909        if const { P::POLICY.check_overflow } {
910            res = a.cmp_eq(Self::ONE).select(Self::INFINITY.copysign(y), res); // erfinv(+-1) = +-inf
911            res = a.cmp_gt(Self::ONE).select(Self::NAN, res); // out of domain
912        }
913
914        res
915    }
916
917    #[inline(always)]
918    fn lgamma_r<P: Policy>(self) -> (Self, Self) {
919        if const { is_reference::<P>() } {
920            // libm hands the sign back as an `i32`; this trait carries it as a float.
921            return map1x2(self, |x| {
922                let (v, s) = libm::lgamma_r(x);
923                (v, s as f64)
924            });
925        }
926
927        generic::gamma::lgamma_r_impl::<P, _, _, _>(self, &crate::tables::gamma::LANCZOS_F64)
928    }
929
930    /// Uses the algorithm from Peter John Acklam, sourced from here:
931    /// <https://web.archive.org/web/20151030215612/http://home.online.no/~pjacklam/notes/invnorm/>
932    #[inline(always)]
933    fn probit<P: Policy>(self) -> Self {
934        const A: [f64; 6] = [
935            2.506628277459239e+00,
936            -3.066479806614716e+01,
937            1.383577518672690e+02,
938            -2.759285104469687e+02,
939            2.209460984245205e+02,
940            -3.969683028665376e+01,
941        ];
942        const B: [f64; 6] = [
943            1.0,
944            -1.328068155288572e+01,
945            6.680131188771972e+01,
946            -1.556989798598866e+02,
947            1.615858368580409e+02,
948            -5.447609879822406e+01,
949        ];
950        const C: [f64; 6] = [
951            2.938163982698783e+00,
952            4.374664141464968e+00,
953            -2.549732539343734e+00,
954            -2.400758277161838e+00,
955            -3.223964580411365e-01,
956            -7.784894002430293e-03,
957        ];
958        const D: [f64; 5] = [
959            1.0,
960            3.754408661907416e+00,
961            2.445134137142996e+00,
962            3.224671290700398e-01,
963            7.784695709041462e-03,
964        ];
965
966        // f64: refine the Acklam estimate with one Halley step (REFINE = true).
967        generic::probit::probit_acklam::<P, _, _, true>(self, &A, &B, &C, &D)
968    }
969
970    #[inline(always)]
971    fn langevin<P: Policy>(self) -> Self {
972        // The Worst/Medium tiers take the short table (see it for its error).
973        if const { P::POLICY.precision.le(PrecisionPolicy::Medium) } {
974            generic::langevin::langevin_primal::<P, _, _, 11, false>(self, &LANGEVIN_SMALL_F64_LO).0
975        } else {
976            generic::langevin::langevin_primal::<P, _, _, 16, false>(self, &LANGEVIN_SMALL_F64).0
977        }
978    }
979
980    #[inline(always)]
981    fn langevin_1m<P: Policy>(self) -> Self {
982        // The Worst/Medium tiers take the short table (see it for its error).
983        if const { P::POLICY.precision.le(PrecisionPolicy::Medium) } {
984            generic::langevin::langevin_primal::<P, _, _, 11, true>(self, &LANGEVIN_SMALL_F64_LO).0
985        } else {
986            generic::langevin::langevin_primal::<P, _, _, 16, true>(self, &LANGEVIN_SMALL_F64).0
987        }
988    }
989
990    // f64 refines with Halley (see the kernel docs).
991    #[inline(always)]
992    fn inv_langevin<P: Policy>(self) -> Self {
993        generic::langevin::inv_langevin::<P, _, _, 16, 9, true, false>(self, &LANGEVIN_SMALL_F64, &LANGEVIN_SEED_F64)
994    }
995
996    #[inline(always)]
997    fn inv_langevin_1m<P: Policy>(self) -> Self {
998        generic::langevin::inv_langevin::<P, _, _, 16, 9, true, true>(self, &LANGEVIN_SMALL_F64, &LANGEVIN_SEED_F64)
999    }
1000
1001    // same form as f32
1002    #[inline(always)]
1003    fn gelu<P: Policy>(self, alpha: Self) -> Self {
1004        let x = self;
1005
1006        let alpha_x = alpha * x;
1007
1008        // GELU(x) = 0.5 * x * (1 + erf(ax / sqrt(2))) = 0.5 * x * erfc(-ax / sqrt(2))
1009        // O = false: skip the exp(-ax^2) byproduct that only the derivative needs.
1010        let mut unused = Self::EMPTY;
1011        let c = erf_d_internal::<V, P, true, false>(alpha_x * -Self::FRAC_1_SQRT_2, &mut unused);
1012
1013        (x * Self::HALF) * c
1014    }
1015}
1016
1017impl<V: FloatVectorWithBits<Element = f64>> SpecializedRealPrimalMath<f64> for V
1018where
1019    V: TranscendentalMathWithPolicy<Element = f64>,
1020    V: SpecializedTranscendentalMath<f64>,
1021    V: thermite::math::PrimalProjection<Primal = V>,
1022{
1023    #[inline(always)]
1024    fn langevin_d<P: Policy>(self) -> (Self, Self) {
1025        generic::langevin::langevin_primal::<P, _, _, 16, false>(self, &LANGEVIN_SMALL_F64)
1026    }
1027
1028    #[inline(always)]
1029    #[allow(clippy::too_many_arguments)]
1030    fn spherical_harmonics_d<P: Policy, const L: usize, const N: usize, const CS: bool>(
1031        x: Self,
1032        y: Self,
1033        z: Self,
1034        out: &mut [Self; N],
1035        ddx: &mut [Self; N],
1036        ddy: &mut [Self; N],
1037        ddz: &mut [Self; N],
1038    ) {
1039        sh_d_impl::<P, f64, Self, L, N, CS>(x, y, z, out, ddx, ddy, ddz);
1040    }
1041
1042    #[inline(always)]
1043    fn gelu_d<P: Policy>(self, alpha: Self) -> (Self, Self) {
1044        let x = self;
1045
1046        let alpha_x = alpha * x;
1047
1048        // 0.5 * x * erfc(-ax / sqrt(2))
1049        let mut exp_neg_ax2 = Self::EMPTY;
1050        let c = erf_d_internal::<V, P, true, true>(alpha_x * -Self::FRAC_1_SQRT_2, &mut exp_neg_ax2);
1051
1052        let half_c = c * Self::HALF; // 0.5 * (1 + erf(ax/sqrt(2)))
1053
1054        let y = x * half_c;
1055        let dy = if matches!(V::HAS_NATIVE_FMA, thermite::tribool::True) {
1056            (alpha_x * Self::FRAC_1_SQRT_TAU).mul_add(exp_neg_ax2, half_c)
1057        } else {
1058            half_c + alpha_x * Self::FRAC_1_SQRT_TAU * exp_neg_ax2
1059        };
1060
1061        (y, dy)
1062    }
1063}
1064
1065#[rustfmt::skip]
1066#[inline(always)]
1067fn erf_d_internal<V: FloatVectorWithBits<Element = f64>, P: Policy, const C: bool, const O: bool>(x0: V, out_exp_neg_x2: &mut V) -> V {
1068    // Extract the sign bit once. abs(x0) = x0 ^ sign, and sign is reused
1069    // for the final operation in every branch, avoiding a redundant bitand.
1070    let sign = x0.signed_zero();
1071    let mut x = (x0 ^ sign).flush_denormals_p::<P>();
1072
1073    // Past |x| = 1.34e154 (and at infinity) x^2 overflows, the rational below is inf/inf
1074    // and the whole thing is NaN. erfc(27.3) has already underflowed and erf(6) is
1075    // exactly 1, so clamping at 32 changes no finite result. Compare-and-select rather
1076    // than `min` so a NaN input stays NaN on every backend.
1077    if const { P::POLICY.check_overflow } {
1078        let cap: V = thermite::const_splat!(f64: 32.0);
1079        x = x.cmp_gt(cap).select(cap, x);
1080    }
1081
1082    // if ignoring denormals (and not clamping), just multiply x0 by itself to save like one
1083    // cycle, instead of waiting on abs(), otherwise use the denormal-flushed x value
1084    let x2 = if const { matches!(P::POLICY.denormal_behavior, DenormalBehavior::Ignore) && !P::POLICY.check_overflow } {
1085        x0 * x0
1086    } else {
1087        x * x
1088    };
1089
1090    // LLVM will still start on exp and interleave it with the below operations.
1091    let e = (-x2).exp_p::<P>();
1092
1093    // `x * x` rounds once, and the exp turns that relative `eps` into a relative `x^2 eps`
1094    // of the result: 47 ulp at x = 14, 237 at 24, the whole of erfc's tail error (measured
1095    // flat at 2.7 ulp over 0..27 once it is gone, `tests/erfc_tail.rs`).
1096    //
1097    // With a hardware FMA the residual `lo = x*x - x2` is exact and
1098    // `e^{-(x2 + lo)} = e^{-x2}(1 - lo)` to first order, `lo` being below `eps * x2 < 1e-13`
1099    // wherever erfc is representable. Two instructions, at every tier. Not for the emulated
1100    // FMA: the residual of an unfused product is meaningless, and the polyfill is not for
1101    // shipped kernels.
1102    //
1103    // Without one, `Best` takes fdlibm's split instead: `z` is `x` with its low 27 mantissa
1104    // bits cleared, so `z * z` is exact, and `e^{-x^2} = e^{-z^2} e^{-t}` with
1105    // `t = (x - z)(x + z)`, both factors exact by Sterbenz, `t < x^2 2^-25 < 2e-5` over the
1106    // whole range, so four Taylor terms of `e^{-t}` are 1e-20. No division and no FMA. It
1107    // replaced six independent divisions that measured no gain in any band.
1108    let e = if const { matches!(V::HAS_NATIVE_FMA, thermite::tribool::True) } {
1109        let lo = x.mul_sub(x, x2);
1110        lo.nmul_add(e, e)
1111    } else if const { P::POLICY.precision.ge(PrecisionPolicy::Best) } {
1112        let z = V::from_bits(x.into_bits::<V::Bits>() & thermite::const_splat!(u64: 0xffff_ffff_f800_0000));
1113        let t = (x - z) * (x + z);
1114        let e = (-(z * z)).exp_p::<P>();
1115        e * t.poly_n_p::<P, _>(&[1.0, -1.0, 0.5, -1.0 / 6.0])
1116    } else {
1117        e
1118    };
1119
1120    let a0: V = thermite::const_splat!(f64: 0.56418958354775629);
1121    let a1 = x + thermite::const_splat!(f64: 2.06955023132914151);
1122
1123    let b0 = x2 + x.mul_adde(thermite::const_splat!(f64: 2.71078540045147805), thermite::const_splat!(f64: 5.80755613130301624));
1124    let b1 = x2 + x.mul_adde(thermite::const_splat!(f64: 3.47954057099518960), thermite::const_splat!(f64: 12.06166887286239555));
1125
1126    let c0 = x2 + x.mul_adde(thermite::const_splat!(f64: 3.47469513777439592), thermite::const_splat!(f64: 12.07402036406381411));
1127    let c1 = x2 + x.mul_adde(thermite::const_splat!(f64: 3.72068443960225092), thermite::const_splat!(f64: 8.44319781003968454));
1128
1129    let d0 = x2 + x.mul_adde(thermite::const_splat!(f64: 4.00561509202259545), thermite::const_splat!(f64: 9.30596659485887898));
1130    let d1 = x2 + x.mul_adde(thermite::const_splat!(f64: 3.90225704029924078), thermite::const_splat!(f64: 6.36161630953880464));
1131
1132    let e0 = x2 + x.mul_adde(thermite::const_splat!(f64: 5.16722705817812584), thermite::const_splat!(f64: 9.12661617673673262));
1133    let e1 = x2 + x.mul_adde(thermite::const_splat!(f64: 4.03296893109262491), thermite::const_splat!(f64: 5.13578530585681539));
1134
1135    let f0 = x2 + x.mul_adde(thermite::const_splat!(f64: 5.95908795446633271), thermite::const_splat!(f64: 9.19435612886969243));
1136    let f1 = x2 + x.mul_adde(thermite::const_splat!(f64: 4.11240942957450885), thermite::const_splat!(f64: 4.48640329523408675));
1137
1138    // One division at every tier. The six independent divisions `Best` used to take here
1139    // measured identical to this in every band on both lowerings (`tests/erfc_tail.rs`,
1140    // 2026-09-01): the product's rounding is not where erfc's error was.
1141    let m = {
1142        let n = (a0 * b0) * (c0 * d0) * (e0 * f0);
1143        let d = (a1 * b1) * (c1 * d1) * (e1 * f1);
1144        n / d
1145    };
1146
1147    if O {
1148        // write this right before we use e normally, so LLVM can interleave exp with the above
1149        *out_exp_neg_x2 = e;
1150    }
1151
1152    if !C {
1153        let y = e.nmul_adde(m, V::ONE) ^ sign;
1154
1155        // `1 - m e^{-x^2}` carries a fixed absolute error of about an ulp of 1, so below
1156        // |x| ~ 1e-3 it has no relative accuracy left (erf(0) itself came out 2.2e-16). At
1157        // `Best` and above, small arguments take fdlibm's erf(x) = x + x R(x^2)/S(x^2) on
1158        // |x| < 0.84375 (libm s_erf.c coefficients, 2^-59 relative), which is exact at zero
1159        // and odd by construction. The f32 kernel has carried the same arm from `Average`.
1160        if const { P::POLICY.precision.ge(PrecisionPolicy::Best) } {
1161            let small = x.cmp_lt(thermite::const_splat!(f64: 0.84375));
1162
1163            if const { P::POLICY.avoid_branching } || small.any() {
1164                let rs = x2.poly_rational_n_p::<P, _, _>(
1165                    &[
1166                        1.28379167095512558561e-01,
1167                        -3.25042107247001499370e-01,
1168                        -2.84817495755985104766e-02,
1169                        -5.77027029648944159157e-03,
1170                        -2.37630166566501626084e-05,
1171                    ],
1172                    &[
1173                        1.0,
1174                        3.97917223959155352819e-01,
1175                        6.50222499887672944485e-02,
1176                        5.08130628187576562776e-03,
1177                        1.32494738004321644526e-04,
1178                        -3.96022827877536812320e-06,
1179                    ],
1180                );
1181
1182                return small.select(x0.mul_adde(rs, x0), y);
1183            }
1184        }
1185
1186        y
1187    } else if const { matches!(V::HAS_NATIVE_FMA, thermite::tribool::True) } {
1188        // exploit instruction-level parallelism if FMA is available
1189        x0.select_negative(m.nmul_add(e, V::TWO), m * e)
1190    } else {
1191        let y = m * e;
1192
1193        x0.select_negative(V::TWO - y, y)
1194    }
1195}
1196
1197/// Every default applies: `expint` on the real line is what they were written for.
1198impl<V: FloatVectorWithBits<Element = f64>> super::ExpIntDetails<f64, V> for V {}
1199
1200/// Minimax fit of `L(x)/x` as a polynomial in `x^2` on `[0, 2]`, relative error
1201/// `6.5e-17` after rounding (`crates/thermite-special/scripts/langevin_coeffs.py`).
1202const LANGEVIN_SMALL_F64: [f64; 16] = [
1203    0.3333333333333333,
1204    -0.022222222222221866,
1205    0.002116402116394456,
1206    -0.0002116402115749962,
1207    2.1377798863187195e-05,
1208    -2.1644034853512783e-06,
1209    2.1925805178692086e-07,
1210    -2.2212830510921946e-08,
1211    2.2491902490670497e-09,
1212    -2.2699964972216279e-10,
1213    2.258858881319397e-11,
1214    -2.149441406576282e-12,
1215    1.8363701272379474e-13,
1216    -1.27083611669092e-14,
1217    6.073730625469803e-16,
1218    -1.4527886936695518e-17,
1219];
1220
1221/// Minimax fit of `L^-1(y) (1 - y^2) / y` as a polynomial in `y^2` on `[0, 0.85^2]`,
1222/// relative error `1.1e-6`. The inverse's Halley seed below the `1/(1-y)` tail. Deg 8
1223/// rather than f32's deg 4 so that one cubic step (constant < 0.07) lands under f64's u.
1224const LANGEVIN_SEED_F64: [f64; 9] = [
1225    3.0000033409892763,
1226    -1.200575454653041,
1227    -0.08655290656271598,
1228    -0.11229572780500324,
1229    0.9977556478992905,
1230    -2.159977066903404,
1231    2.891255158804264,
1232    -0.7894725857798888,
1233    -0.6113193395162252,
1234];
1235
1236/// The `Worst`/`Medium` forward table: same fit as [`LANGEVIN_SMALL_F64`] at degree 10,
1237/// relative error `1.9e-12` (the Medium tier's tolerance is 1e4 eps), five FMAs cheaper.
1238/// The inverse keeps the full table at every tier, since its step is dominated by the
1239/// exp and the division and its Medium tier is documented as full precision.
1240const LANGEVIN_SMALL_F64_LO: [f64; 11] = [
1241    0.3333333333327056,
1242    -0.02222222218374748,
1243    0.002116401724776868,
1244    -0.00021163864787629512,
1245    2.1374572036431745e-05,
1246    -2.160476424670753e-06,
1247    2.162294603438737e-07,
1248    -2.0672900635769838e-08,
1249    1.7227111614831327e-09,
1250    -1.0503475016377067e-10,
1251    3.3022089667780975e-12,
1252];
1253
1254/// Order dispatch for the modified Bessel entry points. `N` is a const parameter, so the
1255/// `if const` collapses to one arm and the unused table is never built.
1256///
1257/// Orders 0 and 1 are closed forms. Everything above seeds from the order-0 form and walks
1258/// the ratio recurrence down. All three arms are selected at compile time, so a call site
1259/// pays for exactly one.
1260#[inline(always)]
1261fn bessel_i_dispatch<P: Policy, V, const N: i32, const SCALED: bool>(x: V) -> V
1262where
1263    V: thermite::vector::FloatVector<Element = f64> + thermite::math::TranscendentalMathWithPolicy,
1264{
1265    if const { N == 0 } {
1266        generic::bessel::ik::bessel_i0_impl::<P, _, _, _, _, SCALED>(x, &crate::tables::bessel::BESSEL_I0_F64)
1267    } else if const { N.unsigned_abs() == 1 } {
1268        generic::bessel::ik::bessel_i1_impl::<P, _, _, _, _, SCALED>(x, &crate::tables::bessel::BESSEL_I1_F64)
1269    } else {
1270        // Orders past 1 seed from the order-0 closed form and walk the ratio recurrence down.
1271        generic::bessel::ik::bessel_in_impl::<P, f64, _, _, _, _, N, SCALED>(x, &crate::tables::bessel::BESSEL_I0_F64)
1272    }
1273}
1274
1275/// Order dispatch for the modified Bessel functions of the second kind.
1276///
1277/// Orders 0 and 1 are closed forms. Above that the recurrence runs **upward**, which is the
1278/// opposite of the `I` family and is stable for exactly that reason: `K` is the dominant
1279/// solution. Both `K` kernels also need the `I` tables, because their small arms are
1280/// `P(x^2) - ln(x) I_0(x)` and `R(x^2) x + 1/x + ln(x) I_1(x)`.
1281#[inline(always)]
1282fn bessel_k_dispatch<P: Policy, V, const N: i32, const SCALED: bool>(x: V) -> V
1283where
1284    V: thermite::vector::FloatVector<Element = f64> + thermite::math::TranscendentalMathWithPolicy,
1285{
1286    use crate::tables::bessel::{BESSEL_I0_F64, BESSEL_I1_F64, BESSEL_K0_F64, BESSEL_K1_F64};
1287    if const { N == 0 } {
1288        generic::bessel::ik::bessel_k0_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(x, &BESSEL_K0_F64, &BESSEL_I0_F64)
1289    } else if const { N.unsigned_abs() == 1 } {
1290        generic::bessel::ik::bessel_k1_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(x, &BESSEL_K1_F64, &BESSEL_I1_F64)
1291    } else {
1292        // The recurrence takes the two seeds, not the tables. Each closed form infers its
1293        // own array lengths here, at the one place that already names them concretely.
1294        let k0 = generic::bessel::ik::bessel_k0_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(
1295            x,
1296            &BESSEL_K0_F64,
1297            &BESSEL_I0_F64,
1298        );
1299        let k1 = generic::bessel::ik::bessel_k1_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(
1300            x,
1301            &BESSEL_K1_F64,
1302            &BESSEL_I1_F64,
1303        );
1304        generic::bessel::ik::bessel_kn_recur::<f64, _, N>(x, k0, k1).1
1305    }
1306}
1307
1308/// `(I_N, I_N prime)`, sharing the order-`N-1` value the recurrence already produces.
1309///
1310/// `I_N' = I_{N-1} - (N/x) I_N`, and at `N = 0` the second term vanishes because
1311/// `I_{-1} = I_1`, so one formula covers every order, with the `N = 0` case written out to
1312/// keep `0/x` from becoming `0/0` at the origin.
1313///
1314/// Scaled adds one term: `d/dx e^{-|x|}f = e^{-|x|}(f' - sgn(x) f)`.
1315#[inline(always)]
1316fn bessel_i_deriv_dispatch<P: Policy, V, const N: i32, const SCALED: bool>(x: V) -> (V, V)
1317where
1318    V: thermite::vector::FloatVector<Element = f64> + thermite::math::TranscendentalMathWithPolicy,
1319{
1320    let (prev, v) = if const { N == 0 } {
1321        (
1322            bessel_i_dispatch::<P, V, 1, SCALED>(x),
1323            bessel_i_dispatch::<P, V, 0, SCALED>(x),
1324        )
1325    } else if const { N.unsigned_abs() == 1 } {
1326        (
1327            bessel_i_dispatch::<P, V, 0, SCALED>(x),
1328            bessel_i_dispatch::<P, V, 1, SCALED>(x),
1329        )
1330    } else {
1331        generic::bessel::ik::bessel_in_pair_impl::<P, f64, _, _, _, _, N, SCALED>(
1332            x,
1333            &crate::tables::bessel::BESSEL_I0_F64,
1334        )
1335    };
1336
1337    let mut d = if const { N == 0 } {
1338        prev
1339    } else {
1340        prev - v * (V::splat(N.unsigned_abs() as f64) / x)
1341    };
1342    if const { SCALED } {
1343        d -= v.copysign(x);
1344    }
1345    (v, d)
1346}
1347
1348/// `(K_N, K_N prime)`. `K_N' = -K_{N-1} - (N/x) K_N`, and `K_{-1} = K_1`.
1349///
1350/// Scaled subtracts rather than adds, since the scaling runs the other way:
1351/// `d/dx e^{x}f = e^{x}(f' + f)`.
1352#[inline(always)]
1353fn bessel_k_deriv_dispatch<P: Policy, V, const N: i32, const SCALED: bool>(x: V) -> (V, V)
1354where
1355    V: thermite::vector::FloatVector<Element = f64> + thermite::math::TranscendentalMathWithPolicy,
1356{
1357    let (prev, v) = if const { N == 0 } {
1358        (
1359            bessel_k_dispatch::<P, V, 1, SCALED>(x),
1360            bessel_k_dispatch::<P, V, 0, SCALED>(x),
1361        )
1362    } else if const { N.unsigned_abs() == 1 } {
1363        (
1364            bessel_k_dispatch::<P, V, 0, SCALED>(x),
1365            bessel_k_dispatch::<P, V, 1, SCALED>(x),
1366        )
1367    } else {
1368        // The upward recurrence walks THROUGH order N-1 on its way to N, so the pair costs
1369        // nothing beyond returning it.
1370        let k0 = generic::bessel::ik::bessel_k0_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(
1371            x,
1372            &crate::tables::bessel::BESSEL_K0_F64,
1373            &crate::tables::bessel::BESSEL_I0_F64,
1374        );
1375        let k1 = generic::bessel::ik::bessel_k1_impl::<P, f64, _, _, _, _, _, _, _, _, SCALED>(
1376            x,
1377            &crate::tables::bessel::BESSEL_K1_F64,
1378            &crate::tables::bessel::BESSEL_I1_F64,
1379        );
1380        generic::bessel::ik::bessel_kn_recur::<f64, _, N>(x, k0, k1)
1381    };
1382
1383    let mut d = if const { N == 0 } {
1384        -prev
1385    } else {
1386        -prev - v * (V::splat(N.unsigned_abs() as f64) / x)
1387    };
1388    if const { SCALED } {
1389        d += v;
1390    }
1391    (v, d)
1392}
thermite_special/specialized/pd.rs

thermite_special/specialized/
pd.rs
