pxfm 0.1.28

/*
 * // 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, f_fmla};
use crate::double_double::DoubleDouble;
use crate::hyperbolic::acosh::{
    ACOSH_ASINH_B, ACOSH_ASINH_LL, ACOSH_ASINH_R1, ACOSH_ASINH_R2, ACOSH_ASINH_REFINE_T2,
    ACOSH_ASINH_REFINE_T4, ACOSH_SINH_REFINE_T1, ACOSH_SINH_REFINE_T3, lpoly_xd_generic,
};

static ATANH_L1: [(u64, u64); 33] = [
    (0x0000000000000000, 0x0000000000000000),
    (0xbe4532c1269e2038, 0x3f862e5000000000),
    (0x3e4ce42d81b54e84, 0x3f962e3c00000000),
    (0xbe525826f815ec3d, 0x3fa0a2ac00000000),
    (0x3e50db1b1e7cee11, 0x3fa62e4a00000000),
    (0xbe51f3a8c6c95003, 0x3fabb9dc00000000),
    (0xbe5774cd4fb8c30d, 0x3fb0a2b200000000),
    (0x3e2452e56c030a0a, 0x3fb3687f00000000),
    (0x3e36b63c4966a79a, 0x3fb62e4100000000),
    (0xbe3b20a21ccb525e, 0x3fb8f40a00000000),
    (0x3e54006cfb3d8f85, 0x3fbbb9d100000000),
    (0xbe5cdb026b310c41, 0x3fbe7f9b00000000),
    (0xbe569124fdc0f16d, 0x3fc0a2b080000000),
    (0xbe5084656cdc2727, 0x3fc2059580000000),
    (0xbe5376fa8b0357fd, 0x3fc3687c00000000),
    (0x3e3e56ae55a47b4a, 0x3fc4cb5e80000000),
    (0x3e5070ff8834eeb4, 0x3fc62e4400000000),
    (0x3e5623516109f4fe, 0x3fc7912900000000),
    (0xbe2ec656b95fbdac, 0x3fc8f40b00000000),
    (0x3e3f0ca2e729f510, 0x3fca56ed80000000),
    (0xbe57d260a858354a, 0x3fcbb9d680000000),
    (0x3e4e7279075503d3, 0x3fcd1cb900000000),
    (0x3e439e1a0a503873, 0x3fce7f9d00000000),
    (0x3e5cd86d7b87c3d6, 0x3fcfe27d80000000),
    (0x3e5060ab88de341e, 0x3fd0a2b240000000),
    (0x3e320a860d3f9390, 0x3fd1542440000000),
    (0xbe4dacee95fc2f10, 0x3fd2059740000000),
    (0x3e545de3a86e0aca, 0x3fd2b70700000000),
    (0x3e4c164cbfb991af, 0x3fd3687b00000000),
    (0x3e5d3f66b24225ef, 0x3fd419ec40000000),
    (0x3e5fc023efa144ba, 0x3fd4cb5f80000000),
    (0x3e3086a8af6f26c0, 0x3fd57cd280000000),
    (0xbe105c610ca86c39, 0x3fd62e4300000000),
];

static ATANH_L2: [(u64, u64); 33] = [
    (0x0000000000000000, 0x0000000000000000),
    (0xbe337e152a129e4e, 0x3f36320000000000),
    (0xbe53f6c916b8be9c, 0x3f46300000000000),
    (0x3e520505936739d5, 0x3f50a24000000000),
    (0xbe523e2e8cb541ba, 0x3f562dc000000000),
    (0xbdfacb7983ac4f5e, 0x3f5bba0000000000),
    (0x3e36f7c7689c63ae, 0x3f60a2a000000000),
    (0x3e1f5ca695b4c58b, 0x3f6368c000000000),
    (0xbe4c6c18bd953226, 0x3f662e6000000000),
    (0x3e57a516c34846bd, 0x3f68f46000000000),
    (0xbe4f3b83dd8b8530, 0x3f6bba0000000000),
    (0xbe0c3459046e4e57, 0x3f6e800000000000),
    (0x3d9b5c7e34cb79f6, 0x3f70a2c000000000),
    (0xbe42487e9af9a692, 0x3f7205c000000000),
    (0x3e5f21bbc4ad79ce, 0x3f73687000000000),
    (0xbe2550ffc857b731, 0x3f74cb7000000000),
    (0x3e487458ec1b7b34, 0x3f762e2000000000),
    (0x3e5103d4fe83ee81, 0x3f77911000000000),
    (0x3e4810483d3b398c, 0x3f78f44000000000),
    (0xbe42085cb340608e, 0x3f7a573000000000),
    (0x3e512698a119c42f, 0x3f7bb9d000000000),
    (0xbe5edb8c172b4c33, 0x3f7d1cc000000000),
    (0xbe58b55b87a5e238, 0x3f7e7fe000000000),
    (0x3e5be5e17763f78a, 0x3f7fe2b000000000),
    (0xbe1c2d496790073e, 0x3f80a2a800000000),
    (0x3e56542f523abeec, 0x3f81541000000000),
    (0xbe5b7fdbe5b193f8, 0x3f8205a000000000),
    (0x3e5fa4d42fe30c7c, 0x3f82b70000000000),
    (0x3e50d46ad04adc86, 0x3f83688800000000),
    (0xbe51c22d02d17c4c, 0x3f8419f000000000),
    (0x3e1a7d1e330dccce, 0x3f84cb7000000000),
    (0x3e0187025e656ba3, 0x3f857cd000000000),
    (0xbe4532c1269e2038, 0x3f862e5000000000),
];

#[cold]
fn as_atanh_zero(x: f64) -> f64 {
    static CH: [(u64, u64); 13] = [
        (0x3c75555555555555, 0x3fd5555555555555),
        (0xbc6999999999611c, 0x3fc999999999999a),
        (0x3c62492490f76b25, 0x3fc2492492492492),
        (0x3c5c71cd5c38a112, 0x3fbc71c71c71c71c),
        (0xbc47556c4165f4ca, 0x3fb745d1745d1746),
        (0xbc4b893c3b36052e, 0x3fb3b13b13b13b14),
        (0x3c44e1afd723ed1f, 0x3fb1111111111105),
        (0xbc4f86ea96fb1435, 0x3fae1e1e1e1e2678),
        (0x3c31e51a6e54fde9, 0x3faaf286bc9f90cc),
        (0xbc2ab913de95c3bf, 0x3fa8618618c779b6),
        (0x3c4632e747641b12, 0x3fa642c84aa383eb),
        (0xbc30c9617e7bcff2, 0x3fa47ae2d205013c),
        (0x3c23adb3e2b7f35e, 0x3fa2f664d60473f9),
    ];

    const CL: [u64; 5] = [
        0x3fa1a9a91fd692af,
        0x3fa06dfbb35e7f44,
        0x3fa037bed4d7588f,
        0x3f95aca6d6d720d6,
        0x3fa99ea5700d53a5,
    ];

    let dx2 = DoubleDouble::from_exact_mult(x, x);

    let yw0 = f_fmla(dx2.hi, f64::from_bits(CL[4]), f64::from_bits(CL[3]));
    let yw1 = f_fmla(dx2.hi, yw0, f64::from_bits(CL[2]));
    let yw2 = f_fmla(dx2.hi, yw1, f64::from_bits(CL[1]));

    let y2 = dx2.hi * f_fmla(dx2.hi, yw2, f64::from_bits(CL[0]));

    let mut y1 = lpoly_xd_generic(dx2, CH, y2);
    y1 = DoubleDouble::quick_mult_f64(y1, x);
    y1 = DoubleDouble::quick_mult(y1, dx2);

    let y0 = DoubleDouble::from_exact_add(x, y1.hi);
    let mut p = DoubleDouble::from_exact_add(y0.lo, y1.lo);

    let mut t = p.hi.to_bits();
    if (t & 0x000fffffffffffff) == 0 {
        let w = p.lo.to_bits();
        if ((w ^ t) >> 63) != 0 {
            t = t.wrapping_sub(1);
        } else {
            t = t.wrapping_add(1);
        }
        p.hi = f64::from_bits(t);
    }
    y0.hi + p.hi
}

#[cold]
fn atanh_refine(x: f64, a: f64, z: DoubleDouble) -> f64 {
    let mut t = z.hi.to_bits();
    let ex: i32 = (t >> 52) as i32;
    let e = ex.wrapping_sub(0x3ff);
    t &= 0x000fffffffffffff;
    t |= 0x3ffu64 << 52;
    let ed = e as f64;
    let v = (a - ed + f64::from_bits(0x3ff0000800000000)).to_bits();
    let i = (v.wrapping_sub(0x3ffu64 << 52)) >> (52 - 16);
    let i1: i32 = ((i >> 12) as i32) & 0x1f;
    let i2 = (i >> 8) & 0xf;
    let i3 = (i >> 4) & 0xf;
    let i4 = i & 0xf;

    const L20: f64 = f64::from_bits(0x3fd62e42fefa3a00);
    const L21: f64 = f64::from_bits(0xbcd0ca86c3898d00);
    const L22: f64 = f64::from_bits(0x397f97b57a079a00);

    let el2 = L22 * ed;
    let el1 = L21 * ed;
    let el0 = L20 * ed;

    let ll0i1 = ACOSH_ASINH_LL[0][i1 as usize];
    let ll1i2 = ACOSH_ASINH_LL[1][i2 as usize];
    let ll2i3 = ACOSH_ASINH_LL[2][i3 as usize];
    let ll3i4 = ACOSH_ASINH_LL[3][i4 as usize];

    let mut dl0 = f64::from_bits(ll0i1.0)
        + f64::from_bits(ll1i2.0)
        + (f64::from_bits(ll2i3.0) + f64::from_bits(ll3i4.0));
    let dl1 = f64::from_bits(ll0i1.1)
        + f64::from_bits(ll1i2.1)
        + (f64::from_bits(ll2i3.1) + f64::from_bits(ll3i4.1));
    let dl2 = f64::from_bits(ll0i1.2)
        + f64::from_bits(ll1i2.2)
        + (f64::from_bits(ll2i3.2) + f64::from_bits(ll3i4.2));
    dl0 += el0;

    let t12 = f64::from_bits(ACOSH_SINH_REFINE_T1[i1 as usize])
        * f64::from_bits(ACOSH_ASINH_REFINE_T2[i2 as usize]);
    let t34 = f64::from_bits(ACOSH_SINH_REFINE_T3[i3 as usize])
        * f64::from_bits(ACOSH_ASINH_REFINE_T4[i4 as usize]);
    let th = t12 * t34;
    let tl = f_fmla(t12, t34, -th);
    let dh = th * f64::from_bits(t);
    let dl = f_fmla(th, f64::from_bits(t), -dh);
    let sh = tl * f64::from_bits(t);
    let sl = f_fmla(tl, f64::from_bits(t), -sh);
    let mut dx = DoubleDouble::from_exact_add(dh - 1., dl);

    t = z.lo.to_bits();
    t = t.wrapping_sub(((e as i64) << 52) as u64);
    dx.lo += th * f64::from_bits(t);

    dx = DoubleDouble::add(dx, DoubleDouble::new(sl, sh));
    const CL: [u64; 3] = [0xbfc0000000000000, 0x3fb9999999a0754f, 0xbfb55555555c3157];

    let slw0 = f_fmla(dx.hi, f64::from_bits(CL[2]), f64::from_bits(CL[1]));

    let sl = dx.hi * f_fmla(dx.hi, slw0, f64::from_bits(CL[0]));

    static CH: [(u64, u64); 3] = [
        (0x39024b67ee516e3b, 0x3fe0000000000000),
        (0xb91932ce43199a8d, 0xbfd0000000000000),
        (0x3c655540c15cf91f, 0x3fc5555555555555),
    ];

    let mut s = lpoly_xd_generic(dx, CH, sl);
    s = DoubleDouble::mult(dx, s);
    s = DoubleDouble::add(s, DoubleDouble::new(el2, el1));
    s = DoubleDouble::add(s, DoubleDouble::new(dl2, dl1));
    let mut v02 = DoubleDouble::from_exact_add(dl0, s.hi);
    let mut v12 = DoubleDouble::from_exact_add(v02.lo, s.lo);
    t = v12.hi.to_bits();
    if (t & 0x000fffffffffffff) == 0 {
        let w = v12.lo.to_bits();
        if ((w ^ t) >> 63) != 0 {
            t = t.wrapping_sub(1);
        } else {
            t = t.wrapping_add(1);
        }
        v12.hi = f64::from_bits(t);
    }
    v02.hi *= f64::copysign(1., x);
    v12.hi *= f64::copysign(1., x);

    v02.hi + v12.hi
}

/// Hyperbolic arc tangent
///
/// Max ULP 0.5
pub fn f_atanh(x: f64) -> f64 {
    let ax = x.abs();
    let ix = ax.to_bits();
    let aix = ix;
    if aix >= 0x3ff0000000000000u64 {
        // |x| >= 1
        if aix == 0x3ff0000000000000u64 {
            // |x| = 1
            return f64::copysign(1., x) / 0.0;
        }
        if aix > 0x7ff0000000000000u64 {
            return x + x;
        } // nan
        return (-1.0f64).sqrt();
    }

    if aix < 0x3fd0000000000000u64 {
        // |x| < 0x1p-2
        // atanh(x) rounds to x to nearest for |x| < 0x1.d12ed0af1a27fp-27
        if aix < 0x3e4d12ed0af1a27fu64 {
            // |x| < 0x1.d12ed0af1a27fp-27
            /* We have underflow exactly when 0 < |x| < 2^-1022:
            for RNDU, atanh(2^-1022-2^-1074) would round to 2^-1022-2^-1075
            with unbounded exponent range */
            return f_fmla(x, f64::from_bits(0x3c80000000000000), x);
        }
        let x2 = x * x;
        const C: [u64; 9] = [
            0x3fc999999999999a,
            0x3fc2492492492244,
            0x3fbc71c71c79715f,
            0x3fb745d16f777723,
            0x3fb3b13ca4174634,
            0x3fb110c9724989bd,
            0x3fae2d17608a5b2e,
            0x3faa0b56308cba0b,
            0x3fafb6341208ad2e,
        ];
        let dx2: f64 = dd_fmla(x, x, -x2);

        let x4 = x2 * x2;
        let x3 = x2 * x;
        let x8 = x4 * x4;

        let zdx3: f64 = dd_fmla(x2, x, -x3);

        let dx3 = f_fmla(dx2, x, zdx3);

        let pw0 = f_fmla(x2, f64::from_bits(C[7]), f64::from_bits(C[6]));
        let pw1 = f_fmla(x2, f64::from_bits(C[5]), f64::from_bits(C[4]));
        let pw2 = f_fmla(x2, f64::from_bits(C[3]), f64::from_bits(C[2]));
        let pw3 = f_fmla(x2, f64::from_bits(C[1]), f64::from_bits(C[0]));

        let pw4 = f_fmla(x4, pw0, pw1);
        let pw5 = f_fmla(x8, f64::from_bits(C[8]), pw4);
        let pw6 = f_fmla(x4, pw2, pw3);

        let p = f_fmla(x8, pw5, pw6);
        let t = f64::from_bits(0x3c75555555555555) + x2 * p;
        let mut p = DoubleDouble::from_exact_add(f64::from_bits(0x3fd5555555555555), t);
        p = DoubleDouble::mult(p, DoubleDouble::new(dx3, x3));
        let z0 = DoubleDouble::from_exact_add(x, p.hi);
        p.lo += z0.lo;
        p.hi = z0.hi;
        let eps = x * f_fmla(
            x4,
            f64::from_bits(0x3cad000000000000),
            f64::from_bits(0x3980000000000000),
        );
        let lb = p.hi + (p.lo - eps);
        let ub = p.hi + (p.lo + eps);
        if lb == ub {
            return lb;
        }
        return as_atanh_zero(x);
    }

    let p = DoubleDouble::from_exact_add(1.0, ax);
    let q = DoubleDouble::from_exact_sub(1.0, ax);
    let iqh = 1.0 / q.hi;
    let th = p.hi * iqh;
    let tl = dd_fmla(p.hi, iqh, -th) + (p.lo + p.hi * (dd_fmla(-q.hi, iqh, 1.) - q.lo * iqh)) * iqh;

    const C: [u64; 5] = [
        0xbff0000000000000,
        0x3ff5555555555530,
        0xbfffffffffffffa0,
        0x40099999e33a6366,
        0xc01555559ef9525f,
    ];

    let mut t = th.to_bits();
    let ex: i32 = (t >> 52) as i32;
    let e = ex.wrapping_sub(0x3ff);
    t &= 0x000fffffffffffff;
    let ed = e as f64;
    let i = t >> (52 - 5);
    let d: i64 = (t & 0x00007fffffffffff) as i64;
    let b_i = ACOSH_ASINH_B[i as usize];
    let j: u64 = t
        .wrapping_add((b_i[0] as u64).wrapping_shl(33))
        .wrapping_add((b_i[1] as i64).wrapping_mul(d >> 16) as u64)
        >> (52 - 10);
    t |= 0x3ffu64 << 52;
    let i1: i32 = (j >> 5) as i32;
    let i2 = j & 0x1f;
    let r = (0.5 * f64::from_bits(ACOSH_ASINH_R1[i1 as usize]))
        * f64::from_bits(ACOSH_ASINH_R2[i2 as usize]);
    let dx = dd_fmla(r, f64::from_bits(t), -0.5);
    let dx2 = dx * dx;
    let rx = r * f64::from_bits(t);
    let dxl = dd_fmla(r, f64::from_bits(t), -rx);

    let fw0 = f_fmla(dx, f64::from_bits(C[3]), f64::from_bits(C[2]));
    let fw2 = f_fmla(dx, f64::from_bits(C[1]), f64::from_bits(C[0]));
    let fw1 = f_fmla(dx2, f64::from_bits(C[4]), fw0);

    let f = dx2 * f_fmla(dx2, fw1, fw2);
    const L2H: f64 = f64::from_bits(0x3fd62e42fefa3a00);
    const L2L: f64 = f64::from_bits(0xbcd0ca86c3898d00);
    let l1i1 = ATANH_L1[i1 as usize];
    let l1i2 = ATANH_L2[i2 as usize];
    let lh = f_fmla(L2H, ed, f64::from_bits(l1i1.1) + f64::from_bits(l1i2.1));
    let mut zl = DoubleDouble::from_exact_add(lh, rx - 0.5);

    zl.lo += f_fmla(L2L, ed, f64::from_bits(l1i1.0) + f64::from_bits(l1i2.0)) + dxl + 0.5 * tl / th;
    zl.lo += f;
    zl.hi *= f64::copysign(1., x);
    zl.lo *= f64::copysign(1., x);
    let eps = 31e-24 + dx2 * f64::from_bits(0x3ce0000000000000);
    let lb = zl.hi + (zl.lo - eps);
    let ub = zl.hi + (zl.lo + eps);
    if lb == ub {
        return lb;
    }
    let t = DoubleDouble::from_exact_add(th, tl);
    atanh_refine(
        x,
        f64::from_bits(0x40071547652b82fe) * (zl.hi + zl.lo).abs(),
        t,
    )
}

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

    #[test]
    fn test_atanh() {
        assert_eq!(f_atanh(-0.5000824928283691), -0.5494161408216048);
        assert_eq!(f_atanh(-0.24218760943040252), -0.24709672810738792);
    }
}