pxfm 0.1.29

/*
 * // Copyright (c) Radzivon Bartoshyk 8/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::is_integer;
use crate::double_double::DoubleDouble;
use crate::gamma::digamma_coeffs::digamma_coeefs;
use crate::logs::fast_log_d_to_dd;
use crate::tangent::cotpi_core;

#[inline]
fn approx_digamma_hard(x: f64) -> DoubleDouble {
    if x <= 12. {
        let x2 = DoubleDouble::from_exact_mult(x, x);
        let x4 = DoubleDouble::quick_mult(x2, x2);
        let x8 = DoubleDouble::quick_mult(x4, x4);
        let (p, q) = digamma_coeefs(x);

        let e0 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[1]),
            x,
            DoubleDouble::from_bit_pair(p[0]),
        );
        let e1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[3]),
            x,
            DoubleDouble::from_bit_pair(p[2]),
        );
        let e2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[5]),
            x,
            DoubleDouble::from_bit_pair(p[4]),
        );
        let e3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[7]),
            x,
            DoubleDouble::from_bit_pair(p[6]),
        );
        let e4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[9]),
            x,
            DoubleDouble::from_bit_pair(p[8]),
        );
        let e5 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(p[11]),
            x,
            DoubleDouble::from_bit_pair(p[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_num = DoubleDouble::mul_add(x8, f2, g0);

        let rcp = DoubleDouble::from_quick_recip(x);

        let e0 = DoubleDouble::mul_f64_add_f64(
            DoubleDouble::from_bit_pair(q[1]),
            x,
            f64::from_bits(0x3ff0000000000000),
        );
        let e1 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(q[3]),
            x,
            DoubleDouble::from_bit_pair(q[2]),
        );
        let e2 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(q[5]),
            x,
            DoubleDouble::from_bit_pair(q[4]),
        );
        let e3 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(q[7]),
            x,
            DoubleDouble::from_bit_pair(q[6]),
        );
        let e4 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(q[9]),
            x,
            DoubleDouble::from_bit_pair(q[8]),
        );
        let e5 = DoubleDouble::mul_f64_add(
            DoubleDouble::from_bit_pair(q[11]),
            x,
            DoubleDouble::from_bit_pair(q[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_den = DoubleDouble::mul_add(x8, f2, g0);

        let q = DoubleDouble::div(p_num, p_den);
        let r = DoubleDouble::quick_dd_sub(q, rcp);
        return r;
    }
    // digamma asymptotic expansion
    // digamma(x) ~ ln(z)+1/(2z)-sum_(n=1)^(infty)(Bernoulli(2n))/(2nz^(2n))
    // Generated in SageMath:
    // var('x')
    //
    // def bernoulli_terms(x, N):
    //     S = 0
    //     S += QQ(1)/QQ(2)/x
    //     for k in range(1, N+1):
    //         B = bernoulli(2*k)
    //         term = (B / QQ(2*k))*x**(-2*k)
    //         S += term
    //     return S
    //
    // terms = bernoulli_terms(x, 9)
    //
    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
    // for k in range(1, 13):
    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
    //     if c == 0:
    //         continue
    //     print("f64::from_bits(" + double_to_hex(c) + "),")
    const C: [(u64, u64); 10] = [
        (0x3c55555555555555, 0x3fb5555555555555),
        (0xbc01111111111111, 0xbf81111111111111),
        (0x3c10410410410410, 0x3f70410410410410),
        (0xbbf1111111111111, 0xbf71111111111111),
        (0xbc0f07c1f07c1f08, 0x3f7f07c1f07c1f08),
        (0x3c39a99a99a99a9a, 0xbf95995995995996),
        (0x3c55555555555555, 0x3fb5555555555555),
        (0xbc77979797979798, 0xbfdc5e5e5e5e5e5e),
        (0xbc69180646019180, 0x40086e7f9b9fe6e8),
        (0x3ccad759ad759ad7, 0xc03a74ca514ca515),
    ];

    let rcp = DoubleDouble::from_quick_recip(x);
    let rcp_sqr = DoubleDouble::quick_mult(rcp, rcp);

    let rx = rcp_sqr;
    let x2 = DoubleDouble::quick_mult(rx, rx);
    let x4 = DoubleDouble::quick_mult(x2, x2);
    let x8 = DoubleDouble::quick_mult(x4, x4);

    let q0 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[1]),
        rx,
        DoubleDouble::from_bit_pair(C[0]),
    );
    let q1 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[3]),
        rx,
        DoubleDouble::from_bit_pair(C[2]),
    );
    let q2 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[5]),
        rx,
        DoubleDouble::from_bit_pair(C[4]),
    );
    let q3 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[7]),
        rx,
        DoubleDouble::from_bit_pair(C[6]),
    );
    let q4 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[9]),
        rx,
        DoubleDouble::from_bit_pair(C[8]),
    );

    let q0 = DoubleDouble::quick_mul_add(x2, q1, q0);
    let q1 = DoubleDouble::quick_mul_add(x2, q3, q2);

    let r0 = DoubleDouble::quick_mul_add(x4, q1, q0);
    let mut p = DoubleDouble::quick_mul_add(x8, q4, r0);
    p = DoubleDouble::quick_mul_f64_add(
        rcp,
        f64::from_bits(0x3fe0000000000000),
        DoubleDouble::quick_mult(p, rcp_sqr),
    );

    let v_log = fast_log_d_to_dd(x);
    DoubleDouble::quick_dd_sub(v_log, p)
}

#[inline]
fn approx_digamma_hard_dd(x: DoubleDouble) -> DoubleDouble {
    if x.hi <= 12. {
        let x2 = DoubleDouble::quick_mult(x, x);
        let x4 = DoubleDouble::quick_mult(x2, x2);
        let x8 = DoubleDouble::quick_mult(x4, x4);
        let (p, q) = digamma_coeefs(x.hi);

        let e0 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[1]),
            x,
            DoubleDouble::from_bit_pair(p[0]),
        );
        let e1 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[3]),
            x,
            DoubleDouble::from_bit_pair(p[2]),
        );
        let e2 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[5]),
            x,
            DoubleDouble::from_bit_pair(p[4]),
        );
        let e3 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[7]),
            x,
            DoubleDouble::from_bit_pair(p[6]),
        );
        let e4 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[9]),
            x,
            DoubleDouble::from_bit_pair(p[8]),
        );
        let e5 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(p[11]),
            x,
            DoubleDouble::from_bit_pair(p[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_num = DoubleDouble::mul_add(x8, f2, g0);

        let rcp = x.recip();

        let e0 = DoubleDouble::mul_add_f64(
            DoubleDouble::from_bit_pair(q[1]),
            x,
            f64::from_bits(0x3ff0000000000000),
        );
        let e1 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(q[3]),
            x,
            DoubleDouble::from_bit_pair(q[2]),
        );
        let e2 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(q[5]),
            x,
            DoubleDouble::from_bit_pair(q[4]),
        );
        let e3 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(q[7]),
            x,
            DoubleDouble::from_bit_pair(q[6]),
        );
        let e4 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(q[9]),
            x,
            DoubleDouble::from_bit_pair(q[8]),
        );
        let e5 = DoubleDouble::mul_add(
            DoubleDouble::from_bit_pair(q[11]),
            x,
            DoubleDouble::from_bit_pair(q[10]),
        );

        let f0 = DoubleDouble::mul_add(x2, e1, e0);
        let f1 = DoubleDouble::mul_add(x2, e3, e2);
        let f2 = DoubleDouble::mul_add(x2, e5, e4);

        let g0 = DoubleDouble::mul_add(x4, f1, f0);

        let p_den = DoubleDouble::mul_add(x8, f2, g0);

        let q = DoubleDouble::div(p_num, p_den);
        let r = DoubleDouble::quick_dd_sub(q, rcp);
        return r;
    }
    // digamma asymptotic expansion
    // digamma(x) ~ ln(z)+1/(2z)-sum_(n=1)^(infty)(Bernoulli(2n))/(2nz^(2n))
    // Generated in SageMath:
    // var('x')
    //
    // def bernoulli_terms(x, N):
    //     S = 0
    //     S += QQ(1)/QQ(2)/x
    //     for k in range(1, N+1):
    //         B = bernoulli(2*k)
    //         term = (B / QQ(2*k))*x**(-2*k)
    //         S += term
    //     return S
    //
    // terms = bernoulli_terms(x, 9)
    //
    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
    // for k in range(1, 13):
    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
    //     if c == 0:
    //         continue
    //     print("f64::from_bits(" + double_to_hex(c) + "),")
    const C: [(u64, u64); 10] = [
        (0x3c55555555555555, 0x3fb5555555555555),
        (0xbc01111111111111, 0xbf81111111111111),
        (0x3c10410410410410, 0x3f70410410410410),
        (0xbbf1111111111111, 0xbf71111111111111),
        (0xbc0f07c1f07c1f08, 0x3f7f07c1f07c1f08),
        (0x3c39a99a99a99a9a, 0xbf95995995995996),
        (0x3c55555555555555, 0x3fb5555555555555),
        (0xbc77979797979798, 0xbfdc5e5e5e5e5e5e),
        (0xbc69180646019180, 0x40086e7f9b9fe6e8),
        (0x3ccad759ad759ad7, 0xc03a74ca514ca515),
    ];

    let rcp = x.recip();
    let rcp_sqr = DoubleDouble::quick_mult(rcp, rcp);

    let rx = rcp_sqr;
    let x2 = DoubleDouble::quick_mult(rx, rx);
    let x4 = DoubleDouble::quick_mult(x2, x2);
    let x8 = DoubleDouble::quick_mult(x4, x4);

    let q0 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[1]),
        rx,
        DoubleDouble::from_bit_pair(C[0]),
    );
    let q1 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[3]),
        rx,
        DoubleDouble::from_bit_pair(C[2]),
    );
    let q2 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[5]),
        rx,
        DoubleDouble::from_bit_pair(C[4]),
    );
    let q3 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[7]),
        rx,
        DoubleDouble::from_bit_pair(C[6]),
    );
    let q4 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[9]),
        rx,
        DoubleDouble::from_bit_pair(C[8]),
    );

    let q0 = DoubleDouble::quick_mul_add(x2, q1, q0);
    let q1 = DoubleDouble::quick_mul_add(x2, q3, q2);

    let r0 = DoubleDouble::quick_mul_add(x4, q1, q0);
    let mut p = DoubleDouble::quick_mul_add(x8, q4, r0);
    p = DoubleDouble::quick_mul_f64_add(
        rcp,
        f64::from_bits(0x3fe0000000000000),
        DoubleDouble::quick_mult(p, rcp_sqr),
    );

    let mut v_log = fast_log_d_to_dd(x.hi);
    v_log.lo += x.lo / x.hi;
    DoubleDouble::quick_dd_sub(v_log, p)
}

/// Computes digamma(x)
///
pub fn f_digamma(x: f64) -> f64 {
    let xb = x.to_bits();
    if !x.is_normal() {
        if x.is_infinite() {
            return if x.is_sign_negative() {
                f64::NAN
            } else {
                f64::INFINITY
            };
        }
        if x.is_nan() {
            return f64::NAN;
        }
        if xb.wrapping_shl(1) == 0 {
            // |x| == 0
            return f64::INFINITY;
        }
    }

    let dx = x;

    if x.abs() <= f64::EPSILON {
        // |x| < ulp(1)
        // digamma near where x -> 1 ~ Digamma[x] = -euler + O(x)
        // considering identity Digamma[x+1] = Digamma[x] + 1/x,
        // hence x < ulp(1) then x+1 ~= 1 that gives
        // Digamma[x] = Digamma[x+1] - 1/x = -euler - 1/x
        // euler is dropped since 1/x >> euler
        // that gives:
        // Digamma[x] = Digamma[x+1] - 1/x = -1/x
        return -1. / x;
    }

    if x < 0. {
        // at negative integers it's inf
        if is_integer(x) {
            return f64::NAN;
        }

        // reflection Gamma(1-x) + Gamma(x) = Pi/tan(PI*x)
        const PI: DoubleDouble =
            DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
        let r = DoubleDouble::from_full_exact_sub(1., x);
        let mut result = PI * cotpi_core(-x);
        let app = approx_digamma_hard_dd(r);
        result = DoubleDouble::quick_dd_add(result, app);
        result.to_f64()
    } else {
        let app = approx_digamma_hard(dx);
        app.to_f64()
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    #[test]
    fn test_digamma() {
        assert_eq!(f_digamma(0.0019531200000040207), -512.5753182109892);
        assert_eq!(f_digamma(-13.999000000012591), -997.3224450000563);
        assert_eq!(f_digamma(13.999000000453323), 2.602844047257257);
        assert_eq!(f_digamma(0.), f64::INFINITY);
        assert_eq!(f_digamma(-0.), f64::INFINITY);
        assert!(f_digamma(f64::NAN).is_nan());
    }
}