pxfm/exponents/

exp2m1f.rs

/*
 * // Copyright (c) Radzivon Bartoshyk 6/2025. All rights reserved.
 * //
 * // Redistribution and use in source and binary forms, with or without modification,
 * // are permitted provided that the following conditions are met:
 * //
 * // 1.  Redistributions of source code must retain the above copyright notice, this
 * // list of conditions and the following disclaimer.
 * //
 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
 * // this list of conditions and the following disclaimer in the documentation
 * // and/or other materials provided with the distribution.
 * //
 * // 3.  Neither the name of the copyright holder nor the names of its
 * // contributors may be used to endorse or promote products derived from
 * // this software without specific prior written permission.
 * //
 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
use crate::common::f_fmla;
use crate::exponents::exp2f::EXP2F_TABLE;
use crate::polyeval::f_polyeval3;
use crate::round::RoundFinite;

/// Computes 2^x - 1
///
/// Max found ULP 0.5
#[inline]
pub fn f_exp2m1f(x: f32) -> f32 {
    let x_u = x.to_bits();
    let x_abs = x_u & 0x7fff_ffffu32;
    if x_abs >= 0x4300_0000u32 || x_abs <= 0x3d00_0000u32 {
        // |x| <= 2^-5
        if x_abs <= 0x3d00_0000u32 {
            // Minimax polynomial generated by Sollya with:
            // > display = hexadecimal;
            // > fpminimax((2^x - 1)/x, 5, [|D...|], [-2^-5, 2^-5]);
            const C: [u64; 6] = [
                0x3fe62e42fefa39f3,
                0x3fcebfbdff82c57b,
                0x3fac6b08d6f2d7aa,
                0x3f83b2ab6fc92f5d,
                0x3f55d897cfe27125,
                0x3f243090e61e6af1,
            ];
            let xd = x as f64;
            let xsq = xd * xd;
            let c0 = f_fmla(xd, f64::from_bits(C[1]), f64::from_bits(C[0]));
            let c1 = f_fmla(xd, f64::from_bits(C[3]), f64::from_bits(C[2]));
            let c2 = f_fmla(xd, f64::from_bits(C[5]), f64::from_bits(C[4]));
            let p = f_polyeval3(xsq, c0, c1, c2);
            return (p * xd) as f32;
        }

        // x >= 128, or x is nan
        if x.is_sign_positive() {
            // x >= 128 and 2^x - 1 rounds to +inf, or x is +inf or nan
            return x + f32::INFINITY;
        }
    }

    if x <= -25.0 {
        // 2^(-inf) - 1 = -1
        if x.is_infinite() {
            return -1.0;
        }
        // 2^nan - 1 = nan
        if x.is_nan() {
            return x;
        }
        return -1.0;
    }

    // For -25 < x < 128, to compute 2^x, we perform the following range
    // reduction: find hi, mid, lo such that:
    //   x = hi + mid + lo, in which:
    //     hi is an integer,
    //     0 <= mid * 2^5 < 32 is an integer,
    //     -2^(-6) <= lo <= 2^(-6).
    // In particular,
    //   hi + mid = round(x * 2^5) * 2^(-5).
    // Then,
    //   2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
    // 2^mid is stored in the lookup table of 32 elements.
    // 2^lo is computed using a degree-4 minimax polynomial generated by Sollya.
    // We perform 2^hi * 2^mid by simply add hi to the exponent field of 2^mid.

    // kf = (hi + mid) * 2^5 = round(x * 2^5)

    let xd = x as f64;

    let kf = (x * 64.0).round_finite();
    let k = unsafe { kf.to_int_unchecked::<i32>() }; // it's already not indeterminate.
    // dx = lo = x - (hi + mid) = x - kf * 2^(-6)
    let dx = f_fmla(f64::from_bits(0xbf90000000000000), kf as f64, xd);

    const TABLE_BITS: u32 = 6;
    const TABLE_MASK: u64 = (1u64 << TABLE_BITS) - 1;

    // hi = floor(kf * 2^(-5))
    // exp_hi = shift hi to the exponent field of double precision.
    let exp_hi: i64 = ((k >> TABLE_BITS) as i64).wrapping_shl(52);

    // mh = 2^hi * 2^mid
    // mh_bits = bit field of mh
    let mh_bits = (EXP2F_TABLE[((k as u64) & TABLE_MASK) as usize] as i64).wrapping_add(exp_hi);
    let mh = f64::from_bits(mh_bits as u64);

    // Degree-4 polynomial approximating (2^x - 1)/x generated by Sollya with:
    // > P = fpminimax((2^y - 1)/y, 4, [|D...|], [-1/64. 1/64]);
    // see ./notes/exp2f.sollya
    const C: [u64; 5] = [
        0x3fe62e42fefa39ef,
        0x3fcebfbdff8131c4,
        0x3fac6b08d7061695,
        0x3f83b2b1bee74b2a,
        0x3f55d88091198529,
    ];
    let dx_sq = dx * dx;
    let c1 = f_fmla(dx, f64::from_bits(C[0]), 1.0);
    let c2 = f_fmla(dx, f64::from_bits(C[2]), f64::from_bits(C[1]));
    let c3 = f_fmla(dx, f64::from_bits(C[4]), f64::from_bits(C[3]));
    let p = f_polyeval3(dx_sq, c1, c2, c3);
    // 2^x = 2^(hi + mid + lo)
    //     = 2^(hi + mid) * 2^lo
    //     ~ mh * (1 + lo * P(lo))
    //     = mh + (mh*lo) * P(lo)
    f_fmla(p, mh, -1.) as f32
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_exp2m1f() {
        assert_eq!(f_exp2m1f(0.432423), 0.34949815);
        assert_eq!(f_exp2m1f(-4.), -0.9375);
        assert_eq!(f_exp2m1f(5.43122), 42.14795);
        assert_eq!(f_exp2m1f(4.), 15.0);
        assert_eq!(f_exp2m1f(3.), 7.);
        assert_eq!(f_exp2m1f(0.1), 0.07177346);
        assert_eq!(f_exp2m1f(0.0543432432), 0.038386293);
        assert!(f_exp2m1f(f32::NAN).is_nan());
        assert_eq!(f_exp2m1f(f32::INFINITY), f32::INFINITY);
        assert_eq!(f_exp2m1f(f32::NEG_INFINITY), -1.0);
    }
}