/*

 * // 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::bits::{biased_exponent_f64, get_exponent_f64, mantissa_f64};

use crate::common::{dd_fmla, dyad_fmla, f_fmla};

use crate::double_double::DoubleDouble;

use crate::dyadic_float::{DyadicFloat128, DyadicSign};

use crate::logs::log2p1_dyadic_tables::{LOG2P1_F128_POLY, LOG2P1_INVERSE_2, LOG2P1_LOG_INV_2};

use crate::logs::log2p1_tables::{LOG2P1_EXACT, LOG2P1_INVERSE, LOG2P1_LOG_DD_INVERSE};

/* put in h+l a double-double approximation of log(z)-z for

|z| < 0.03125, with absolute error bounded by 2^-67.14

(see analyze_p1a(-0.03125,0.03125) from log1p.sage) */

#[inline]

pub(crate) fn log_p_1a(z: f64) -> DoubleDouble {

    let z2: DoubleDouble = if z.abs() >= f64::from_bits(0x3000000000000000) {

        DoubleDouble::from_exact_mult(z, z)

    } else {

        // avoid spurious underflow

        DoubleDouble::default()

    };

    let z4h = z2.hi * z2.hi;

    /* The following is a degree-11 polynomial generated by Sollya

    approximating log(1+x) for |x| < 0.03125,

    with absolute error < 2^-73.441 and relative error < 2^-67.088

    (see file Pabs_a.sollya).

    The polynomial is P[0]*x + P[1]*x^2 + ... + P[10]*x^11.

    The algorithm assumes that the degree-1 coefficient P[0] is 1

    and the degree-2 coefficient P[1] is -0.5. */

    const PA: [u64; 11] = [

        0x3ff0000000000000,

        0xbfe0000000000000,

        0x3fd5555555555555,

        0xbfcffffffffffe5f,

        0x3fc999999999aa82,

        0xbfc555555583a8c8,

        0x3fc2492491c359e6,

        0xbfbffffc728edeea,

        0x3fbc71c961f34980,

        0xbfb9a82ac77c05f4,

        0x3fb74b40dd1707d3,

    ];

    let p910 = dd_fmla(f64::from_bits(PA[10]), z, f64::from_bits(PA[9]));

    let p78 = dd_fmla(f64::from_bits(PA[8]), z, f64::from_bits(PA[7]));

    let p56 = dd_fmla(f64::from_bits(PA[6]), z, f64::from_bits(PA[5]));

    let p34 = dd_fmla(f64::from_bits(PA[4]), z, f64::from_bits(PA[3]));

    let p710 = dd_fmla(p910, z2.hi, p78);

    let p36 = dd_fmla(p56, z2.hi, p34);

    let mut ph = dd_fmla(p710, z4h, p36);

    ph = dd_fmla(ph, z, f64::from_bits(PA[2]));

    ph *= z2.hi;

    let mut p = DoubleDouble::from_exact_add(-0.5 * z2.hi, ph * z);

    p.lo += -0.5 * z2.lo;

    p

}

/* put in h+l a double-double approximation of log(z)-z for

|z| < 0.00212097167968735, with absolute error bounded by 2^-78.25

(see analyze_p1(-0.00212097167968735,0.00212097167968735)

from accompanying file log1p.sage, which also yields |l| < 2^-69.99) */

#[inline]

fn p_1(z: f64) -> DoubleDouble {

    const P: [u64; 7] = [

        0x3ff0000000000000,

        0xbfe0000000000000,

        0x3fd5555555555550,

        0xbfcfffffffff572d,

        0x3fc999999a2d7868,

        0xbfc5555c0d31b08e,

        0x3fc2476b9058e396,

    ];

    let z2 = DoubleDouble::from_exact_mult(z, z);

    let p56 = dd_fmla(f64::from_bits(P[6]), z, f64::from_bits(P[5]));

    let p34 = dd_fmla(f64::from_bits(P[4]), z, f64::from_bits(P[3]));

    let mut ph = dd_fmla(p56, z2.hi, p34);

    /* ph approximates P[3]+P[4]*z+P[5]*z^2+P[6]*z^3 */

    ph = dd_fmla(ph, z, f64::from_bits(P[2]));

    /* ph approximates P[2]+P[3]*z+P[4]*z^2+P[5]*z^3+P[6]*z^4 */

    ph *= z2.hi;

    /* ph approximates P[2]*z^2+P[3]*z^3+P[4]*z^4+P[5]*z^5+P[6]*z^6 */

    let mut p = DoubleDouble::from_exact_add(-0.5 * z2.hi, ph * z);

    p.lo += -0.5 * z2.lo;

    p

}

#[inline]

pub(crate) fn log_fast(e: i32, v_u: u64) -> DoubleDouble {

    let m: u64 = 0x10000000000000u64.wrapping_add(v_u & 0xfffffffffffff);

    /* x = m/2^52 */

    /* if x > sqrt(2), we divide it by 2 to avoid cancellation */

    let c: i32 = if m >= 0x16a09e667f3bcd { 1 } else { 0 };

    let e = e.wrapping_add(c); /* now -1074 <= e <= 1024 */

    static CY: [f64; 2] = [1.0, 0.5];

    static CM: [u32; 2] = [43, 44];

    let i: i32 = (m >> CM[c as usize]) as i32;

    let y = f64::from_bits(v_u) * CY[c as usize];

    const OFFSET: i32 = 362;

    let r = f64::from_bits(LOG2P1_INVERSE[(i - OFFSET) as usize]);

    let log2_inv_dd = LOG2P1_LOG_DD_INVERSE[(i - OFFSET) as usize];

    let l1 = f64::from_bits(log2_inv_dd.1);

    let l2 = f64::from_bits(log2_inv_dd.0);

    let z = dd_fmla(r, y, -1.0); /* exact */

    /* evaluate P(z), for |z| < 0.00212097167968735 */

    let p = p_1(z);

    /* Add e*log(2) to (h,l), where -1074 <= e <= 1023, thus e has at most

    11 bits. log2_h is an integer multiple of 2^-42, so that e*log2_h

    is exact. */

    const LOG2_H: f64 = f64::from_bits(0x3fe62e42fefa3800);

    const LOG2_L: f64 = f64::from_bits(0x3d2ef35793c76730);

    let ee = e as f64;

    let mut vl = DoubleDouble::from_exact_add(f_fmla(ee, LOG2_H, l1), z);

    vl.lo = p.hi + (vl.lo + (l2 + p.lo));

    vl.lo = dd_fmla(ee, LOG2_L, vl.lo);

    vl

}

const INV_LOG2_DD: DoubleDouble = DoubleDouble::new(

    f64::from_bits(0x3c7777d0ffda0d24),

    f64::from_bits(0x3ff71547652b82fe),

);

fn log2p1_accurate_small(x: f64) -> f64 {

    static P_ACC: [u64; 24] = [

        0x3ff71547652b82fe,

        0x3c7777d0ffda0d24,

        0xbfe71547652b82fe,

        0xbc6777d0ffd9ddb8,

        0x3fdec709dc3a03fd,

        0x3c7d27f055481523,

        0xbfd71547652b82fe,

        0xbc5777d1456a14c4,

        0x3fd2776c50ef9bfe,

        0x3c7e4b2a04f81513,

        0xbfcec709dc3a03fd,

        0xbc6d2072e751087a,

        0x3fca61762a7aded9,

        0x3c5f90f4895378ac,

        0xbfc71547652b8301,

        0x3fc484b13d7c02ae,

        0xbfc2776c50ef7591,

        0x3fc0c9a84993cabb,

        0xbfbec709de7b1612,

        0x3fbc68f56ba73fd1,

        0xbfba616c83da87e7,

        0x3fb89f3042097218,

        0xbfb72b376930a3fa,

        0x3fb5d0211d5ab530,

    ];

    /* for degree 11 or more, ulp(c[d]*x^d) < 2^-105.5*|log2p1(x)|

    where c[d] is the degree-d coefficient of Pacc, thus we can compute

    with a double only */

    let mut h = dd_fmla(f64::from_bits(P_ACC[23]), x, f64::from_bits(P_ACC[22])); // degree 16

    for i in (11..=15).rev() {

        h = dd_fmla(h, x, f64::from_bits(P_ACC[(i + 6) as usize])); // degree i

    }

    let mut l = 0.;

    for i in (8..=10).rev() {

        let mut p = DoubleDouble::quick_f64_mult(x, DoubleDouble::new(l, h));

        l = p.lo;

        p = DoubleDouble::from_exact_add(f64::from_bits(P_ACC[(i + 6) as usize]), p.hi);

        h = p.hi;

        l += p.lo;

    }

    for i in (1..=7).rev() {

        let mut p = DoubleDouble::quick_f64_mult(x, DoubleDouble::new(l, h));

        l = p.lo;

        p = DoubleDouble::from_exact_add(f64::from_bits(P_ACC[(2 * i - 2) as usize]), p.hi);

        h = p.hi;

        l += p.lo + f64::from_bits(P_ACC[(2 * i - 1) as usize]);

    }

    let pz = DoubleDouble::quick_f64_mult(x, DoubleDouble::new(l, h));

    pz.to_f64()

}

/* deal with |x| < 2^-900, then log2p1(x) ~ x/log(2) */

#[cold]

fn log2p1_accurate_tiny(x: f64) -> f64 {

    // exceptional values

    if x.abs() == f64::from_bits(0x0002c316a14459d8) {

        return if x > 0. {

            dd_fmla(

                f64::from_bits(0x1a70000000000000),

                f64::from_bits(0x1a70000000000000),

                f64::from_bits(0x0003fc1ce8b1583f),

            )

        } else {

            dd_fmla(

                f64::from_bits(0x9a70000000000000),

                f64::from_bits(0x1a70000000000000),

                f64::from_bits(0x8003fc1ce8b1583f),

            )

        };

    }

    /* first scale x to avoid truncation of l in the underflow region */

    let sx = x * f64::from_bits(0x4690000000000000);

    let mut zh = DoubleDouble::quick_f64_mult(sx, INV_LOG2_DD);

    let res = zh.to_f64() * f64::from_bits(0x3950000000000000); // expected result

    zh.lo += dd_fmla(-res, f64::from_bits(0x4690000000000000), zh.hi);

    // the correction to apply to res is l*2^-106

    /* For all rounding modes, we have underflow

    for |x| <= 0x1.62e42fefa39eep-1023 */

    dyad_fmla(zh.lo, f64::from_bits(0x3950000000000000), res)

}

/* Given x > -1, put in (h,l) a double-double approximation of log2(1+x),

   and return a bound err on the maximal absolute error so that:

   |h + l - log2(1+x)| < err.

   We have x = m*2^e with 1 <= m < 2 (m = v.f) and -1074 <= e <= 1023.

   This routine is adapted from cr_log1p_fast.

*/

#[inline]

fn log2p1_fast(x: f64, e: i32) -> (DoubleDouble, f64) {

    if e < -5

    /* e <= -6 thus |x| < 2^-5 */

    {

        if e <= -969 {

            /* then |x| might be as small as 2^-969, thus h=x/log(2) might in the

            binade [2^-969,2^-968), with ulp(h) = 2^-1021, and if |l| < ulp(h),

            then l.ulp() might be smaller than 2^-1074. We defer that case to

            the accurate path. */

            // *h = *l = 0;

            // return 1;

            let ax = x.abs();

            let result = if ax < f64::from_bits(0x3960000000000000) {

                log2p1_accurate_tiny(x)

            } else {

                log2p1_accurate_small(x)

            };

            return (DoubleDouble::new(0.0, result), 0.0);

        }

        let mut p = log_p_1a(x);

        let p_lo = p.lo;

        p = DoubleDouble::from_exact_add(x, p.hi);

        p.lo += p_lo;

        p = DoubleDouble::quick_mult(p, INV_LOG2_DD);

        return (p, f64::from_bits(0x3c1d400000000000) * p.hi); /* 2^-61.13 < 0x1.d4p-62 */

    }

    /* (xh,xl) <- 1+x */

    let zx = DoubleDouble::from_full_exact_add(1.0, x);

    let mut v_u = zx.hi.to_bits();

    let e = ((v_u >> 52) as i32).wrapping_sub(0x3ff);

    v_u = (0x3ffu64 << 52) | (v_u & 0xfffffffffffff);

    let mut p = log_fast(e, v_u);

    /* log(xh+xl) = log(xh) + log(1+xl/xh) */

    let c = if zx.hi <= f64::from_bits(0x7fd0000000000000) || zx.lo.abs() >= 4.0 {

        zx.lo / zx.hi

    } else {

        0.

    }; // avoid spurious underflow

    /* Since |xl| < ulp(xh), we have |xl| < 2^-52 |xh|,

    thus |c| < 2^-52, and since |log(1+x)-x| < x^2 for |x| < 0.5,

    we have |log(1+c)-c)| < c^2 < 2^-104. */

    p.lo += c;

    /* now multiply h+l by 1/log(2) */

    p = DoubleDouble::quick_mult(p, INV_LOG2_DD);

    (p, f64::from_bits(0x3bb2300000000000)) /* 2^-67.82 < 0x1.23p-68 */

}

fn log_dyadic_taylor_poly(x: DyadicFloat128) -> DyadicFloat128 {

    let mut r = LOG2P1_F128_POLY[12];

    for i in (0..12).rev() {

        r = x * r + LOG2P1_F128_POLY[i];

    }

    r * x

}

pub(crate) fn log2_dyadic(d: DyadicFloat128, x: f64) -> DyadicFloat128 {

    let biased_exp = biased_exponent_f64(x);

    let e = get_exponent_f64(x);

    let base_mant = mantissa_f64(x);

    let mant = base_mant + if biased_exp != 0 { 1u64 << 52 } else { 0 };

    let lead = mant.leading_zeros();

    let kk = e - (if lead > 11 { lead - 12 } else { 0 }) as i64;

    let mut fe: i16 = kk as i16;

    let adjusted_mant = mant << lead;

    // Find the lookup index

    let mut i: i16 = (adjusted_mant >> 55) as i16;

    if adjusted_mant > 0xb504f333f9de6484 {

        fe = fe.wrapping_add(1);

        i >>= 1;

    }

    let mut x = d;

    x.exponent = x.exponent.wrapping_sub(fe);

    let inverse_2 = LOG2P1_INVERSE_2[(i - 128) as usize];

    let mut z = x * inverse_2;

    const F128_MINUS_ONE: DyadicFloat128 = DyadicFloat128 {

        sign: DyadicSign::Neg,

        exponent: -127,

        mantissa: 0x8000_0000_0000_0000_0000_0000_0000_0000_u128,

    };

    z = z + F128_MINUS_ONE;

    const LOG2: DyadicFloat128 = DyadicFloat128 {

        sign: DyadicSign::Pos,

        exponent: -128,

        mantissa: 0xb172_17f7_d1cf_79ab_c9e3_b398_03f2_f6af_u128,

    };

    // E·log(2)

    let r = LOG2.mul_int64(fe as i64);

    let mut p = log_dyadic_taylor_poly(z);

    p = LOG2P1_LOG_INV_2[(i - 128) as usize] + p;

    p + r

}

#[cold]

fn log2p1_accurate(x: f64) -> f64 {

    let ax = x.abs();

    if ax < f64::from_bits(0x3fa0000000000000) {

        return if ax < f64::from_bits(0x3960000000000000) {

            log2p1_accurate_tiny(x)

        } else {

            log2p1_accurate_small(x)

        };

    }

    let dx = if x > 1.0 {

        DoubleDouble::from_exact_add(x, 1.0)

    } else {

        DoubleDouble::from_exact_add(1.0, x)

    };

    /* log2p1(x) is exact when 1+x = 2^e, thus when 2^e-1 is exactly

    representable. This can only occur when xl=0 here. */

    let mut t: u64 = x.to_bits();

    if dx.lo == 0. {

        /* check if xh is a power of two */

        t = dx.hi.to_bits();

        if (t.wrapping_shl(12)) == 0 {

            let e = ((t >> 52) as i32).wrapping_sub(0x3ff);

            return e as f64;

        }

    }

    /* if x=2^e, the accurate path will fail for directed roundings */

    if (t.wrapping_shl(12)) == 0 {

        let e: i32 = ((t >> 52) as i32).wrapping_sub(0x3ff); // x = 2^e

        /* for e >= 49, log2p1(x) rounds to e for rounding to nearest;

        for e >= 48, log2p1(x) rounds to e for rounding toward zero;

        for e >= 48, log2p1(x) rounds to nextabove(e) for rounding up;

        for e >= 48, log2p1(x) rounds to e for rounding down. */

        if e >= 49 {

            return e as f64 + f64::from_bits(0x3cf0000000000000); // 0x1p-48 = 1/2 ulp(49)

        }

    }

    let x_d = DyadicFloat128::new_from_f64(dx.hi);

    let mut y = log2_dyadic(x_d, dx.hi);

    let mut c = DyadicFloat128::from_div_f64(dx.lo, dx.hi);

    let mut bx = c * c;

    /* multiply X by -1/2 */

    bx.exponent -= 1;

    bx.sign = DyadicSign::Neg;

    /* C <- C - C^2/2 */

    c = c + bx;

    /* |C-log(1+xl/xh)| ~ 2e-64 */

    y = y + c;

    const LOG2_INV: DyadicFloat128 = DyadicFloat128 {

        sign: DyadicSign::Pos,

        exponent: -115,

        mantissa: 0xb8aa_3b29_5c17_f0bb_be87_fed0_691d_3e89_u128,

    };

    y = y * LOG2_INV;

    y.exponent -= 12;

    y.fast_as_f64()

}

/// Computes log2(x+1)

///

/// Max ULP 0.5

pub fn f_log2p1(x: f64) -> f64 {

    let x_u = x.to_bits();

    let e = (((x_u >> 52) & 0x7ff) as i32).wrapping_sub(0x3ff);

    if e == 0x400 || x == 0. || x <= -1.0 {

        /* case NaN/Inf, +/-0 or x <= -1 */

        if e == 0x400 && x.to_bits() != 0xfffu64 << 52 {

            /* NaN or + Inf*/

            return x + x;

        }

        if x <= -1.0

        /* we use the fact that NaN < -1 is false */

        {

            /* log2p(x<-1) is NaN, log2p(-1) is -Inf and raises DivByZero */

            return if x < -1.0 {

                f64::NAN

            } else {

                // x=-1

                f64::NEG_INFINITY

            };

        }

        return x + x; /* +/-0 */

    }

    /* now x > -1 */

    /* check x=2^n-1 for 0 <= n <= 53, where log2p1(x) is exact,

    and we shouldn't raise the inexact flag */

    if 0 <= e && e <= 52 {

        /* T[e]=2^(e+1)-1, i.e., the unique value of the form 2^n-1

        in the interval [2^e, 2^(e+1)). */

        if x == f64::from_bits(LOG2P1_EXACT[e as usize]) {

            return (e + 1) as f64;

        }

    }

    /* For x=2^k-1, -53 <= k <= -1, log2p1(x) = k is also exact. */

    if e == -1 && x < 0. {

        // -1 < x <= -1/2

        let w = (1.0 + x).to_bits(); // 1+x is exact

        if w.wrapping_shl(12) == 0 {

            // 1+x = 2^k

            let k: i32 = ((w >> 52) as i32).wrapping_sub(0x3ff);

            return k as f64;

        }

    }

    /* now x = m*2^e with 1 <= m < 2 (m = v.f) and -1074 <= e <= 1023 */

    let (p, err) = log2p1_fast(x, e);

    let left = p.hi + (p.lo - err);

    let right = p.hi + (p.lo + err);

    if left == right {

        return left;

    }

    log2p1_accurate(x)

}

#[cfg(test)]

mod tests {

    use super::*;

    #[test]

    fn test_log2p1() {

        assert_eq!(f_log2p1(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000008344095884546873),

                   0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012037985753337781);

        assert_eq!(f_log2p1(0.00006669877554532304), 0.00009622278377734607);

        assert_eq!(f_log2p1(1.00006669877554532304), 1.0000481121941047);

        assert_eq!(f_log2p1(-0.90006669877554532304), -3.322890675865049);

    }

}