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::*;
use crate::double_double::DoubleDouble;
use crate::dyadic_float::{DyadicFloat128, DyadicSign};
use crate::logs::log2p1::{log_fast, log_p_1a, log2_dyadic};
use crate::logs::log10p1_tables::{LOG10P1_EXACT_INT_S_TABLE, LOG10P1_EXACT_INT_TABLE};

const INV_LOG10_DD: DoubleDouble = DoubleDouble::new(
    f64::from_bits(0x3c695355baaafad3),
    f64::from_bits(0x3fdbcb7b1526e50e),
);

/* deal with |x| < 2^-900, then log10p1(x) ~ x/log(10) */
#[cold]
fn log10p1_accurate_tiny(x: f64) -> f64 {
    /* first scale x to avoid truncation of l in the underflow region */
    let sx = x * f64::from_bits(0x4690000000000000);
    let mut px = DoubleDouble::quick_f64_mult(sx, INV_LOG10_DD);

    let res = px.to_f64() * f64::from_bits(0x3950000000000000); // expected result
    px.lo += dd_fmla(-res, f64::from_bits(0x4690000000000000), px.hi);
    // the correction to apply to res is l*2^-106
    /* For RNDN, we have underflow for |x| <= 0x1.26bb1bbb55515p-1021,
    and for rounding away, for |x| < 0x1.26bb1bbb55515p-1021. */

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

fn log10p1_accurate_small(x: f64) -> f64 {
    /* the following is a degree-17 polynomial approximating log10p1(x) for
    |x| <= 2^-5 with relative error < 2^-105.067*/

    static P_ACC: [u64; 25] = [
        0x3fdbcb7b1526e50e,
        0x3c695355baaafad4,
        0xbfcbcb7b1526e50e,
        0xbc595355baaae078,
        0x3fc287a7636f435f,
        0xbc59c871838f83ac,
        0xbfbbcb7b1526e50e,
        0xbc495355e23285f2,
        0x3fb63c62775250d8,
        0x3c4442abd5831422,
        0xbfb287a7636f435f,
        0x3c49d116f225c4e4,
        0x3fafc3fa615105c7,
        0x3c24e1d7b4790510,
        0xbfabcb7b1526e512,
        0x3c49f884199ab0ce,
        0x3fa8b4df2f3f0486,
        0xbfa63c6277522391,
        0x3fa436e526a79e5c,
        0xbfa287a764c5a762,
        0x3fa11ac1e784daec,
        0xbf9fc3eedc920817,
        0x3f9da5cac3522edb,
        0xbf9be5ca1f9a97cd,
        0x3f9a44b64ca06e9b,
    ];

    /* for degree 11 or more, ulp(c[d]*x^d) < 2^-105.7*|log10p1(x)|
    where c[d] is the degree-d coefficient of Pacc, thus we can compute
    with a double only, and even with degree 10 (this does not increase
    the number of exceptional cases) */

    let mut h = dd_fmla(f64::from_bits(P_ACC[24]), x, f64::from_bits(P_ACC[23])); // degree 16
    for i in (10..=15).rev() {
        h = dd_fmla(h, x, f64::from_bits(P_ACC[(i + 7) as usize])); // degree i
    }

    // degree 9
    let px = DoubleDouble::from_exact_mult(x, h);
    let hl = DoubleDouble::from_exact_add(f64::from_bits(P_ACC[9 + 7]), px.hi);
    h = hl.hi;
    let mut l = px.lo + hl.lo;

    for i in (1..=8).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()
}

#[cold]
fn log10p1_accurate(x: f64) -> f64 {
    let ax = x.abs();

    if ax < f64::from_bits(0x3fa0000000000000) {
        return if ax < f64::from_bits(0x07b0000000000000) {
            log10p1_accurate_tiny(x)
        } else {
            log10p1_accurate_small(x)
        };
    }
    let dx = if x > 1.0 {
        DoubleDouble::from_exact_add(x, 1.0)
    } else {
        DoubleDouble::from_exact_add(1.0, x)
    };
    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;

    // Sage Math:
    // from sage.all import *
    //
    // def format_hex2(value):
    //     l = hex(value)[2:]
    //     n = 4
    //     x = [l[i:i + n] for i in range(0, len(l), n)]
    //     return "0x" + "_".join(x) + "_u128"
    // (s, m, e) = (RealField(128)(1)/RealField(128)(10)).log().sign_mantissa_exponent();
    // print(format_hex2(m));
    const LOG10_INV: DyadicFloat128 = DyadicFloat128 {
        sign: DyadicSign::Pos,
        exponent: -129,
        mantissa: 0xde5b_d8a9_3728_7195_355b_aaaf_ad33_dc32_u128,
    };
    y = y * LOG10_INV;
    y.fast_as_f64()
}

#[inline]
fn log10p1_fast(x: f64, e: i32) -> (DoubleDouble, f64) {
    if e < -5
    /* e <= -6 thus |x| < 2^-5 */
    {
        if e <= -968 {
            /* then |x| might be as small as 2^-968, thus h=x/log(10) might in the
            binade [2^-970,2^-969), with ulp(h) = 2^-1022, and if |l| < ulp(h),
            then l.ulp() might be smaller than 2^-1074. We defer that case to
            the accurate path. */
            let ax = x.abs();
            let result = if ax < f64::from_bits(0x07b0000000000000) {
                log10p1_accurate_tiny(x)
            } else {
                log10p1_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_LOG10_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(x, 1.0);

    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_LOG10_DD);
    (p, f64::from_bits(0x3bb0a00000000000)) /* 2^-67.92 < 0x1.0ap-68 */
}

/// Computes log10(x+1)
///
/// Max ULP 0.5
pub fn f_log10p1(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 */
    }

    /* check x=10^n-1 for 1 <= n <= 15, where log10p1(x) is exact,
    and we shouldn't raise the inexact flag */
    if 3 <= e && e <= 49 && x == f64::from_bits(LOG10P1_EXACT_INT_TABLE[e as usize]) {
        return LOG10P1_EXACT_INT_S_TABLE[e as usize] as f64;
    }

    /* now x = m*2^e with 1 <= m < 2 (m = v.f) and -1074 <= e <= 1023 */
    let (p, err) = log10p1_fast(x, e);
    let left = p.hi + (p.lo - err);
    let right = p.hi + (p.lo + err);
    if left == right {
        return left;
    }

    log10p1_accurate(x)
}

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

    #[test]
    fn test_log10p1() {
        assert_eq!(f_log10p1(0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000013904929147106097),
                   0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000006038833999843867 );
        assert!(f_log10p1(-2.0).is_nan());
        assert_eq!(f_log10p1(9.0), 1.0);
        assert_eq!(f_log10p1(2.0), 0.47712125471966244);
        assert_eq!(f_log10p1(-0.5), -0.3010299956639812);
    }
}