/*

 * // 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::bits::biased_exponent_f64;

use crate::common::*;

use crate::double_double::DoubleDouble;

use crate::exponents::expf;

use crate::logf;

use crate::logs::LOG2_R;

use crate::pow_tables::EXP2_MID1;

use crate::powf_tables::{LOG2_R_TD, LOG2_R2_DD, POWF_R2};

use crate::rounding::CpuRound;

/// Power function for given value for const context.

/// This is simplified version just to make a good approximation on const context.

pub const fn powf(d: f32, n: f32) -> f32 {

    let value = d.abs();

    let c = expf(n * logf(value));

    if n == 1. {

        return d;

    }

    if d < 0.0 {

        let y = n as i32;

        if y % 2 == 0 { c } else { -c }

    } else {

        c

    }

}

pub(crate) trait PowfBackend {

    fn fmaf(&self, x: f32, y: f32, z: f32) -> f32;

    fn fma(&self, x: f64, y: f64, z: f64) -> f64;

    fn polyeval3(&self, x: f64, a0: f64, a1: f64, a2: f64) -> f64;

    fn integerf(&self, x: f32) -> bool;

    fn odd_integerf(&self, x: f32) -> bool;

    fn round(&self, x: f64) -> f64;

    fn quick_mult(&self, x: DoubleDouble, y: DoubleDouble) -> DoubleDouble;

    fn quick_mult_f64(&self, x: DoubleDouble, y: f64) -> DoubleDouble;

    fn dd_polyeval6(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

    ) -> DoubleDouble;

    fn dd_polyeval10(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

        a6: DoubleDouble,

        a7: DoubleDouble,

        a8: DoubleDouble,

        a9: DoubleDouble,

    ) -> DoubleDouble;

    const HAS_FMA: bool;

    const ERR: u64;

}

pub(crate) struct GenPowfBackend {}

impl PowfBackend for GenPowfBackend {

    #[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 {

        f_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 integerf(&self, x: f32) -> bool {

        is_integerf(x)

    }

    #[inline(always)]

    fn odd_integerf(&self, x: f32) -> bool {

        is_odd_integerf(x)

    }

    #[inline(always)]

    fn round(&self, x: f64) -> f64 {

        x.cpu_round()

    }

    #[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 dd_polyeval6(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

    ) -> DoubleDouble {

        use crate::polyeval::dd_quick_polyeval6;

        dd_quick_polyeval6(x, a0, a1, a2, a3, a4, a5)

    }

    #[inline(always)]

    fn dd_polyeval10(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

        a6: DoubleDouble,

        a7: DoubleDouble,

        a8: DoubleDouble,

        a9: DoubleDouble,

    ) -> DoubleDouble {

        use crate::polyeval::dd_quick_polyeval10;

        dd_quick_polyeval10(x, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9)

    }

    #[cfg(any(

        all(

            any(target_arch = "x86", target_arch = "x86_64"),

            target_feature = "fma"

        ),

        target_arch = "aarch64"

    ))]

    const HAS_FMA: bool = true;

    #[cfg(not(any(

        all(

            any(target_arch = "x86", target_arch = "x86_64"),

            target_feature = "fma"

        ),

        target_arch = "aarch64"

    )))]

    const HAS_FMA: bool = false;

    #[cfg(any(

        all(

            any(target_arch = "x86", target_arch = "x86_64"),

            target_feature = "fma"

        ),

        target_arch = "aarch64"

    ))]

    const ERR: u64 = 64;

    #[cfg(not(any(

        all(

            any(target_arch = "x86", target_arch = "x86_64"),

            target_feature = "fma"

        ),

        target_arch = "aarch64"

    )))]

    const ERR: u64 = 128;

}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

pub(crate) struct FmaPowfBackend {}

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

impl PowfBackend for FmaPowfBackend {

    #[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 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 integerf(&self, x: f32) -> bool {

        x.round_ties_even() == x

    }

    #[inline(always)]

    fn odd_integerf(&self, x: f32) -> bool {

        use crate::common::is_odd_integerf_fast;

        is_odd_integerf_fast(x)

    }

    #[inline(always)]

    fn round(&self, x: f64) -> f64 {

        x.round()

    }

    #[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 dd_polyeval6(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

    ) -> DoubleDouble {

        use crate::polyeval::dd_quick_polyeval6_fma;

        dd_quick_polyeval6_fma(x, a0, a1, a2, a3, a4, a5)

    }

    #[inline(always)]

    fn dd_polyeval10(

        &self,

        x: DoubleDouble,

        a0: DoubleDouble,

        a1: DoubleDouble,

        a2: DoubleDouble,

        a3: DoubleDouble,

        a4: DoubleDouble,

        a5: DoubleDouble,

        a6: DoubleDouble,

        a7: DoubleDouble,

        a8: DoubleDouble,

        a9: DoubleDouble,

    ) -> DoubleDouble {

        use crate::polyeval::dd_quick_polyeval10_fma;

        dd_quick_polyeval10_fma(x, a0, a1, a2, a3, a4, a5, a6, a7, a8, a9)

    }

    const HAS_FMA: bool = true;

    const ERR: u64 = 64;

}

#[inline]

const fn larger_exponent(a: f64, b: f64) -> bool {

    biased_exponent_f64(a) >= biased_exponent_f64(b)

}

// Calculate 2^(y * log2(x)) in double-double precision.

// At this point we can reuse the following values:

//   idx_x: index for extra precision of log2 for the middle part of log2(x).

//   dx: the reduced argument for log2(x)

//   y6: 2^6 * y.

//   lo6_hi: the high part of 2^6 * (y - (hi + mid))

//   exp2_hi_mid: high part of 2^(hi + mid)

#[cold]

#[inline(always)]

fn powf_dd<B: PowfBackend>(

    idx_x: i32,

    dx: f64,

    y6: f64,

    lo6_hi: f64,

    exp2_hi_mid: DoubleDouble,

    backend: &B,

) -> f64 {

    // Perform a second range reduction step:

    //   idx2 = round(2^14 * (dx  + 2^-8)) = round ( dx * 2^14 + 2^6)

    //   dx2 = (1 + dx) * r2 - 1

    // Output range:

    //   -0x1.3ffcp-15 <= dx2 <= 0x1.3e3dp-15

    let idx2 = backend.round(backend.fma(

        dx,

        f64::from_bits(0x40d0000000000000),

        f64::from_bits(0x4050000000000000),

    )) as usize;

    let dx2 = backend.fma(1.0 + dx, f64::from_bits(POWF_R2[idx2]), -1.0); // Exact

    const COEFFS: [(u64, u64); 6] = [

        (0x3c7777d0ffda25e0, 0x3ff71547652b82fe),

        (0xbc6777d101cf0a84, 0xbfe71547652b82fe),

        (0x3c7ce04b5140d867, 0x3fdec709dc3a03fd),

        (0x3c7137b47e635be5, 0xbfd71547652b82fb),

        (0xbc5b5a30b3bdb318, 0x3fd2776c516a92a2),

        (0x3c62d2fbd081e657, 0xbfcec70af1929ca6),

    ];

    let dx_dd = DoubleDouble::new(0.0, dx2);

    let p = backend.dd_polyeval6(

        dx_dd,

        DoubleDouble::from_bit_pair(COEFFS[0]),

        DoubleDouble::from_bit_pair(COEFFS[1]),

        DoubleDouble::from_bit_pair(COEFFS[2]),

        DoubleDouble::from_bit_pair(COEFFS[3]),

        DoubleDouble::from_bit_pair(COEFFS[4]),

        DoubleDouble::from_bit_pair(COEFFS[5]),

    );

    // log2(1 + dx2) ~ dx2 * P(dx2)

    let log2_x_lo = backend.quick_mult_f64(p, dx2);

    // Lower parts of (e_x - log2(r1)) of the first range reduction constant

    let log2_r_td = LOG2_R_TD[idx_x as usize];

    let log2_x_mid = DoubleDouble::new(f64::from_bits(log2_r_td.0), f64::from_bits(log2_r_td.1));

    // -log2(r2) + lower part of (e_x - log2(r1))

    let log2_x_m = DoubleDouble::add(DoubleDouble::from_bit_pair(LOG2_R2_DD[idx2]), log2_x_mid);

    // log2(1 + dx2) - log2(r2) + lower part of (e_x - log2(r1))

    // Since we don't know which one has larger exponent to apply Fast2Sum

    // algorithm, we need to check them before calling double-double addition.

    let log2_x = if larger_exponent(log2_x_m.hi, log2_x_lo.hi) {

        DoubleDouble::add(log2_x_m, log2_x_lo)

    } else {

        DoubleDouble::add(log2_x_lo, log2_x_m)

    };

    let lo6_hi_dd = DoubleDouble::new(0.0, lo6_hi);

    // 2^6 * y * (log2(1 + dx2) - log2(r2) + lower part of (e_x - log2(r1)))

    let prod = backend.quick_mult_f64(log2_x, y6);

    // 2^6 * (y * log2(x) - (hi + mid)) = 2^6 * lo

    let lo6 = if larger_exponent(prod.hi, lo6_hi) {

        DoubleDouble::add(prod, lo6_hi_dd)

    } else {

        DoubleDouble::add(lo6_hi_dd, prod)

    };

    const EXP2_COEFFS: [(u64, u64); 10] = [

        (0x0000000000000000, 0x3ff0000000000000),

        (0x3c1abc9e3b398024, 0x3f862e42fefa39ef),

        (0xbba5e43a5429bddb, 0x3f0ebfbdff82c58f),

        (0xbb2d33162491268f, 0x3e8c6b08d704a0c0),

        (0x3a94fb32d240a14e, 0x3e03b2ab6fba4e77),

        (0x39ee84e916be83e0, 0x3d75d87fe78a6731),

        (0xb989a447bfddc5e6, 0x3ce430912f86bfb8),

        (0xb8e31a55719de47f, 0x3c4ffcbfc588ded9),

        (0xb850ba57164eb36b, 0x3bb62c034beb8339),

        (0xb7b8483eabd9642d, 0x3b1b5251ff97bee1),

    ];

    let pp = backend.dd_polyeval10(

        lo6,

        DoubleDouble::from_bit_pair(EXP2_COEFFS[0]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[1]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[2]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[3]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[4]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[5]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[6]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[7]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[8]),

        DoubleDouble::from_bit_pair(EXP2_COEFFS[9]),

    );

    let rr = backend.quick_mult(exp2_hi_mid, pp);

    // Make sure the sum is normalized:

    let r = DoubleDouble::from_exact_add(rr.hi, rr.lo);

    const FRACTION_MASK: u64 = (1u64 << 52) - 1;

    let mut r_bits = r.hi.to_bits();

    if ((r_bits & 0xfff_ffff) == 0) && (r.lo != 0.0) {

        let hi_sign = r.hi.to_bits() >> 63;

        let lo_sign = r.lo.to_bits() >> 63;

        if hi_sign == lo_sign {

            r_bits = r_bits.wrapping_add(1);

        } else if (r_bits & FRACTION_MASK) > 0 {

            r_bits = r_bits.wrapping_sub(1);

        }

    }

    f64::from_bits(r_bits)

}

Unexecuted instantiation: pxfm::powf::powf_dd::<pxfm::powf::FmaPowfBackend>

Unexecuted instantiation: pxfm::powf::powf_dd::<pxfm::powf::GenPowfBackend>

#[inline(always)]

fn powf_gen<B: PowfBackend>(x: f32, y: f32, backend: B) -> f32 {

    let mut x_u = x.to_bits();

    let x_abs = x_u & 0x7fff_ffff;

    let mut y = y;

    let y_u = y.to_bits();

    let y_abs = y_u & 0x7fff_ffff;

    let mut x = x;

    if ((y_abs & 0x0007_ffff) == 0) || (y_abs > 0x4f170000) {

        // y is signaling NaN

        if x.is_nan() || y.is_nan() {

            if y.abs() == 0. {

                return 1.;

            }

            if x == 1. {

                return 1.;

            }

            return f32::NAN;

        }

        // Exceptional exponents.

        if y == 0.0 {

            return 1.0;

        }

        match y_abs {

            0x7f80_0000 => {

                if x_abs > 0x7f80_0000 {

                    // pow(NaN, +-Inf) = NaN

                    return x;

                }

                if x_abs == 0x3f80_0000 {

                    // pow(+-1, +-Inf) = 1.0f

                    return 1.0;

                }

                if x == 0.0 && y_u == 0xff80_0000 {

                    // pow(+-0, -Inf) = +inf and raise FE_DIVBYZERO

                    return f32::INFINITY;

                }

                // pow (|x| < 1, -inf) = +inf

                // pow (|x| < 1, +inf) = 0.0f

                // pow (|x| > 1, -inf) = 0.0f

                // pow (|x| > 1, +inf) = +inf

                return if (x_abs < 0x3f80_0000) == (y_u == 0xff80_0000) {

                    f32::INFINITY

                } else {

                    0.

                };

            }

            _ => {

                match y_u {

                    0x3f00_0000 => {

                        // pow(x, 1/2) = sqrt(x)

                        if x == 0.0 || x_u == 0xff80_0000 {

                            // pow(-0, 1/2) = +0

                            // pow(-inf, 1/2) = +inf

                            // Make sure it is correct for FTZ/DAZ.

                            return x * x;

                        }

                        let r = x.sqrt();

                        return if r.to_bits() != 0x8000_0000 { r } else { 0.0 };

                    }

                    0x3f80_0000 => {

                        return x;

                    } // y = 1.0f

                    0x4000_0000 => return x * x, // y = 2.0f

                    _ => {

                        let is_int = backend.integerf(y);

                        if is_int && (y_u > 0x4000_0000) && (y_u <= 0x41c0_0000) {

                            // Check for exact cases when 2 < y < 25 and y is an integer.

                            let mut msb: i32 = if x_abs == 0 {

                                32 - 2

                            } else {

                                x_abs.leading_zeros() as i32

                            };

                            msb = if msb > 8 { msb } else { 8 };

                            let mut lsb: i32 = if x_abs == 0 {

                                0

                            } else {

                                x_abs.trailing_zeros() as i32

                            };

                            lsb = if lsb > 23 { 23 } else { lsb };

                            let extra_bits: i32 = 32 - 2 - lsb - msb;

                            let iter = y as i32;

                            if extra_bits * iter <= 23 + 2 {

                                // The result is either exact or exactly half-way.

                                // But it is exactly representable in double precision.

                                let x_d = x as f64;

                                let mut result = x_d;

                                for _ in 1..iter {

                                    result *= x_d;

                                }

                                return result as f32;

                            }

                        }

                        if y_abs > 0x4f17_0000 {

                            // if y is NaN

                            if y_abs > 0x7f80_0000 {

                                if x_u == 0x3f80_0000 {

                                    // x = 1.0f

                                    // pow(1, NaN) = 1

                                    return 1.0;

                                }

                                // pow(x, NaN) = NaN

                                return y;

                            }

                            // x^y will be overflow / underflow in single precision.  Set y to a

                            // large enough exponent but not too large, so that the computations

                            // won't be overflow in double precision.

                            y = f32::from_bits((y_u & 0x8000_0000).wrapping_add(0x4f800000u32));

                        }

                    }

                }

            }

        }

    }

    const E_BIAS: u32 = (1u32 << (8 - 1u32)) - 1u32;

    let mut ex = -(E_BIAS as i32);

    let mut sign: u64 = 0;

    if ((x_u & 0x801f_ffffu32) == 0) || x_u >= 0x7f80_0000u32 || x_u < 0x0080_0000u32 {

        if x.is_nan() {

            return f32::NAN;

        }

        if x_u == 0x3f80_0000 {

            return 1.;

        }

        let x_is_neg = x.to_bits() > 0x8000_0000;

        if x == 0.0 {

            let out_is_neg = x_is_neg && backend.odd_integerf(f32::from_bits(y_u));

            if y_u > 0x8000_0000u32 {

                // pow(0, negative number) = inf

                return if x_is_neg {

                    f32::NEG_INFINITY

                } else {

                    f32::INFINITY

                };

            }

            // pow(0, positive number) = 0

            return if out_is_neg { -0.0 } else { 0.0 };

        }

        if x_abs == 0x7f80_0000u32 {

            // x = +-Inf

            let out_is_neg = x_is_neg && backend.odd_integerf(f32::from_bits(y_u));

            if y_u >= 0x7fff_ffff {

                return if out_is_neg { -0.0 } else { 0.0 };

            }

            return if out_is_neg {

                f32::NEG_INFINITY

            } else {

                f32::INFINITY

            };

        }

        if x_abs > 0x7f80_0000 {

            // x is NaN.

            // pow (aNaN, 0) is already taken care above.

            return x;

        }

        // Normalize denormal inputs.

        if x_abs < 0x0080_0000u32 {

            ex = ex.wrapping_sub(64);

            x *= f32::from_bits(0x5f800000);

        }

        // x is finite and negative, and y is a finite integer.

        if x.is_sign_negative() {

            if backend.integerf(y) {

                x = -x;

                if backend.odd_integerf(y) {

                    sign = 0x8000_0000_0000_0000u64;

                }

            } else {

                // pow( negative, non-integer ) = NaN

                return f32::NAN;

            }

        }

    }

    // x^y = 2^( y * log2(x) )

    //     = 2^( y * ( e_x + log2(m_x) ) )

    // First we compute log2(x) = e_x + log2(m_x)

    x_u = x.to_bits();

    // Extract exponent field of x.

    ex = ex.wrapping_add((x_u >> 23) as i32);

    let e_x = ex as f64;

    // Use the highest 7 fractional bits of m_x as the index for look up tables.

    let x_mant = x_u & ((1u32 << 23) - 1);

    let idx_x = (x_mant >> (23 - 7)) as i32;

    // Add the hidden bit to the mantissa.

    // 1 <= m_x < 2

    let m_x = f32::from_bits(x_mant | 0x3f800000);

    // Reduced argument for log2(m_x):

    //   dx = r * m_x - 1.

    // The computation is exact, and -2^-8 <= dx < 2^-7.

    // Then m_x = (1 + dx) / r, and

    //   log2(m_x) = log2( (1 + dx) / r )

    //             = log2(1 + dx) - log2(r).

    let dx = if B::HAS_FMA {

        use crate::logs::LOG_REDUCTION_F32;

        backend.fmaf(

            m_x,

            f32::from_bits(LOG_REDUCTION_F32.0[idx_x as usize]),

            -1.0,

        ) as f64 // Exact.

    } else {

        use crate::logs::LOG_RANGE_REDUCTION;

        backend.fma(

            m_x as f64,

            f64::from_bits(LOG_RANGE_REDUCTION[idx_x as usize]),

            -1.0,

        ) // Exact

    };

    // Degree-5 polynomial approximation:

    //   dx * P(dx) ~ log2(1 + dx)

    // Generated by Sollya with:

    // > P = fpminimax(log2(1 + x)/x, 5, [|D...|], [-2^-8, 2^-7]);

    // > dirtyinfnorm(log2(1 + x)/x - P, [-2^-8, 2^-7]);

    //   0x1.653...p-52

    const COEFFS: [u64; 6] = [

        0x3ff71547652b82fe,

        0xbfe71547652b7a07,

        0x3fdec709dc458db1,

        0xbfd715479c2266c9,

        0x3fd2776ae1ddf8f0,

        0xbfce7b2178870157,

    ];

    let dx2 = dx * dx; // Exact

    let c0 = backend.fma(dx, f64::from_bits(COEFFS[1]), f64::from_bits(COEFFS[0]));

    let c1 = backend.fma(dx, f64::from_bits(COEFFS[3]), f64::from_bits(COEFFS[2]));

    let c2 = backend.fma(dx, f64::from_bits(COEFFS[5]), f64::from_bits(COEFFS[4]));

    let p = backend.polyeval3(dx2, c0, c1, c2);

    // s = e_x - log2(r) + dx * P(dx)

    // Approximation errors:

    //   |log2(x) - s| < ulp(e_x) + (bounds on dx) * (error bounds of P(dx))

    //                 = ulp(e_x) + 2^-7 * 2^-51

    //                 < 2^8 * 2^-52 + 2^-7 * 2^-43

    //                 ~ 2^-44 + 2^-50

    let s = backend.fma(dx, p, f64::from_bits(LOG2_R[idx_x as usize]) + e_x);

    // To compute 2^(y * log2(x)), we break the exponent into 3 parts:

    //   y * log(2) = hi + mid + lo, where

    //   hi is an integer

    //   mid * 2^6 is an integer

    //   |lo| <= 2^-7

    // Then:

    //   x^y = 2^(y * log2(x)) = 2^hi * 2^mid * 2^lo,

    // In which 2^mid is obtained from a look-up table of size 2^6 = 64 elements,

    // and 2^lo ~ 1 + lo * P(lo).

    // Thus, we have:

    //   hi + mid = 2^-6 * round( 2^6 * y * log2(x) )

    // If we restrict the output such that |hi| < 150, (hi + mid) uses (8 + 6)

    // bits, hence, if we use double precision to perform

    //   round( 2^6 * y * log2(x))

    // the lo part is bounded by 2^-7 + 2^(-(52 - 14)) = 2^-7 + 2^-38

    // In the following computations:

    //   y6  = 2^6 * y

    //   hm  = 2^6 * (hi + mid) = round(2^6 * y * log2(x)) ~ round(y6 * s)

    //   lo6 = 2^6 * lo = 2^6 * (y - (hi + mid)) = y6 * log2(x) - hm.

    let y6 = (y * f32::from_bits(0x42800000)) as f64; // Exact.

    let hm = backend.round(s * y6);

    // let log2_rr = LOG2_R2_DD[idx_x as usize];

    // // lo6 = 2^6 * lo.

    // let lo6_hi = f_fmla(y6, e_x + f64::from_bits(log2_rr.1), -hm); // Exact

    // // Error bounds:

    // //   | (y*log2(x) - hm * 2^-6 - lo) / y| < err(dx * p) + err(LOG2_R_DD.lo)

    // //                                       < 2^-51 + 2^-75

    // let lo6 = f_fmla(y6, f_fmla(dx, p, f64::from_bits(log2_rr.0)), lo6_hi);

    // lo6 = 2^6 * lo.

    let lo6_hi = backend.fma(y6, e_x + f64::from_bits(LOG2_R_TD[idx_x as usize].2), -hm); // Exact

    // Error bounds:

    //   | (y*log2(x) - hm * 2^-6 - lo) / y| < err(dx * p) + err(LOG2_R_DD.lo)

    //                                       < 2^-51 + 2^-75

    let lo6 = backend.fma(

        y6,

        backend.fma(dx, p, f64::from_bits(LOG2_R_TD[idx_x as usize].1)),

        lo6_hi,

    );

    // |2^(hi + mid) - exp2_hi_mid| <= ulp(exp2_hi_mid) / 2

    // Clamp the exponent part into smaller range that fits double precision.

    // For those exponents that are out of range, the final conversion will round

    // them correctly to inf/max float or 0/min float accordingly.

    let mut hm_i = unsafe { hm.to_int_unchecked::<i64>() };

    hm_i = if hm_i > (1i64 << 15) {

        1 << 15

    } else if hm_i < (-(1i64 << 15)) {

        -(1 << 15)

    } else {

        hm_i

    };

    let idx_y = hm_i & 0x3f;

    // 2^hi

    let exp_hi_i = (hm_i >> 6).wrapping_shl(52);

    // 2^mid

    let exp_mid_i = EXP2_MID1[idx_y as usize].1;

    // (-1)^sign * 2^hi * 2^mid

    // Error <= 2^hi * 2^-53

    let exp2_hi_mid_i = (exp_hi_i.wrapping_add(exp_mid_i as i64) as u64).wrapping_add(sign);

    let exp2_hi_mid = f64::from_bits(exp2_hi_mid_i);

    // Degree-5 polynomial approximation P(lo6) ~ 2^(lo6 / 2^6) = 2^(lo).

    // Generated by Sollya with:

    // > P = fpminimax(2^(x/64), 5, [|1, D...|], [-2^-1, 2^-1]);

    // > dirtyinfnorm(2^(x/64) - P, [-0.5, 0.5]);

    // 0x1.a2b77e618f5c4c176fd11b7659016cde5de83cb72p-60

    const EXP2_COEFFS: [u64; 6] = [

        0x3ff0000000000000,

        0x3f862e42fefa39ef,

        0x3f0ebfbdff82a23a,

        0x3e8c6b08d7076268,

        0x3e03b2ad33f8b48b,

        0x3d75d870c4d84445,

    ];

    let lo6_sqr = lo6 * lo6;

    let d0 = backend.fma(

        lo6,

        f64::from_bits(EXP2_COEFFS[1]),

        f64::from_bits(EXP2_COEFFS[0]),

    );

    let d1 = backend.fma(

        lo6,

        f64::from_bits(EXP2_COEFFS[3]),

        f64::from_bits(EXP2_COEFFS[2]),

    );

    let d2 = backend.fma(

        lo6,

        f64::from_bits(EXP2_COEFFS[5]),

        f64::from_bits(EXP2_COEFFS[4]),

    );

    let pp = backend.polyeval3(lo6_sqr, d0, d1, d2);

    let r = pp * exp2_hi_mid;

    let r_u = r.to_bits();

    let r_upper = f64::from_bits(r_u + B::ERR) as f32;

    let r_lower = f64::from_bits(r_u - B::ERR) as f32;

    if r_upper == r_lower {

        return r_upper;

    }

    // Scale lower part of 2^(hi + mid)

    let exp2_hi_mid_dd = DoubleDouble {

        lo: if idx_y != 0 {

            f64::from_bits((exp_hi_i as u64).wrapping_add(EXP2_MID1[idx_y as usize].0))

        } else {

            0.

        },

        hi: exp2_hi_mid,

    };

    let r_dd = powf_dd(idx_x, dx, y6, lo6_hi, exp2_hi_mid_dd, &backend);

    r_dd as f32

}

Unexecuted instantiation: pxfm::powf::powf_gen::<pxfm::powf::FmaPowfBackend>

Unexecuted instantiation: pxfm::powf::powf_gen::<pxfm::powf::GenPowfBackend>

#[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

#[target_feature(enable = "avx", enable = "fma")]

unsafe fn powf_fma_impl(x: f32, y: f32) -> f32 {

    powf_gen(x, y, FmaPowfBackend {})

}

/// Power function

///

/// Max found ULP 0.5

pub fn f_powf(x: f32, y: f32) -> f32 {

    #[cfg(not(any(target_arch = "x86", target_arch = "x86_64")))]

    {

        powf_gen(x, y, GenPowfBackend {})

    }

    #[cfg(any(target_arch = "x86", target_arch = "x86_64"))]

    {

        use std::sync::OnceLock;

        static EXECUTOR: OnceLock<unsafe fn(f32, 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")

            {

                powf_fma_impl

            } else {

                fn def_powf(x: f32, y: f32) -> f32 {

                    powf_gen(x, y, GenPowfBackend {})

                }

                def_powf

            }

        });

        unsafe { q(x, y) }

    }

}

/// Dirty fast pow

#[inline]

pub fn dirty_powf(d: f32, n: f32) -> f32 {

    use crate::exponents::dirty_exp2f;

    use crate::logs::dirty_log2f;

    let value = d.abs();

    let lg = dirty_log2f(value);

    let c = dirty_exp2f(n * lg);

    if d < 0.0 {

        let y = n as i32;

        if y % 2 == 0 { c } else { -c }

    } else {

        c

    }

}

Unexecuted instantiation: pxfm::powf::dirty_powf

Unexecuted instantiation: pxfm::powf::dirty_powf

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn powf_test() {

        assert!(

            (powf(2f32, 3f32) - 8f32).abs() < 1e-6,

            "Invalid result {}",

            powf(2f32, 3f32)

        );

        assert!(

            (powf(0.5f32, 2f32) - 0.25f32).abs() < 1e-6,

            "Invalid result {}",

            powf(0.5f32, 2f32)

        );

    }

    #[test]

    fn f_powf_test() {

        assert!(

            (f_powf(2f32, 3f32) - 8f32).abs() < 1e-6,

            "Invalid result {}",

            f_powf(2f32, 3f32)

        );

        assert!(

            (f_powf(0.5f32, 2f32) - 0.25f32).abs() < 1e-6,

            "Invalid result {}",

            f_powf(0.5f32, 2f32)

        );

        assert_eq!(f_powf(0.5f32, 1.5432f32), 0.34312353);

        assert_eq!(

            f_powf(f32::INFINITY, 0.00000000000000000000000000000000038518824),

            f32::INFINITY

        );

        assert_eq!(f_powf(f32::NAN, 0.0), 1.);

        assert_eq!(f_powf(1., f32::NAN), 1.);

    }

    #[test]

    fn dirty_powf_test() {

        assert!(

            (dirty_powf(2f32, 3f32) - 8f32).abs() < 1e-6,

            "Invalid result {}",

            dirty_powf(2f32, 3f32)

        );

        assert!(

            (dirty_powf(0.5f32, 2f32) - 0.25f32).abs() < 1e-6,

            "Invalid result {}",

            dirty_powf(0.5f32, 2f32)

        );

    }

}