pxfm/bessel/
bessel_exp.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::{dd_fmla, f_fmla};
30use crate::double_double::DoubleDouble;
31use crate::exponents::{EXP_REDUCE_T0, EXP_REDUCE_T1};
32use crate::round::RoundFinite;
33
34#[inline(always)]
35fn exp_poly(z: f64) -> DoubleDouble {
36    /* The following is a degree-4 polynomial generated by Sollya for exp(x)
37    over [-2^-12.905,2^-12.905]
38    with absolute error < 2^-74.34 (see sollya/Q_1.sollya). */
39    const Q_1: [u64; 5] = [
40        0x3ff0000000000000,
41        0x3ff0000000000000,
42        0x3fe0000000000000,
43        0x3fc5555555997996,
44        0x3fa5555555849d8d,
45    ];
46    let mut q = dd_fmla(f64::from_bits(Q_1[4]), z, f64::from_bits(Q_1[3]));
47    q = dd_fmla(q, z, f64::from_bits(Q_1[2]));
48    let h0 = dd_fmla(q, z, f64::from_bits(Q_1[1]));
49
50    let v1 = DoubleDouble::from_exact_mult(z, h0);
51    DoubleDouble::f64_add(f64::from_bits(Q_1[0]), v1)
52}
53
54#[inline]
55pub(crate) fn i0_exp(r: f64) -> DoubleDouble {
56    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);
57
58    let k = (r * INVLOG2).round_finite();
59
60    const LOG_2E: DoubleDouble = DoubleDouble::new(
61        f64::from_bits(0x3d0718432a1b0e26),
62        f64::from_bits(0x3f262e42ff000000),
63    );
64
65    let zh = f_fmla(LOG_2E.lo, k, f_fmla(-LOG_2E.hi, k, r));
66
67    let bk = unsafe {
68        k.to_int_unchecked::<i64>() // k is already integer, this is just a conversion
69    };
70    let mk = (bk >> 12) + 0x3ff;
71    let i2 = (bk >> 6) & 0x3f;
72    let i1 = bk & 0x3f;
73
74    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i2 as usize]);
75    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
76    let mut de = DoubleDouble::quick_mult(t1, t0);
77    let q = exp_poly(zh);
78    de = DoubleDouble::quick_mult(de, q);
79
80    let mut du = (mk as u64).wrapping_shl(52);
81    du = f64::from_bits(du).to_bits();
82    DoubleDouble::quick_mult_f64(de, f64::from_bits(du))
83}
84
85#[inline(always)]
86fn exp_poly_dd(z: DoubleDouble) -> DoubleDouble {
87    // Generated by Sollya:
88    // d = [-2^-12.905,2^-12.905];
89    // f = exp(x);
90    // w = 1;
91    // p = remez(f, 6, d, w);
92    // pf = fpminimax(f, [|0,1,2,3,4,5,6|], [|1, 107...|], d, absolute, floating, 0, p);
93    // err_p = -log2(dirtyinfnorm(pf*w-f, d));
94    // display = decimal;
95    const Q_1: [(u64, u64); 7] = [
96        (0x0000000000000000, 0x3ff0000000000000),
97        (0x3a20e40000000000, 0x3ff0000000000000),
98        (0x3a04820000000000, 0x3fe0000000000000),
99        (0xbc756423c5338a66, 0x3fc5555555555556),
100        (0xbc5560f74db5556c, 0x3fa5555555555556),
101        (0x3c3648eca89bc6ac, 0x3f8111111144fbee),
102        (0xbbd53d924ae90c8c, 0x3f56c16c16ffeecc),
103    ];
104    let mut p = DoubleDouble::quick_mul_add(
105        z,
106        DoubleDouble::from_bit_pair(Q_1[6]),
107        DoubleDouble::from_bit_pair(Q_1[5]),
108    );
109    p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[4]));
110    p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[3]));
111    p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[2]));
112    p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[1]));
113    DoubleDouble::quick_mul_add_f64(z, p, f64::from_bits(0x3ff0000000000000))
114}
115
116#[cold]
117pub(crate) fn i0_exp_accurate(r: f64) -> DoubleDouble {
118    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);
119
120    let k = (r * INVLOG2).round_finite();
121
122    const L2: DoubleDouble = DoubleDouble::new(
123        f64::from_bits(0x3d0718432a1b0e26),
124        f64::from_bits(0x3f262e42ff000000),
125    );
126    const L2LL: f64 = f64::from_bits(0x3999ff0342542fc3);
127    let dx = f_fmla(-L2.hi, k, r);
128    let dx_dd = DoubleDouble::quick_mult_f64(DoubleDouble::new(L2LL, L2.lo), k);
129    let dz = DoubleDouble::full_add_f64(dx_dd, dx);
130
131    let bk = unsafe {
132        k.to_int_unchecked::<i64>() // k is already integer, this is just a conversion
133    };
134    let mk = (bk >> 12) + 0x3ff;
135    let i2 = (bk >> 6) & 0x3f;
136    let i1 = bk & 0x3f;
137
138    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i2 as usize]);
139    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
140    let mut de = DoubleDouble::quick_mult(t1, t0);
141    let q = exp_poly_dd(dz);
142    de = DoubleDouble::quick_mult(de, q);
143
144    let mut du = (mk as u64).wrapping_shl(52);
145    du = f64::from_bits(du).to_bits();
146    DoubleDouble::quick_mult_f64(de, f64::from_bits(du))
147}
148
149#[cfg(test)]
150mod tests {
151    use super::*;
152    #[test]
153    fn test_i0() {
154        assert_eq!(i0_exp(0.5).to_f64(), 1.6487212707001282);
155        assert_eq!(i0_exp_accurate(0.5).to_f64(), 1.6487212707001282);
156    }
157}
pxfm/bessel/bessel_exp.rs

pxfm/bessel/
bessel_exp.rs
