pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 4/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.
 */
#![allow(clippy::too_many_arguments)]
use crate::common::{dd_fmla, dyad_fmla, f_fmla, f_fmlaf, fmlaf, pow2if, rintfk};
use crate::double_double::DoubleDouble;
use crate::rounding::{CpuFloor, CpuRound};

const L2U_F: f32 = 0.693_145_751_953_125;
const L2L_F: f32 = 1.428_606_765_330_187_045_e-6;
const R_LN2_F: f32 = std::f32::consts::LOG2_E;

/// Exp for given value for const context.
/// This is simplified version just to make a good approximation on const context.
#[inline]
pub const fn expf(d: f32) -> f32 {
    const EXP_POLY_1_S: f32 = 2f32;
    const EXP_POLY_2_S: f32 = 0.16666707f32;
    const EXP_POLY_3_S: f32 = -0.002775669f32;
    let qf = rintfk(d * R_LN2_F);
    let q = qf as i32;
    let r = fmlaf(qf, -L2U_F, d);
    let r = fmlaf(qf, -L2L_F, r);

    let f = r * r;
    // Poly for u = r*(exp(r)+1)/(exp(r)-1)
    let mut u = EXP_POLY_3_S;
    u = fmlaf(u, f, EXP_POLY_2_S);
    u = fmlaf(u, f, EXP_POLY_1_S);
    let u = 1f32 + 2f32 * r / (u - r);
    let i2 = pow2if(q);
    u * i2
    // if d < -87f32 {
    //     r = 0f32;
    // }
    // if d > 88f32 {
    //     r = f32::INFINITY;
    // }
}

// Lookup table for exp(m) with m = -104, ..., 102.
//   -104 = floor(log(single precision's min denormal))
//    103 = ceil(log(single precision's max bessel K(n) that will be used))
// Table is generated with SageMath as follows:
// for r in range(-104, 103):
//     print(double_to_hex(RealField(180)(r).exp()) + ",")
pub(crate) static EXP_M1: [u64; 207] = [
    0x368f1e6b68529e33,
    0x36a525be4e4e601d,
    0x36bcbe0a45f75eb1,
    0x36d3884e838aea68,
    0x36ea8c1f14e2af5d,
    0x37020a717e64a9bd,
    0x3718851d84118908,
    0x3730a9bdfb02d240,
    0x3746a5bea046b42e,
    0x375ec7f3b269efa8,
    0x3774eafb87eab0f2,
    0x378c6e2d05bbc000,
    0x37a35208867c2683,
    0x37ba425b317eeacd,
    0x37d1d8508fa8246a,
    0x37e840fbc08fdc8a,
    0x38007b7112bc1ffe,
    0x381666d0dad2961d,
    0x382e726c3f64d0fe,
    0x3844b0dc07cabf98,
    0x385c1f2daf3b6a46,
    0x38731c5957a47de2,
    0x3889f96445648b9f,
    0x38a1a6baeadb4fd1,
    0x38b7fd974d372e45,
    0x38d04da4d1452919,
    0x38e62891f06b3450,
    0x38fe1dd273aa8a4a,
    0x3914775e0840bfdd,
    0x392bd109d9d94bda,
    0x3942e73f53fba844,
    0x3959b138170d6bfe,
    0x397175af0cf60ec5,
    0x3987baee1bffa80b,
    0x39a02057d1245ceb,
    0x39b5eafffb34ba31,
    0x39cdca23bae16424,
    0x39e43e7fc88b8056,
    0x39fb83bf23a9a9eb,
    0x3a12b2b8dd05b318,
    0x3a2969d47321e4cc,
    0x3a41452b7723aed2,
    0x3a5778fe2497184c,
    0x3a6fe7116182e9cc,
    0x3a85ae191a99585a,
    0x3a9d775d87da854d,
    0x3ab4063f8cc8bb98,
    0x3acb374b315f87c1,
    0x3ae27ec458c65e3c,
    0x3af923372c67a074,
    0x3b11152eaeb73c08,
    0x3b2737c5645114b5,
    0x3b3f8e6c24b5592e,
    0x3b5571db733a9d61,
    0x3b6d257d547e083f,
    0x3b83ce9b9de78f85,
    0x3b9aebabae3a41b5,
    0x3bb24b6031b49bda,
    0x3bc8dd5e1bb09d7e,
    0x3be0e5b73d1ff53d,
    0x3bf6f741de1748ec,
    0x3c0f36bd37f42f3e,
    0x3c2536452ee2f75c,
    0x3c3cd480a1b74820,
    0x3c539792499b1a24,
    0x3c6aa0de4bf35b38,
    0x3c82188ad6ae3303,
    0x3c9898471fca6055,
    0x3cb0b6c3afdde064,
    0x3cc6b7719a59f0e0,
    0x3cdee001eed62aa0,
    0x3cf4fb547c775da8,
    0x3d0c8464f7616468,
    0x3d236121e24d3bba,
    0x3d3a56e0c2ac7f75,
    0x3d51e642baeb84a0,
    0x3d6853f01d6d53ba,
    0x3d80885298767e9a,
    0x3d967852a7007e42,
    0x3dae8a37a45fc32e,
    0x3dc4c1078fe9228a,
    0x3ddc3527e433fab1,
    0x3df32b48bf117da2,
    0x3e0a0db0d0ddb3ec,
    0x3e21b48655f37267,
    0x3e381056ff2c5772,
    0x3e505a628c699fa1,
    0x3e6639e3175a689d,
    0x3e7e355bbaee85cb,
    0x3e94875ca227ec38,
    0x3eabe6c6fdb01612,
    0x3ec2f6053b981d98,
    0x3ed9c54c3b43bc8b,
    0x3ef18354238f6764,
    0x3f07cd79b5647c9b,
    0x3f202cf22526545a,
    0x3f35fc21041027ad,
    0x3f4de16b9c24a98f,
    0x3f644e51f113d4d6,
    0x3f7b993fe00d5376,
    0x3f92c155b8213cf4,
    0x3fa97db0ccceb0af,
    0x3fc152aaa3bf81cc,
    0x3fd78b56362cef38,
    0x3ff0000000000000,
    0x4005bf0a8b145769,
    0x401d8e64b8d4ddae,
    0x403415e5bf6fb106,
    0x404b4c902e273a58,
    0x40628d389970338f,
    0x407936dc5690c08f,
    0x409122885aaeddaa,
    0x40a749ea7d470c6e,
    0x40bfa7157c470f82,
    0x40d5829dcf950560,
    0x40ed3c4488ee4f7f,
    0x4103de1654d37c9a,
    0x411b00b5916ac955,
    0x413259ac48bf05d7,
    0x4148f0ccafad2a87,
    0x4160f2ebd0a80020,
    0x417709348c0ea4f9,
    0x418f4f22091940bd,
    0x41a546d8f9ed26e1,
    0x41bceb088b68e804,
    0x41d3a6e1fd9eecfd,
    0x41eab5adb9c43600,
    0x420226af33b1fdc1,
    0x4218ab7fb5475fb7,
    0x4230c3d3920962c9,
    0x4246c932696a6b5d,
    0x425ef822f7f6731d,
    0x42750bba3796379a,
    0x428c9aae4631c056,
    0x42a370470aec28ed,
    0x42ba6b765d8cdf6d,
    0x42d1f43fcc4b662c,
    0x42e866f34a725782,
    0x4300953e2f3a1ef7,
    0x431689e221bc8d5b,
    0x432ea215a1d20d76,
    0x4344d13fbb1a001a,
    0x435c4b334617cc67,
    0x43733a43d282a519,
    0x438a220d397972eb,
    0x43a1c25c88df6862,
    0x43b8232558201159,
    0x43d0672a3c9eb871,
    0x43e64b41c6d37832,
    0x43fe4cf766fe49be,
    0x44149767bc0483e3,
    0x442bfc951eb8bb76,
    0x444304d6aeca254b,
    0x4459d97010884251,
    0x44719103e4080b45,
    0x4487e013cd114461,
    0x44a03996528e074c,
    0x44b60d4f6fdac731,
    0x44cdf8c5af17ba3b,
    0x44e45e3076d61699,
    0x44fbaed16a6e0da7,
    0x4512cffdfebde1a1,
    0x4529919cabefcb69,
    0x454160345c9953e3,
    0x45579dbc9dc53c66,
    0x45700c810d464097,
    0x4585d009394c5c27,
    0x459da57de8f107a8,
    0x45b425982cf597cd,
    0x45cb61e5ca3a5e31,
    0x45e29bb825dfcf87,
    0x45f94a90db0d6fe2,
    0x46112fec759586fd,
    0x46275c1dc469e3af,
    0x463fbfd219c43b04,
    0x4655936d44e1a146,
    0x466d531d8a7ee79c,
    0x4683ed9d24a2d51b,
    0x469b15cfe5b6e17b,
    0x46b268038c2c0e00,
    0x46c9044a73545d48,
    0x46e1002ab6218b38,
    0x46f71b3540cbf921,
    0x470f6799ea9c414a,
    0x47255779b984f3eb,
    0x473d01a210c44aa4,
    0x4753b63da8e91210,
    0x476aca8d6b0116b8,
    0x478234de9e0c74e9,
    0x4798bec7503ca477,
    0x47b0d0eda9796b90,
    0x47c6db0118477245,
    0x47df1056dc7bf22d,
    0x47f51c2cc3433801,
    0x480cb108ffbec164,
    0x48237f780991b584,
    0x483a801c0ea8ac4d,
    0x48520247cc4c46c1,
    0x48687a0553328015,
    0x4880a233dee4f9bb,
    0x48969b7f55b808ba,
    0x48aeba064644060a,
    0x48c4e184933d9364,
    0x48dc614fe2531841,
    0x48f3494a9b171bf5,
    0x490a36798b9d969b,
    0x4921d03d8c0c04af,
];

// Lookup table for exp(m * 2^(-7)) with m = 0, ..., 127.
// Table is generated with Sollya as follows:
// > display = hexadecimal;
// > for i from 0 to 127 do { D(exp(i / 128)); };
pub(crate) static EXP_M2: [u64; 128] = [
    0x3ff0000000000000,
    0x3ff0202015600446,
    0x3ff04080ab55de39,
    0x3ff06122436410dd,
    0x3ff08205601127ed,
    0x3ff0a32a84e9c1f6,
    0x3ff0c49236829e8c,
    0x3ff0e63cfa7ab09d,
    0x3ff1082b577d34ed,
    0x3ff12a5dd543ccc5,
    0x3ff14cd4fc989cd6,
    0x3ff16f9157587069,
    0x3ff192937074e0cd,
    0x3ff1b5dbd3f68122,
    0x3ff1d96b0eff0e79,
    0x3ff1fd41afcba45e,
    0x3ff2216045b6f5cd,
    0x3ff245c7613b8a9b,
    0x3ff26a7793f60164,
    0x3ff28f7170a755fd,
    0x3ff2b4b58b372c79,
    0x3ff2da4478b620c7,
    0x3ff3001ecf601af7,
    0x3ff32645269ea829,
    0x3ff34cb8170b5835,
    0x3ff373783a722012,
    0x3ff39a862bd3c106,
    0x3ff3c1e2876834aa,
    0x3ff3e98deaa11dcc,
    0x3ff41188f42c3e32,
    0x3ff439d443f5f159,
    0x3ff462707b2bac21,
    0x3ff48b5e3c3e8186,
    0x3ff4b49e2ae5ac67,
    0x3ff4de30ec211e60,
    0x3ff50817263c13cd,
    0x3ff5325180cfacf7,
    0x3ff55ce0a4c58c7c,
    0x3ff587c53c5a7af0,
    0x3ff5b2fff3210fd9,
    0x3ff5de9176045ff5,
    0x3ff60a7a734ab0e8,
    0x3ff636bb9a983258,
    0x3ff663559cf1bc7c,
    0x3ff690492cbf9433,
    0x3ff6bd96fdd034a2,
    0x3ff6eb3fc55b1e76,
    0x3ff719443a03acb9,
    0x3ff747a513dbef6a,
    0x3ff776630c678bc1,
    0x3ff7a57ede9ea23e,
    0x3ff7d4f946f0ba8d,
    0x3ff804d30347b546,
    0x3ff8350cd30ac390,
    0x3ff865a7772164c5,
    0x3ff896a3b1f66a0e,
    0x3ff8c802477b0010,
    0x3ff8f9c3fd29beaf,
    0x3ff92be99a09bf00,
    0x3ff95e73e6b1b75e,
    0x3ff99163ad4b1dcc,
    0x3ff9c4b9b995509b,
    0x3ff9f876d8e8c566,
    0x3ffa2c9bda3a3e78,
    0x3ffa61298e1e069c,
    0x3ffa9620c6cb3374,
    0x3ffacb82581eee54,
    0x3ffb014f179fc3b8,
    0x3ffb3787dc80f95f,
    0x3ffb6e2d7fa5eb18,
    0x3ffba540dba56e56,
    0x3ffbdcc2cccd3c85,
    0x3ffc14b431256446,
    0x3ffc4d15e873c193,
    0x3ffc85e8d43f7cd0,
    0x3ffcbf2dd7d490f2,
    0x3ffcf8e5d84758a9,
    0x3ffd3311bc7822b4,
    0x3ffd6db26d16cd67,
    0x3ffda8c8d4a66969,
    0x3ffde455df80e3c0,
    0x3ffe205a7bdab73e,
    0x3ffe5cd799c6a54e,
    0x3ffe99ce2b397649,
    0x3ffed73f240dc142,
    0x3fff152b7a07bb76,
    0x3fff539424d90f5e,
    0x3fff927a1e24bb76,
    0x3fffd1de6182f8c9,
    0x400008e0f64294ab,
    0x40002912df5ce72a,
    0x400049856cd84339,
    0x40006a39207f0a09,
    0x40008b2e7d2035cf,
    0x4000ac6606916501,
    0x4000cde041b0e9ae,
    0x4000ef9db467dcf8,
    0x4001119ee5ac36b6,
    0x400133e45d82e952,
    0x4001566ea50201d7,
    0x4001793e4652cc50,
    0x40019c53ccb3fc6b,
    0x4001bfafc47bda73,
    0x4001e352bb1a74ad,
    0x4002073d3f1bd518,
    0x40022b6fe02a3b9c,
    0x40024feb2f105cb8,
    0x400274afbdbba4a6,
    0x400299be1f3e7f1c,
    0x4002bf16e7d2a38c,
    0x4002e4baacdb6614,
    0x40030aaa04e80d05,
    0x400330e587b62b28,
    0x4003576dce33fead,
    0x40037e437282d4ee,
    0x4003a5670ff972ed,
    0x4003ccd9432682b4,
    0x4003f49aa9d30590,
    0x40041cabe304cb34,
    0x4004450d8f00edd4,
    0x40046dc04f4e5338,
    0x400496c4c6b832da,
    0x4004c01b9950a111,
    0x4004e9c56c731f5d,
    0x400513c2e6c731d7,
    0x40053e14b042f9ca,
    0x400568bb722dd593,
    0x400593b7d72305bb,
];

pub(crate) trait ExpfBackend {
    fn fmaf(&self, x: f32, y: f32, z: f32) -> f32;
    fn fma(&self, x: f64, y: f64, z: f64) -> f64;
    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64;
    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64;
    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64;
    fn polyeval5(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64) -> f64;
    fn polyeval6(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64, a5: f64) -> f64;
    fn polyeval7(
        &self,
        x: f64,
        a0: f64,
        a1: f64,
        a2: f64,
        a3: f64,
        a4: f64,
        a5: f64,
        a6: f64,
    ) -> f64;
    fn roundf(&self, x: f32) -> f32;
    fn round(&self, x: f64) -> f64;
    fn floor(&self, x: f64) -> f64;
    fn round_ties_even(&self, x: f64) -> f64;
    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble;
    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble;
    fn quick_f64_mult(&self, x: f64, y: DoubleDouble) -> DoubleDouble;
    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble;
}

pub(crate) struct GenericExpfBackend {}

impl ExpfBackend for GenericExpfBackend {
    #[inline(always)]
    fn fmaf(&self, x: f32, y: f32, z: f32) -> f32 {
        f_fmlaf(x, y, z)
    }

    #[inline(always)]
    fn fma(&self, x: f64, y: f64, z: f64) -> f64 {
        use crate::common::f_fmla;
        f_fmla(x, y, z)
    }
    #[inline(always)]
    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64 {
        dd_fmla(x, y, z)
    }
    #[inline(always)]
    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64 {
        dyad_fmla(x, y, z)
    }

    #[inline(always)]
    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64 {
        use crate::polyeval::f_polyeval3;
        f_polyeval3(x, a0, a1, a2)
    }

    #[inline(always)]
    fn polyeval5(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64) -> f64 {
        use crate::polyeval::f_polyeval5;
        f_polyeval5(x, a0, a1, a2, a3, a4)
    }

    #[inline(always)]
    fn polyeval6(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64, a5: f64) -> f64 {
        use crate::polyeval::f_polyeval6;
        f_polyeval6(x, a0, a1, a2, a3, a4, a5)
    }

    #[inline(always)]
    fn polyeval7(
        &self,
        x: f64,
        a0: f64,
        a1: f64,
        a2: f64,
        a3: f64,
        a4: f64,
        a5: f64,
        a6: f64,
    ) -> f64 {
        use crate::polyeval::f_polyeval7;
        f_polyeval7(x, a0, a1, a2, a3, a4, a5, a6)
    }

    #[inline(always)]
    fn roundf(&self, x: f32) -> f32 {
        x.cpu_round()
    }

    #[inline(always)]
    fn round(&self, x: f64) -> f64 {
        x.cpu_round()
    }

    #[inline(always)]
    fn floor(&self, x: f64) -> f64 {
        x.cpu_floor()
    }

    #[inline(always)]
    fn round_ties_even(&self, x: f64) -> f64 {
        use crate::rounding::CpuRoundTiesEven;
        x.cpu_round_ties_even()
    }

    #[inline(always)]
    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
        DoubleDouble::quick_mult(x, y)
    }

    #[inline(always)]
    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble {
        DoubleDouble::quick_mult_f64(x, y)
    }

    #[inline(always)]
    fn quick_f64_mult(&self, x: f64, y: DoubleDouble) -> DoubleDouble {
        DoubleDouble::quick_mult_f64(y, x)
    }

    #[inline(always)]
    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble {
        DoubleDouble::from_exact_mult(x, y)
    }
}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
pub(crate) struct FmaBackend {}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
impl ExpfBackend for FmaBackend {
    #[inline(always)]
    fn fmaf(&self, x: f32, y: f32, z: f32) -> f32 {
        f32::mul_add(x, y, z)
    }

    #[inline(always)]
    fn fma(&self, x: f64, y: f64, z: f64) -> f64 {
        f64::mul_add(x, y, z)
    }

    #[inline(always)]
    fn dd_fma(&self, x: f64, y: f64, z: f64) -> f64 {
        f64::mul_add(x, y, z)
    }
    #[inline(always)]
    fn dyad_fma(&self, x: f64, y: f64, z: f64) -> f64 {
        f64::mul_add(x, y, z)
    }

    #[inline(always)]
    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64 {
        use crate::polyeval::d_polyeval3;
        d_polyeval3(x, a0, a1, a2)
    }

    #[inline(always)]
    fn polyeval5(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64) -> f64 {
        use crate::polyeval::d_polyeval5;
        d_polyeval5(x, a0, a1, a2, a3, a4)
    }

    #[inline(always)]
    fn polyeval6(&self, x: f64, a0: f64, a1: f64, a2: f64, a3: f64, a4: f64, a5: f64) -> f64 {
        use crate::polyeval::d_polyeval6;
        d_polyeval6(x, a0, a1, a2, a3, a4, a5)
    }

    #[inline(always)]
    fn polyeval7(
        &self,
        x: f64,
        a0: f64,
        a1: f64,
        a2: f64,
        a3: f64,
        a4: f64,
        a5: f64,
        a6: f64,
    ) -> f64 {
        use crate::polyeval::d_polyeval7;
        d_polyeval7(x, a0, a1, a2, a3, a4, a5, a6)
    }

    #[inline(always)]
    fn roundf(&self, x: f32) -> f32 {
        x.round()
    }

    #[inline(always)]
    fn round(&self, x: f64) -> f64 {
        x.round()
    }

    #[inline(always)]
    fn floor(&self, x: f64) -> f64 {
        x.floor()
    }

    #[inline(always)]
    fn round_ties_even(&self, x: f64) -> f64 {
        x.round_ties_even()
    }

    #[inline(always)]
    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble {
        DoubleDouble::quick_mult_fma(x, y)
    }

    #[inline(always)]
    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble {
        DoubleDouble::quick_mult_f64_fma(x, y)
    }

    #[inline(always)]
    fn quick_f64_mult(&self, x: f64, y: DoubleDouble) -> DoubleDouble {
        DoubleDouble::quick_mult_f64_fma(y, x)
    }

    #[inline(always)]
    fn exact_mult(&self, x: f64, y: f64) -> DoubleDouble {
        DoubleDouble::from_exact_mult_fma(x, y)
    }
}

#[inline(always)]
fn expf_gen_impl<B: ExpfBackend>(x: f32, backend: B) -> f32 {
    let x_u = x.to_bits();
    let x_abs = x_u & 0x7fff_ffffu32;
    if x_abs >= 0x42b2_0000u32 || x_abs <= 0x3280_0000u32 {
        let exp = ((x_u >> 23) & 0xFF) as i32;
        // |x| < 2^-25
        if exp <= 101i32 {
            return 1.0 + x;
        }

        // When x < log(2^-150) or nan
        if x_u >= 0xc2cf_f1b5u32 {
            // exp(-Inf) = 0
            if x.is_infinite() {
                return 0.0;
            }
            // exp(nan) = nan
            if x.is_nan() {
                return x;
            }
            return 0.0;
        }
        // x >= 89 or nan
        if x.is_sign_positive() && (x_u >= 0x42b2_0000) {
            // x is +inf or nan
            return x + f32::INFINITY;
        }
    }

    // For -104 < x < 89, to compute exp(x), we perform the following range
    // reduction: find hi, mid, lo such that:
    //   x = hi + mid + lo, in which
    //     hi is an integer,
    //     mid * 2^7 is an integer
    //     -2^(-8) <= lo < 2^-8.
    // In particular,
    //   hi + mid = round(x * 2^7) * 2^(-7).
    // Then,
    //   exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
    // We store exp(hi) and exp(mid) in the lookup tables EXP_M1 and EXP_M2
    // respectively.  exp(lo) is computed using a degree-4 minimax polynomial
    // generated by Sollya.

    // x_hi = (hi + mid) * 2^7 = round(x * 2^7).
    let kf = backend.roundf(x * 128.);
    // Subtract (hi + mid) from x to get lo.
    let xd = backend.fmaf(kf, -0.0078125 /* - 1/128 */, x) as f64;
    let mut x_hi = unsafe { kf.to_int_unchecked::<i32>() }; // it's already not indeterminate.
    x_hi += 104 << 7;
    // hi = x_hi >> 7
    let exp_hi = f64::from_bits(EXP_M1[(x_hi >> 7) as usize]);
    // mid * 2^7 = x_hi & 0x0000'007fU;
    let exp_mid = f64::from_bits(EXP_M2[(x_hi & 0x7f) as usize]);
    // Degree-4 minimax polynomial generated by Sollya with the following
    // commands:
    // d = [-2^-8, 2^-8];
    // f_exp = expm1(x)/x;
    // Q = fpminimax(f_exp, 3, [|D...|], [-2^-8, 2^-8]);
    let p = backend.polyeval5(
        xd,
        1.,
        f64::from_bits(0x3feffffffffff777),
        f64::from_bits(0x3fe000000000071c),
        f64::from_bits(0x3fc555566668e5e7),
        f64::from_bits(0x3fa55555555ef243),
    );
    (p * exp_hi * exp_mid) as f32
}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
#[target_feature(enable = "avx", enable = "fma")]
unsafe fn expf_fma_impl(x: f32) -> f32 {
    expf_gen_impl(x, FmaBackend {})
}

/// Computes exp
///
/// Max found ULP 0.5
#[inline]
pub fn f_expf(x: f32) -> f32 {
    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]
    {
        expf_gen_impl(x, GenericExpfBackend {})
    }
    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]
    {
        use std::sync::OnceLock;
        static EXECUTOR: OnceLock<unsafe fn(f32) -> f32> = OnceLock::new();
        let q = EXECUTOR.get_or_init(|| {
            if std::arch::is_x86_feature_detected!("avx")
                && std::arch::is_x86_feature_detected!("fma")
            {
                expf_fma_impl
            } else {
                fn def_expf(x: f32) -> f32 {
                    expf_gen_impl(x, GenericExpfBackend {})
                }
                def_expf
            }
        });
        unsafe { q(x) }
    }
}

#[inline]
pub(crate) fn core_expf(x: f32) -> f64 {
    // x_hi = (hi + mid) * 2^7 = round(x * 2^7).
    let kf = (x * 128.).cpu_round();
    // Subtract (hi + mid) from x to get lo.
    let xd = f_fmlaf(kf, -0.0078125 /* - 1/128 */, x) as f64;
    let mut x_hi = unsafe { kf.to_int_unchecked::<i32>() }; // it's already not indeterminate.
    x_hi += 104 << 7;
    // hi = x_hi >> 7
    let exp_hi = f64::from_bits(EXP_M1[(x_hi >> 7) as usize]);
    // mid * 2^7 = x_hi & 0x0000'007fU;
    let exp_mid = f64::from_bits(EXP_M2[(x_hi & 0x7f) as usize]);
    // Degree-4 minimax polynomial generated by Sollya with the following
    // commands:
    // d = [-2^-8, 2^-8];
    // f_exp = expm1(x)/x;
    // Q = fpminimax(f_exp, 3, [|D...|], [-2^-8, 2^-8]);
    use crate::polyeval::f_polyeval5;
    let p = f_polyeval5(
        xd,
        1.,
        f64::from_bits(0x3feffffffffff777),
        f64::from_bits(0x3fe000000000071c),
        f64::from_bits(0x3fc555566668e5e7),
        f64::from_bits(0x3fa55555555ef243),
    );
    p * exp_hi * exp_mid
}

#[inline]
pub(crate) fn core_expdf(x: f64) -> f64 {
    // x_hi = (hi + mid) * 2^7 = round(x * 2^7).
    let kf = (x * 128.).cpu_round();
    // Subtract (hi + mid) from x to get lo.
    let xd = f_fmla(kf, -0.0078125 /* - 1/128 */, x);
    let mut x_hi = unsafe { kf.to_int_unchecked::<i32>() }; // it's already not indeterminate.
    x_hi += 104 << 7;
    // hi = x_hi >> 7
    let exp_hi = f64::from_bits(EXP_M1[(x_hi >> 7) as usize]);
    // mid * 2^7 = x_hi & 0x0000'007fU;
    let exp_mid = f64::from_bits(EXP_M2[(x_hi & 0x7f) as usize]);
    // Degree-4 minimax polynomial generated by Sollya with the following
    // commands:
    // d = [-2^-8, 2^-8];
    // f_exp = expm1(x)/x;
    // Q = fpminimax(f_exp, 3, [|D...|], [-2^-8, 2^-8]);
    use crate::polyeval::f_polyeval5;
    let p = f_polyeval5(
        xd,
        1.,
        f64::from_bits(0x3feffffffffff777),
        f64::from_bits(0x3fe000000000071c),
        f64::from_bits(0x3fc555566668e5e7),
        f64::from_bits(0x3fa55555555ef243),
    );
    p * exp_hi * exp_mid
}

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

    #[test]
    fn expf_test() {
        assert!(
            (expf(0f32) - 1f32).abs() < 1e-6,
            "Invalid result {}",
            expf(0f32)
        );
        assert!(
            (expf(5f32) - 148.4131591025766f32).abs() < 1e-6,
            "Invalid result {}",
            expf(5f32)
        );
    }

    #[test]
    fn f_expf_test() {
        assert_eq!(f_expf(-103.971596), 1e-45);
        assert!(
            (f_expf(0f32) - 1f32).abs() < 1e-6,
            "Invalid result {}",
            f_expf(0f32)
        );
        assert!(
            (f_expf(5f32) - 148.4131591025766f32).abs() < 1e-6,
            "Invalid result {}",
            f_expf(5f32)
        );
        assert_eq!(f_expf(f32::INFINITY), f32::INFINITY);
        assert_eq!(f_expf(f32::NEG_INFINITY), 0.);
        assert!(f_expf(f32::NAN).is_nan());
    }
}