pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 7/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::dd_fmla;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::exponents::{EXP_REDUCE_T0, EXP_REDUCE_T1, ldexp};
use crate::exponents::{EXPM1_T0, EXPM1_T1};
use crate::polyeval::f_polyeval6;
use crate::pow_tables::{EXP_T1_2_DYADIC, EXP_T2_2_DYADIC, POW_INVERSE, POW_LOG_INV};
use crate::rounding::CpuRound;
use crate::rounding::CpuRoundTiesEven;

#[inline(always)]
pub(crate) fn log_poly_1(z: f64) -> DoubleDouble {
    /* The following is a degree-8 polynomial generated by Sollya for
    log(1+x)-x+x^2/2 over [-0.0040283203125,0.0040283203125]
    with absolute error < 2^-81.63
    and relative error < 2^-72.423 (see sollya/P_1.sollya).
    The relative error is for x - x^2/2 + P(x) with respect to log(1+x). */
    const P_1: [u64; 6] = [
        0x3fd5555555555558,
        0xbfd0000000000003,
        0x3fc999999981f535,
        0xbfc55555553d1eb4,
        0x3fc2494526fd4a06,
        0xbfc0001f0c80e8ce,
    ];
    let w = DoubleDouble::from_exact_mult(z, z);
    let t = dd_fmla(f64::from_bits(P_1[5]), z, f64::from_bits(P_1[4]));
    let mut u = dd_fmla(f64::from_bits(P_1[3]), z, f64::from_bits(P_1[2]));
    let mut v = dd_fmla(f64::from_bits(P_1[1]), z, f64::from_bits(P_1[0]));
    u = dd_fmla(t, w.hi, u);
    v = dd_fmla(u, w.hi, v);
    u = v * w.hi;
    DoubleDouble::new(dd_fmla(u, z, -0.5 * w.lo), -0.5 * w.hi)
}

/* Given 2^-1074 <= x <= 0x1.fffffffffffffp+1023, this routine puts in h+l
   an approximation of log(x) such that |l| < 2^-23.89*|h| and

   | h + l - log(x) | <= elog * |log x|

   with elog = 2^-73.527  if x < 1/sqrt(2) or sqrt(2) < x,
   and  elog = 2^-67.0544 if 1/sqrt(2) < x < sqrt(2)
   (note that x cannot equal 1/sqrt(2) nor sqrt(2)).
*/
#[inline]
pub(crate) fn pow_log_1(x: f64) -> (DoubleDouble, bool) {
    /* for 181 <= i <= 362, r[i] = _INVERSE[i-181] is a 9-bit approximation of
    1/x[i], where i*2^-8 <= x[i] < (i+1)*2^-8.
    More precisely r[i] is a 9-bit value such that r[i]*y-1 is representable
    exactly on 53 bits for for any y, i*2^-8 <= y < (i+1)*2^-8.
    Moreover |r[i]*y-1| < 0.0040283203125.
    Table generated with the accompanying pow.sage file,
    with l=inverse_centered(k=8,prec=9,maxbits=53,verbose=false) */
    let x_u = x.to_bits();
    let mut m = x_u & 0xfffffffffffff;
    let mut e: i64 = ((x_u >> 52) & 0x7ff) as i64;

    let t;
    if e != 0 {
        t = m | (0x3ffu64 << 52);
        m = m.wrapping_add(1u64 << 52);
        e -= 0x3ff;
    } else {
        /* x is a subnormal double  */
        let k = m.leading_zeros() - 11;

        e = -0x3fei64 - k as i64;
        m = m.wrapping_shl(k);
        t = m | (0x3ffu64 << 52);
    }

    /* now |x| = 2^_e*_t = 2^(_e-52)*m with 1 <= _t < 2,
    and 2^52 <= _m < 2^53 */

    //   log(x) = log(t) + E · log(2)
    let mut t = f64::from_bits(t);

    // If m > sqrt(2) we divide it by 2 so ensure 1/sqrt(2) < t < sqrt(2)
    let c: usize = (m >= 0x16a09e667f3bcd) as usize;
    static CY: [f64; 2] = [1.0, 0.5];
    static CM: [u64; 2] = [44, 45];

    e = e.wrapping_add(c as i64);
    let be = e;
    let i = m >> CM[c]; /* i/2^8 <= t < (i+1)/2^8 */
    /* when c=1, we have 0x16a09e667f3bcd <= m < 2^53, thus 90 <= i <= 127;
    when c=0, we have 2^52 <= m < 0x16a09e667f3bcd, thus 128 <= i <= 181 */
    t *= CY[c];
    /* now 0x1.6a09e667f3bcdp-1 <= t < 0x1.6a09e667f3bcdp+0,
    and log(x) = E * log(2) + log(t) */

    let r = f64::from_bits(POW_INVERSE[(i - 181) as usize]);
    let l1 = f64::from_bits(POW_LOG_INV[(i - 181) as usize].1);
    let l2 = f64::from_bits(POW_LOG_INV[(i - 181) as usize].0);

    let z = dd_fmla(r, t, -1.0);

    const LOG2_DD: DoubleDouble = DoubleDouble::new(
        f64::from_bits(0x3d2ef35793c76730),
        f64::from_bits(0x3fe62e42fefa3800),
    );

    let th = dd_fmla(be as f64, LOG2_DD.hi, l1);
    let tl = dd_fmla(be as f64, LOG2_DD.lo, l2);

    let mut v = DoubleDouble::f64_add(th, DoubleDouble::new(tl, z));
    let p = log_poly_1(z);
    v = DoubleDouble::f64_add(v.hi, DoubleDouble::new(v.lo + p.lo, p.hi));

    if e == 0 && v.lo.abs() > (v.hi.abs()) * f64::from_bits(0x3e70000000000000) {
        v = DoubleDouble::from_exact_add(v.hi, v.lo);
        return (v, true);
    }

    (v, false)
}

/* Given z such that |z| < 2^-12.905,
   this routine puts in qh+ql an approximation of exp(z) such that

   | (qh+ql) / exp(z) - 1 | < 2^-64.902632

   and |ql| <= 2^-51.999.
*/
#[inline(always)]
fn exp_poly_1(z: f64) -> DoubleDouble {
    /* The following is a degree-4 polynomial generated by Sollya for exp(x)
    over [-2^-12.905,2^-12.905]
    with absolute error < 2^-74.34 (see sollya/Q_1.sollya). */
    const Q_1: [u64; 5] = [
        0x3ff0000000000000,
        0x3ff0000000000000,
        0x3fe0000000000000,
        0x3fc5555555997996,
        0x3fa5555555849d8d,
    ];
    let mut q = dd_fmla(f64::from_bits(Q_1[4]), z, f64::from_bits(Q_1[3]));
    q = dd_fmla(q, z, f64::from_bits(Q_1[2]));
    let h0 = dd_fmla(q, z, f64::from_bits(Q_1[1]));

    let v1 = DoubleDouble::from_exact_mult(z, h0);
    DoubleDouble::f64_add(f64::from_bits(Q_1[0]), v1)
}

/* Given z such that |z| < 2^-12.905,
   this routine puts in qh+ql an approximation of exp(z) such that

   | (qh+ql) / exp(z) - 1 | < 2^-64.902632

   and |ql| <= 2^-51.999.
*/

// #[inline(always)]
// fn exp_poly_dd(z: DoubleDouble) -> DoubleDouble {
//     /* The following is a degree-4 polynomial generated by Sollya for exp(x)
//     over [-2^-12.905,2^-12.905] */
//     const Q_1: [(u64, u64); 7] = [
//         (0x0000000000000000, 0x3ff0000000000000),
//         (0x3a20e40000000000, 0x3ff0000000000000),
//         (0x3a04820000000000, 0x3fe0000000000000),
//         (0xbc756423c5338a66, 0x3fc5555555555556),
//         (0xbc5560f74db5556c, 0x3fa5555555555556),
//         (0x3c3648eca89bc6ac, 0x3f8111111144fbee),
//         (0xbbd53d924ae90c8c, 0x3f56c16c16ffeecc),
//     ];
//     let mut p = DoubleDouble::mult(z, DoubleDouble::from_bit_pair(Q_1[6]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[5]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[4]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[3]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[2]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[1]));
//     p = DoubleDouble::quick_mul_add(z, p, DoubleDouble::from_bit_pair(Q_1[0]));
//     p
// }

#[inline]
pub(crate) fn pow_exp_1(r: DoubleDouble, s: f64) -> DoubleDouble {
    const RHO0: f64 = f64::from_bits(0xc0874910ee4e8a27);
    // #define RHO1 -0x1.577453f1799a6p+9
    /* We increase the initial value of RHO1 to avoid spurious underflow in
    the result value el. However, it is not possible to obtain a lower
    bound on |el| from the input value rh, thus this modified value of RHO1
    is obtained experimentally. */
    const RHO1: f64 = f64::from_bits(0xc08483b8cca421af);
    const RHO2: f64 = f64::from_bits(0x40862e42e709a95b);
    const RHO3: f64 = f64::from_bits(0x40862e4316ea5df9);

    // use !(rh <= RHO2) instead of rh < RHO2 to catch rh = NaN too
    #[allow(clippy::neg_cmp_op_on_partial_ord)]
    if !(r.hi <= RHO2) {
        return if r.hi > RHO3 {
            /* If rh > RHO3, we are sure there is overflow,
            For s=1 we return eh = el = DBL_MAX, which yields
            res_min = res_max = +Inf for rounding up or to nearest,
            and res_min = res_max = DBL_MAX for rounding down or toward zero,
            which will yield the correct rounding.
            For s=-1 we return eh = el = -DBL_MAX, which similarly gives
            res_min = res_max = -Inf or res_min = res_max = -DBL_MAX,
            which is the correct rounding. */
            DoubleDouble::new(
                f64::from_bits(0x7fefffffffffffff) * s,
                f64::from_bits(0x7fefffffffffffff) * s,
            )
        } else {
            /* If RHO2 < rh <= RHO3, we are in the intermediate region
            where there might be overflow or not, thus we set eh = el = NaN,
            which will set res_min = res_max = NaN, the comparison
            res_min == res_max will fail: we defer to the 2nd phase. */
            DoubleDouble::new(f64::NAN, f64::NAN)
        };
    }

    if r.hi < RHO1 {
        return if r.hi < RHO0 {
            /* For s=1, we have eh=el=+0 except for rounding up,
               thus res_min=+0 or -0, res_max=+0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=1 and rounding up, we have eh=+0, el=2^-1074,
               thus res_min = res_max = 2^-1074, which is the expected result too.
               For s=-1, we have eh=el=-0 except for rounding down,
               thus res_min=-0 or +0, res_max=-0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=-1 and rounding down, we have eh=-0, el=-2^-1074,
               thus res_min = res_max = -2^-1074, which is the expected result too.
            */
            DoubleDouble::new(f64::from_bits(0x0000000000000001) * (0.5 * s), 0.0 * s)
        } else {
            /* RHO0 <= rh < RHO1 or s < 0: we defer to the 2nd phase */
            DoubleDouble::new(f64::NAN, f64::NAN)
        };
    }
    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);

    let k = (r.hi * INVLOG2).cpu_round();

    const LOG2H: f64 = f64::from_bits(0x3f262e42fefa39ef);
    const LOG2L: f64 = f64::from_bits(0x3bbabc9e3b39803f);

    let zh = dd_fmla(LOG2H, -k, r.hi);
    let zl = dd_fmla(LOG2L, -k, r.lo);

    let bk = unsafe { k.to_int_unchecked::<i64>() }; /* Note: k is an integer, this is just a conversion. */
    let mk = (bk >> 12) + 0x3ff;
    let i2 = (bk >> 6) & 0x3f;
    let i1 = bk & 0x3f;

    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i2 as usize]);
    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
    let mut de = DoubleDouble::quick_mult(t1, t0);
    let q = exp_poly_1(zh + zl);
    de = DoubleDouble::quick_mult(de, q);
    /* we should have 1 < M < 2047 here, since we filtered out
    potential underflow/overflow cases at the beginning of this function */

    let mut du = (mk as u64).wrapping_shl(52);
    du = (f64::from_bits(du) * s).to_bits();
    de.hi *= f64::from_bits(du);
    de.lo *= f64::from_bits(du);
    de
}

#[inline]
pub(crate) fn exp_dd_fast(r: DoubleDouble) -> DoubleDouble {
    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);

    let k = (r.hi * INVLOG2).cpu_round();

    const LOG2H: f64 = f64::from_bits(0x3f262e42fefa39ef);
    const LOG2L: f64 = f64::from_bits(0x3bbabc9e3b39803f);

    let mut z = DoubleDouble::mul_f64_add(DoubleDouble::new(LOG2L, LOG2H), -k, r);
    z = DoubleDouble::from_exact_add(z.hi, z.lo);

    let bk = unsafe { k.to_int_unchecked::<i64>() }; /* Note: k is an integer, this is just a conversion. */
    let mk = (bk >> 12) + 0x3ff;
    let i2 = (bk >> 6) & 0x3f;
    let i1 = bk & 0x3f;

    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i2 as usize]);
    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
    let mut de = DoubleDouble::quick_mult(t1, t0);
    // exp(hi + lo) = exp(hi) * exp(lo)
    let q_hi = exp_poly_1(z.hi);
    // Taylor series exp(x) ~ 1 + x since z.lo < ulp(z.h)
    let q_lo = DoubleDouble::from_exact_add(1., z.lo);
    let q = DoubleDouble::quick_mult(q_hi, q_lo);
    de = DoubleDouble::quick_mult(de, q);
    /* we should have 1 < M < 2047 here, since we filtered out
    potential underflow/overflow cases at the beginning of this function */

    let du = (mk as u64).wrapping_shl(52);
    de.hi *= f64::from_bits(du);
    de.lo *= f64::from_bits(du);
    de
}

#[inline]
pub(crate) fn pow_exp_dd(r: DoubleDouble, s: f64) -> DoubleDouble {
    const RHO0: f64 = f64::from_bits(0xc0874910ee4e8a27);
    // #define RHO1 -0x1.577453f1799a6p+9
    /* We increase the initial value of RHO1 to avoid spurious underflow in
    the result value el. However, it is not possible to obtain a lower
    bound on |el| from the input value rh, thus this modified value of RHO1
    is obtained experimentally. */
    const RHO1: f64 = f64::from_bits(0xc08483b8cca421af);
    const RHO2: f64 = f64::from_bits(0x40862e42e709a95b);
    const RHO3: f64 = f64::from_bits(0x40862e4316ea5df9);

    // use !(rh <= RHO2) instead of rh < RHO2 to catch rh = NaN too
    #[allow(clippy::neg_cmp_op_on_partial_ord)]
    if !(r.hi <= RHO2) {
        return if r.hi > RHO3 {
            /* If rh > RHO3, we are sure there is overflow,
            For s=1 we return eh = el = DBL_MAX, which yields
            res_min = res_max = +Inf for rounding up or to nearest,
            and res_min = res_max = DBL_MAX for rounding down or toward zero,
            which will yield the correct rounding.
            For s=-1 we return eh = el = -DBL_MAX, which similarly gives
            res_min = res_max = -Inf or res_min = res_max = -DBL_MAX,
            which is the correct rounding. */
            DoubleDouble::new(
                f64::from_bits(0x7fefffffffffffff) * s,
                f64::from_bits(0x7fefffffffffffff) * s,
            )
        } else {
            /* If RHO2 < rh <= RHO3, we are in the intermediate region
            where there might be overflow or not, thus we set eh = el = NaN,
            which will set res_min = res_max = NaN, the comparison
            res_min == res_max will fail: we defer to the 2nd phase. */
            DoubleDouble::new(f64::NAN, f64::NAN)
        };
    }

    if r.hi < RHO1 {
        return if r.hi < RHO0 {
            /* For s=1, we have eh=el=+0 except for rounding up,
               thus res_min=+0 or -0, res_max=+0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=1 and rounding up, we have eh=+0, el=2^-1074,
               thus res_min = res_max = 2^-1074, which is the expected result too.
               For s=-1, we have eh=el=-0 except for rounding down,
               thus res_min=-0 or +0, res_max=-0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=-1 and rounding down, we have eh=-0, el=-2^-1074,
               thus res_min = res_max = -2^-1074, which is the expected result too.
            */
            DoubleDouble::new(f64::from_bits(0x0000000000000001) * (0.5 * s), 0.0 * s)
        } else {
            /* RHO0 <= rh < RHO1 or s < 0: we defer to the 2nd phase */
            DoubleDouble::new(f64::NAN, f64::NAN)
        };
    }
    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);

    let k = (r.hi * INVLOG2).cpu_round();

    const LOG2H: f64 = f64::from_bits(0x3f262e42fefa39ef);
    const LOG2L: f64 = f64::from_bits(0x3bbabc9e3b39803f);

    let z = DoubleDouble::mul_f64_add(DoubleDouble::new(LOG2L, LOG2H), -k, r);

    let bk = unsafe { k.to_int_unchecked::<i64>() }; /* Note: k is an integer, this is just a conversion. */
    let mk = (bk >> 12) + 0x3ff;
    let i2 = (bk >> 6) & 0x3f;
    let i1 = bk & 0x3f;

    let t0 = DoubleDouble::from_bit_pair(EXP_REDUCE_T0[i2 as usize]);
    let t1 = DoubleDouble::from_bit_pair(EXP_REDUCE_T1[i1 as usize]);
    let mut de = DoubleDouble::quick_mult(t1, t0);
    let q = exp_poly_1(z.to_f64());
    de = DoubleDouble::quick_mult(de, q);
    /* we should have 1 < M < 2047 here, since we filtered out
    potential underflow/overflow cases at the beginning of this function */

    let mut du = (mk as u64).wrapping_shl(52);
    du = (f64::from_bits(du) * s).to_bits();
    de.hi *= f64::from_bits(du);
    de.lo *= f64::from_bits(du);
    de
}
/*
#[inline(always)]
pub(crate) fn expm1_poly_dd(z: DoubleDouble) -> DoubleDouble {
    /*
       Sollya:
       pretty = proc(u) {
         return ~(floor(u*1000)/1000);
       };

       d = [-2^-12.905,2^-12.905];
       f = expm1(x);
       w = 1;
       pf = fpminimax(f, [|1,2,3,4,5,6,7|], [|1, 1, 107...|], d, absolute, floating);
       err_p = -log2(dirtyinfnorm(pf*w-f, d));
       display = decimal;

       for i from 1 to degree(pf) do print(coeff(pf, i));

       print (pf);
       display = decimal;
       print ("absolute error:",pretty(err_p));
       f = 1;
       w = 1/expm1(x);
       err_p = -log2(dirtyinfnorm(pf*w-f, d));
       print ("relative error:",pretty(err_p));
    */
    const Q: [(u64, u64); 7] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x0000000000000000, 0x3fe0000000000000),
        (0xbc75555554d7c48c, 0x3fc5555555555556),
        (0xbc555a40ffb472d9, 0x3fa5555555555556),
        (0x3c24866314c38093, 0x3f8111111111110e),
        (0x3be34665978dddb8, 0x3f56c16c16efac90),
        (0x3baeab43b813ef24, 0x3f2a01a1e12d253c),
    ];
    let z2 = z * z;
    let z4 = z2 * z2;

    let b0 = DoubleDouble::quick_mul_add(
        z,
        DoubleDouble::from_bit_pair(Q[1]),
        DoubleDouble::from_bit_pair(Q[0]),
    );
    let b1 = DoubleDouble::quick_mul_add(
        z,
        DoubleDouble::from_bit_pair(Q[3]),
        DoubleDouble::from_bit_pair(Q[2]),
    );
    let b2 = DoubleDouble::quick_mul_add(
        z,
        DoubleDouble::from_bit_pair(Q[5]),
        DoubleDouble::from_bit_pair(Q[4]),
    );

    let c0 = DoubleDouble::quick_mul_add(z2, b1, b0);
    let c1 = DoubleDouble::quick_mul_add(z2, DoubleDouble::from_bit_pair(Q[6]), b2);

    let p = DoubleDouble::quick_mul_add(z4, c1, c0);
    DoubleDouble::quick_mult(p, z)
}
*/
#[inline(always)]
pub(crate) fn expm1_poly_fast(z: DoubleDouble) -> DoubleDouble {
    // Polynomial generated by Sollya:
    // d = [-2^-12.905,2^-12.905];
    // f = expm1(x);
    // w = 1;
    // pf = fpminimax(f, [|1,2,3,4,5,6,7|], [|1, 1, D...|], d, absolute, floating);
    // See ./notes/compound_m1_expm1_fast.sollya
    let p = f_polyeval6(
        z.hi,
        f64::from_bits(0x3fe0000000000000),
        f64::from_bits(0x3fc5555555555555),
        f64::from_bits(0x3fa55555555553de),
        f64::from_bits(0x3f81111144995a9a),
        f64::from_bits(0x3f56c241f9a791c5),
        f64::from_bits(0xbfad9209c6d8b9e1),
    );
    let px = DoubleDouble::quick_mult_f64(DoubleDouble::from_exact_mult(z.hi, p), z.hi);
    // expm1(hi + lo) = expm1(hi) + expm1(lo)(1 + expm1(hi)) = expm1(hi) + expm1(lo)expm1(hi) + expm1(lo)
    // expm1(lo) ~ lo
    let expm1_hi = DoubleDouble::f64_add(z.hi, px);
    let mut lowest_part = DoubleDouble::quick_mult_f64(expm1_hi, z.lo);
    lowest_part = DoubleDouble::full_add_f64(lowest_part, z.lo);
    DoubleDouble::quick_dd_add(expm1_hi, lowest_part)
}

/// |z.hi| < 2^-7
#[inline(always)]
pub(crate) fn expm1_poly_dd_tiny(z: DoubleDouble) -> DoubleDouble {
    // Polynomial generated in Sollya
    // d = [-2^-7,2^-7];
    // f = expm1(x);
    // w = 1;
    // pf = fpminimax(f, [|1,2,3,4,5,6,7,8,9|], [|1, 1, 107...|], d, absolute, floating);
    // See ./notes/compound_expm1_tiny.sollya
    const Q: [(u64, u64); 9] = [
        (0x0000000000000000, 0x3ff0000000000000),
        (0x0000000000000000, 0x3fe0000000000000),
        (0x3c6555564150ff16, 0x3fc5555555555555),
        (0x3c4586275c26f8a5, 0x3fa5555555555555),
        (0xbc19e6193ac658a6, 0x3f81111111111111),
        (0xbbf025e72dc21051, 0x3f56c16c16c1500a),
        (0x3bc2d641a7b7b9b8, 0x3f2a01a01a07dc46),
        (0xbb42cc8aaeeb3d00, 0x3efa01a29fef3e6f),
        (0x3b52b1589125ce82, 0x3ec71db6af553255),
    ];
    let z = DoubleDouble::from_exact_add(z.hi, z.lo);
    let mut d = DoubleDouble::quick_mul_add(
        z,
        DoubleDouble::from_bit_pair(Q[8]),
        DoubleDouble::from_bit_pair(Q[7]),
    );
    d = DoubleDouble::quick_mul_add(z, d, DoubleDouble::from_bit_pair(Q[6]));
    d = DoubleDouble::quick_mul_add(z, d, DoubleDouble::from_bit_pair(Q[5]));
    d = DoubleDouble::quick_mul_add(z, d, DoubleDouble::from_bit_pair(Q[4]));
    d = DoubleDouble::quick_mul_add(z, d, DoubleDouble::from_bit_pair(Q[3]));
    d = DoubleDouble::quick_mul_add(z, d, DoubleDouble::from_bit_pair(Q[2]));
    d = DoubleDouble::quick_mul_add_f64(z, d, f64::from_bits(0x3fe0000000000000));
    d = DoubleDouble::quick_mul_add_f64(z, d, f64::from_bits(0x3ff0000000000000));
    DoubleDouble::quick_mult(d, z)
}

#[inline]
pub(crate) fn pow_expm1_1(r: DoubleDouble, s: f64) -> DoubleDouble {
    const RHO0: f64 = f64::from_bits(0xc0874910ee4e8a27);
    // #define RHO1 -0x1.577453f1799a6p+9
    /* We increase the initial value of RHO1 to avoid spurious underflow in
    the result value el. However, it is not possible to obtain a lower
    bound on |el| from the input value rh, thus this modified value of RHO1
    is obtained experimentally. */
    const RHO1: f64 = f64::from_bits(0xc08483b8cca421af);
    const RHO2: f64 = f64::from_bits(0x40862e42e709a95b);
    const RHO3: f64 = f64::from_bits(0x40862e4316ea5df9);

    // use !(rh <= RHO2) instead of rh < RHO2 to catch rh = NaN too
    #[allow(clippy::neg_cmp_op_on_partial_ord)]
    if !(r.hi <= RHO2) {
        return if r.hi > RHO3 {
            /* If rh > RHO3, we are sure there is overflow,
            For s=1 we return eh = el = DBL_MAX, which yields
            res_min = res_max = +Inf for rounding up or to nearest,
            and res_min = res_max = DBL_MAX for rounding down or toward zero,
            which will yield the correct rounding.
            For s=-1 we return eh = el = -DBL_MAX, which similarly gives
            res_min = res_max = -Inf or res_min = res_max = -DBL_MAX,
            which is the correct rounding. */
            DoubleDouble::new(
                f64::from_bits(0x7fefffffffffffff) * s,
                f64::from_bits(0x7fefffffffffffff) * s,
            )
        } else {
            /* If RHO2 < rh <= RHO3, we are in the intermediate region
            where there might be overflow or not, thus we set eh = el = NaN,
            which will set res_min = res_max = NaN, the comparison
            res_min == res_max will fail: we defer to the 2nd phase. */
            DoubleDouble::new(f64::NAN, f64::NAN)
        };
    }

    if r.hi < RHO1 {
        if r.hi < RHO0 {
            /* For s=1, we have eh=el=+0 except for rounding up,
               thus res_min=+0 or -0, res_max=+0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=1 and rounding up, we have eh=+0, el=2^-1074,
               thus res_min = res_max = 2^-1074, which is the expected result too.
               For s=-1, we have eh=el=-0 except for rounding down,
               thus res_min=-0 or +0, res_max=-0 in the main code,
               the rounding test succeeds, and we return res_max which is the
               expected result in the underflow case.
               For s=-1 and rounding down, we have eh=-0, el=-2^-1074,
               thus res_min = res_max = -2^-1074, which is the expected result too.
            */
            return DoubleDouble::full_add_f64(
                DoubleDouble::new(f64::from_bits(0x0000000000000001) * (0.5 * s), 0.0 * s),
                -1.0,
            );
        } else {
            /* RHO0 <= rh < RHO1 or s < 0: we return -1 */
            return DoubleDouble::new(0., -1.);
        };
    }

    let ax = r.hi.to_bits() & 0x7fffffffffffffffu64;

    const LOG2H: f64 = f64::from_bits(0x3f262e42fefa39ef);
    const LOG2L: f64 = f64::from_bits(0x3bbabc9e3b39803f);

    if ax <= 0x3f80000000000000 {
        // |x| < 2^-7
        if ax < 0x3970000000000000 {
            // |x| < 2^-104
            return r;
        }
        let d = expm1_poly_dd_tiny(r);
        return d;
    }

    const INVLOG2: f64 = f64::from_bits(0x40b71547652b82fe);

    let k = (r.hi * INVLOG2).cpu_round_ties_even();

    let z = DoubleDouble::mul_f64_add(DoubleDouble::new(LOG2L, LOG2H), -k, r);

    let bk = unsafe { k.to_int_unchecked::<i64>() }; /* Note: k is an integer, this is just a conversion. */
    let mk = (bk >> 12) + 0x3ff;
    let i2 = (bk >> 6) & 0x3f;
    let i1 = bk & 0x3f;

    let t0 = DoubleDouble::from_bit_pair(EXPM1_T0[i2 as usize]);
    let t1 = DoubleDouble::from_bit_pair(EXPM1_T1[i1 as usize]);
    let tbh = DoubleDouble::quick_mult(t1, t0);
    let mut de = tbh;
    // exp(k)=2^k*exp(r) + (2^k - 1)
    let q = expm1_poly_fast(z);
    de = DoubleDouble::quick_mult(de, q);
    de = DoubleDouble::add(tbh, de);

    let ie = mk - 0x3ff;

    let off: f64 = f64::from_bits((2048i64 + 1023i64).wrapping_sub(ie).wrapping_shl(52) as u64);

    let e: f64;
    if ie < 53 {
        let fhz = DoubleDouble::from_exact_add(off, de.hi);
        de.hi = fhz.hi;
        e = fhz.lo;
    } else if ie < 104 {
        let fhz = DoubleDouble::from_exact_add(de.hi, off);
        de.hi = fhz.hi;
        e = fhz.lo;
    } else {
        e = 0.;
    }
    de.lo += e;
    de.hi = ldexp(de.to_f64(), ie as i32);
    de.lo = 0.;
    de
}

fn exp_dyadic_poly(x: DyadicFloat128) -> DyadicFloat128 {
    const Q_2: [DyadicFloat128; 8] = [
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -128,
            mantissa: 0xffff_ffff_ffff_ffff_ffff_ffff_ffff_ffd0_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -127,
            mantissa: 0x8000_0000_0000_0000_0000_0000_0000_0088_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -128,
            mantissa: 0x8000_0000_0000_0000_0000_000c_06f3_cd29_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -130,
            mantissa: 0xaaaa_aaaa_aaaa_aaaa_aaaa_aa6a_1e07_76ae_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -132,
            mantissa: 0xaaaa_aaaa_aaaa_aaa3_0000_0000_0000_0000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -134,
            mantissa: 0x8888_8888_8888_8897_0000_0000_0000_0000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -137,
            mantissa: 0xb60b_60b9_3214_6a54_0000_0000_0000_0000_u128,
        },
        DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: -140,
            mantissa: 0xd00d_00cd_9841_6862_0000_0000_0000_0000_u128,
        },
    ];
    let mut p = Q_2[7];
    for i in (0..7).rev() {
        p = x * p + Q_2[i];
    }
    p
}

/* put in r an approximation of exp(x), for |x| < 744.45,
with relative error < 2^-121.70 */
#[inline]
pub(crate) fn exp_dyadic(x: DyadicFloat128) -> DyadicFloat128 {
    let ex = x.exponent + 127;
    if ex >= 10
    // underflow or overflow
    {
        return DyadicFloat128 {
            sign: DyadicSign::Pos,
            exponent: if x.sign == DyadicSign::Neg {
                -1076
            } else {
                1025
            },
            mantissa: x.mantissa,
        };
    }

    const LOG2_INV: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -115,
        mantissa: 0xb8aa_3b29_5c17_f0bc_0000_0000_0000_0000_u128,
    };

    const LOG2: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -128,
        mantissa: 0xb172_17f7_d1cf_79ab_c9e3_b398_03f2_f6af_u128,
    };

    let mut bk = x * LOG2_INV;
    let k = bk.trunc_to_i64(); /* k = trunc(K) [rounded towards zero, exact] */
    /* The rounding error of mul_dint_int64() is bounded by 6 ulps, thus since
    |K| <= 4399162*log(2) < 3049267, the error on K is bounded by 2^-103.41.
    This error is divided by 2^12 below, thus yields < 2^-115.41. */
    bk = LOG2.mul_int64(k);
    bk.exponent -= 12;
    bk.sign = bk.sign.negate();
    let y = x + bk;

    let bm = k >> 12;
    let i2 = (k >> 6) & 0x3f;
    let i1 = k & 0x3f;
    let mut r = exp_dyadic_poly(y);
    r = EXP_T1_2_DYADIC[i2 as usize] * r;
    r = EXP_T2_2_DYADIC[i1 as usize] * r;
    r.exponent += bm as i16; /* exact */
    r
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::f_expm1;
    #[test]
    fn test_log() {
        let k = DyadicFloat128::new_from_f64(2.5);
        assert_eq!(exp_dyadic(k).fast_as_f64(), 12.182493960703473);
    }

    #[test]
    fn test_exp() {
        let k = pow_expm1_1(DoubleDouble::new(0., 2.543543543543), 1.);
        println!("{}", k.to_f64());
        println!("{}", f_expm1(2.543543543543));
    }
}