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::sincospi::f_fast_sinpi_dd;

// Generated in Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},11,11},WorkingPrecision->75,MaxIterations->100]
// num=Numerator[approx][[1]];
// poly=num;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static P_1: [(u64, u64); 12] = [
    (0x3c81873d89121fec, 0x3ffa51a6625307d3),
    (0x3cb78bf7af507504, 0x4016a65cbac476ca),
    (0xbc94179c65e021c2, 0x40218234a0a79582),
    (0x3cb842a8ab5e0994, 0x401fde32175f8515),
    (0xbc8768b33f5776b7, 0x4012de6bde49abff),
    (0x3c8e06354a27f081, 0x3ffe5d6ef3a2eac6),
    (0xbc8b391d09fed17a, 0x3fe0d186688252cf),
    (0xbc59a5c46bb8b8cc, 0x3fb958f9a0f156b7),
    (0x3c2ac44c6a197244, 0x3f88e605f24e1a89),
    (0x3bd2d05fa8be27f2, 0x3f4cd369f5d68104),
    (0x3b84ad0a748fdd22, 0x3efde955ebb17874),
    (0xb96aa8b9a65e0899, 0xbce053d04459ead7),
];

// Generated by Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx}=MiniMaxApproximation[f[z],{z,{0,3},11,11},WorkingPrecision->75,MaxIterations->100]
// num=Numerator[approx][[1]];
// poly=num;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static P_2: [(u64, u64); 12] = [
    (0x3c81873d8912236c, 0x3ffa51a6625307d3),
    (0xbcbc731a62342288, 0x40176c23f7ea51e6),
    (0xbcc45a6fd00e67a8, 0x4022cb5ae67657ef),
    (0x3cc0876fde7fe4e6, 0x4021d766062b9550),
    (0xbcaec4a4859cba1d, 0x401629f91cd4f291),
    (0x3c76184014e4d7e3, 0x4002d43da3352004),
    (0x3c812c7609483e0e, 0x3fe62e3266eef8c7),
    (0xbc5f991047f52d2b, 0x3fc1eacb910b951c),
    (0x3c28b9f38d603f2f, 0x3f930960a301df34),
    (0x3bf9b620eb930504, 0x3f5814f8e057b14b),
    (0xbb990860b88b54e4, 0x3f0b9f67c71aa3bf),
    (0x38e5cb6acfbaab77, 0xbc4194b8c01afe9a),
];

// Generated by Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{3,9},11,11},WorkingPrecision->75,MaxIterations->100]
// num=Numerator[approx];
// poly=num;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static P_3: [(u64, u64); 12] = [
    (0x3c9cd56dbb295efc, 0x3ffa51a662556db9),
    (0x3c9f4ee74f5f9daf, 0x4018ff913088cb34),
    (0x3ccf08737350609c, 0x402593d55686b8b1),
    (0xbcc6cd4ed33afebb, 0x402641d10de4def5),
    (0xbcb24d1957c1303c, 0x401e682c37e8e2cf),
    (0x3ca30ac79162ceb2, 0x400ccfc7c4566f55),
    (0x3c9efea5ff293dc9, 0x3ff33eb2c6e89d0b),
    (0x3c74670a11068abc, 0x3fd1fbf456e5c6f0),
    (0x3c47b5dcdea19c36, 0x3fa6a6a2148c482c),
    (0xbc14642012a1cc1e, 0x3f71851e927f52e7),
    (0x3bc7db88a4ec5478, 0x3f29a45059a43475),
    (0xb7bc31e55271eab0, 0xbb375518529c52fb),
];

// Generated in Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},11,11},WorkingPrecision->75,MaxIterations->100]
// den=Denominator[approx][[1]];
// poly=den;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static Q_1: [(u64, u64); 12] = [
    (0x0000000000000000, 0x3ff0000000000000),
    (0xbcb84c43a11fc28a, 0x40139d9587da0fb5),
    (0x3ca1cf3dbcddbb57, 0x402507cb6225f0f0),
    (0x3cb01aa6ddcc3cfd, 0x402a1b416d0ed4e6),
    (0xbcbc31c216b5ff66, 0x4024ec8829e535d4),
    (0xbcb335c23022f43e, 0x4016d2ba6d1a18e6),
    (0x3cafbfffc03ad28a, 0x400158c4611ed51f),
    (0xbc8d1fb10a031a27, 0x3fe26f15bb52f89b),
    (0x3c56a9fea160eecb, 0x3fbaec13f663049d),
    (0xbc2f4ee869ba9364, 0x3f89cde0500cd68f),
    (0x3be3f23afc9398b6, 0x3f4d4b0f4dcf3eb8),
    (0x3b67a3ed4795d33e, 0x3efde955eafac9c2),
];

// Generated by Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx}=MiniMaxApproximation[f[z],{z,{0,3},11,11},WorkingPrecision->75,MaxIterations->100]
// den=Denominator[approx][[1]];
// poly=den;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static Q_2: [(u64, u64); 12] = [
    (0x0000000000000000, 0x3ff0000000000000),
    (0xbc81a3e4e026b7b1, 0x401415d1a20a9339),
    (0x3cc279576dfe3ec9, 0x402627c1a95d33d2),
    (0x3c9a94b5cf0cae88, 0x402c724dc5cf4577),
    (0xbc8aa1fa0c3820a8, 0x4027b7a332bb07f4),
    (0xbc96968367088d66, 0x401b14376177bdd7),
    (0x3ca2d3dfa5847f4d, 0x4005b0511cd98f2c),
    (0xbc8cfad394d41dd1, 0x3fe877bc2d02c7f3),
    (0xbc51592b8ec81a92, 0x3fc31f52afc72b95),
    (0x3c2cbef277d587e9, 0x3f93cb2f0e574376),
    (0xbbfbb670fd94f6ba, 0x3f5883767f745a92),
    (0xbb931b04d74e5893, 0x3f0b9f67c71a60f3),
];

// Generated by Wolfram Mathematica:
// <<FunctionApproximations`
// ClearAll["Global`*"]
// f[x_]:=PolyGamma[1, x+1]
// {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{3,9},11,11},WorkingPrecision->75,MaxIterations->100]
// den=Denominator[approx];
// poly=den;
// coeffs=CoefficientList[poly,z];
// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
static Q_3: [(u64, u64); 12] = [
    (0x0000000000000000, 0x3ff0000000000000),
    (0x3cbafcb4b6d646d9, 0x40150b12a79fc9cf),
    (0x3c989ef814b8dd2a, 0x40288c1d26ffdca5),
    (0x3cb0282cfea9c473, 0x4030d737c893cd5f),
    (0x3cc955b8aaadb37d, 0x402e5c289b6de3e0),
    (0x3cb377161f8861d2, 0x4022fb66d87bd522),
    (0xbcb4b0e4cff46ad6, 0x4010e5d13c2a5907),
    (0xbc8824539e4b1bd6, 0x3ff58c8fe8f26fca),
    (0xbc7d34220d810ea0, 0x3fd36c1351f43e66),
    (0xbc4cbdbe85570017, 0x3fa7c1170466605e),
    (0xbc0c3afb98775c53, 0x3f71ebafd3e5e3b9),
    (0x3bc0b0b7f16afd0a, 0x3f29a45059a43475),
];

#[inline]
fn approx_trigamma(x: f64) -> DoubleDouble {
    if x <= 10. {
        let (p, q) = if x <= 1. {
            (&P_1, &Q_1)
        } else if x <= 4. {
            (&P_2, &Q_2)
        } else {
            (&P_3, &Q_3)
        };
        let x2 = DoubleDouble::from_exact_mult(x, x);
        let x4 = DoubleDouble::quick_mult(x2, x2);
        let x8 = DoubleDouble::quick_mult(x4, x4);

        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 rcp2 = DoubleDouble::quick_mult(rcp, rcp);

        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_add(q, rcp2);
        return r;
    }
    // asymptotic expansion Trigamma[x] = 1/x + 1/x^2 + sum(Bernoulli(2*k)/x^(2*k + 1))
    // Generated in SageMath:
    // var('x')
    // def bernoulli_terms(x, N):
    //     S = 0
    //     for k in range(1, N+1):
    //         B = bernoulli(2*k)
    //         term = B*x**(-(2*k+1))
    //         S += term
    //     return S
    //
    // terms = bernoulli_terms(x, 10)
    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
    // for k in range(0, 14):
    //     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] = [
        (0x3c65555555555555, 0x3fc5555555555555),
        (0xbc21111111111111, 0xbfa1111111111111),
        (0x3c38618618618618, 0x3f98618618618618),
        (0xbc21111111111111, 0xbfa1111111111111),
        (0xbc4364d9364d9365, 0x3fb364d9364d9365),
        (0xbc6981981981981a, 0xbfd0330330330330),
        (0xbc95555555555555, 0x3ff2aaaaaaaaaaab),
        (0xbcb7979797979798, 0xc01c5e5e5e5e5e5e),
        (0xbcac3b070ec1c3b0, 0x404b7c4f8f13e3c5),
        (0x3cc8d3018d3018d3, 0xc08088fe72cfe72d),
    ];

    let rcp = DoubleDouble::from_quick_recip(x);

    let q = DoubleDouble::quick_mult(rcp, rcp);

    let q2 = DoubleDouble::quick_mult(q, q);
    let q4 = q2 * q2;
    let q8 = q4 * q4;

    let e0 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[1]),
        q,
        DoubleDouble::from_bit_pair(C[0]),
    );
    let e1 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[3]),
        q,
        DoubleDouble::from_bit_pair(C[2]),
    );
    let e2 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[5]),
        q,
        DoubleDouble::from_bit_pair(C[4]),
    );
    let e3 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[7]),
        q,
        DoubleDouble::from_bit_pair(C[6]),
    );
    let e4 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[9]),
        q,
        DoubleDouble::from_bit_pair(C[8]),
    );

    let q0 = DoubleDouble::quick_mul_add(q2, e1, e0);
    let q1 = DoubleDouble::quick_mul_add(q2, e3, e2);

    let r0 = DoubleDouble::quick_mul_add(q4, q1, q0);
    let mut p = DoubleDouble::quick_mul_add(q8, e4, r0);

    let q_over_2 = DoubleDouble::quick_mult_f64(q, 0.5);
    p = DoubleDouble::quick_mult(p, q);
    p = DoubleDouble::quick_mult(p, rcp);
    p = DoubleDouble::quick_dd_add(q_over_2, p);
    p = DoubleDouble::quick_dd_add(p, rcp);
    p
}

#[inline]
fn approx_trigamma_dd(x: DoubleDouble) -> DoubleDouble {
    if x.hi <= 10. {
        let (p, q) = if x.hi <= 1. {
            (&P_1, &Q_1)
        } else if x.hi <= 4. {
            (&P_2, &Q_2)
        } else {
            (&P_3, &Q_3)
        };
        let x2 = DoubleDouble::quick_mult(x, x);
        let x4 = DoubleDouble::quick_mult(x2, x2);
        let x8 = DoubleDouble::quick_mult(x4, x4);

        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 rcp2 = DoubleDouble::quick_mult(rcp, rcp);

        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_add(q, rcp2);
        return r;
    }
    // asymptotic expansion Trigamma[x] = 1/x + 1/x^2 + sum(Bernoulli(2*k)/x^(2*k + 1))
    // Generated in SageMath:
    // var('x')
    // def bernoulli_terms(x, N):
    //     S = 0
    //     for k in range(1, N+1):
    //         B = bernoulli(2*k)
    //         term = B*x**(-(2*k+1))
    //         S += term
    //     return S
    //
    // terms = bernoulli_terms(x, 10)
    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
    // for k in range(0, 14):
    //     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] = [
        (0x3c65555555555555, 0x3fc5555555555555),
        (0xbc21111111111111, 0xbfa1111111111111),
        (0x3c38618618618618, 0x3f98618618618618),
        (0xbc21111111111111, 0xbfa1111111111111),
        (0xbc4364d9364d9365, 0x3fb364d9364d9365),
        (0xbc6981981981981a, 0xbfd0330330330330),
        (0xbc95555555555555, 0x3ff2aaaaaaaaaaab),
        (0xbcb7979797979798, 0xc01c5e5e5e5e5e5e),
        (0xbcac3b070ec1c3b0, 0x404b7c4f8f13e3c5),
        (0x3cc8d3018d3018d3, 0xc08088fe72cfe72d),
    ];

    let rcp = x.recip();

    let q = DoubleDouble::quick_mult(rcp, rcp);

    let q2 = DoubleDouble::quick_mult(q, q);
    let q4 = q2 * q2;
    let q8 = q4 * q4;

    let e0 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[1]),
        q,
        DoubleDouble::from_bit_pair(C[0]),
    );
    let e1 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[3]),
        q,
        DoubleDouble::from_bit_pair(C[2]),
    );
    let e2 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[5]),
        q,
        DoubleDouble::from_bit_pair(C[4]),
    );
    let e3 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[7]),
        q,
        DoubleDouble::from_bit_pair(C[6]),
    );
    let e4 = DoubleDouble::quick_mul_add(
        DoubleDouble::from_bit_pair(C[9]),
        q,
        DoubleDouble::from_bit_pair(C[8]),
    );

    let q0 = DoubleDouble::quick_mul_add(q2, e1, e0);
    let q1 = DoubleDouble::quick_mul_add(q2, e3, e2);

    let r0 = DoubleDouble::quick_mul_add(q4, q1, q0);
    let mut p = DoubleDouble::quick_mul_add(q8, e4, r0);

    let q_over_2 = DoubleDouble::quick_mult_f64(q, 0.5);
    p = DoubleDouble::quick_mult(p, q);
    p = DoubleDouble::quick_mult(p, rcp);
    p = DoubleDouble::quick_dd_add(q_over_2, p);
    p = DoubleDouble::quick_dd_add(p, rcp);
    p
}

/// Computes the trigamma function ψ₁(x).
///
/// The trigamma function is the second derivative of the logarithm of the gamma function.
pub fn f_trigamma(x: f64) -> f64 {
    let xb = x.to_bits();
    if !x.is_normal() {
        if x.is_infinite() {
            return if x.is_sign_negative() {
                f64::NEG_INFINITY
            } else {
                0.
            };
        }
        if x.is_nan() {
            return f64::NAN;
        }
        if xb == 0 {
            return f64::INFINITY;
        }
    }

    let x_e = (x.to_bits() >> 52) & 0x7ff;

    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;

    if x_e < E_BIAS - 52 {
        // |x| < 2^-52
        let dx = x;
        return 1. / (dx * dx);
    }

    if x < 0. {
        if is_integer(x) {
            return f64::INFINITY;
        }
        // reflection formula
        // Trigamma[1-x] + Trigamma[x] = PI^2 / sinpi^2(x)
        const SQR_PI: DoubleDouble =
            DoubleDouble::from_bit_pair((0x3cc692b71366cc05, 0x4023bd3cc9be45de)); // pi^2
        let sinpi_ax = f_fast_sinpi_dd(-x);
        let dx = DoubleDouble::from_full_exact_sub(1., x);
        let result = DoubleDouble::div(SQR_PI, DoubleDouble::quick_mult(sinpi_ax, sinpi_ax));
        let trigamma_x = approx_trigamma_dd(dx);
        return DoubleDouble::quick_dd_sub(result, trigamma_x).to_f64();
    }

    approx_trigamma(x).to_f64()
}

#[cfg(test)]
mod tests {

    use super::*;

    #[test]
    fn test_trigamma() {
        assert_eq!(f_trigamma(-27.058018), 300.35629698636757);
        assert_eq!(f_trigamma(27.058018), 0.037648965757704725);
        assert_eq!(f_trigamma(8.058018), 0.13211796975281037);
        assert_eq!(f_trigamma(-8.058018), 300.2758629255111);
        assert_eq!(f_trigamma(2.23432), 0.5621320243666134);
        assert_eq!(f_trigamma(-2.4653), 9.653674003034206);
        assert_eq!(f_trigamma(0.123541), 66.91128231455282);
        assert_eq!(f_trigamma(-0.54331), 9.154415950366596);
        assert_eq!(f_trigamma(-5.), f64::INFINITY);
        assert_eq!(f_trigamma(f64::INFINITY), 0.0);
        assert_eq!(f_trigamma(f64::NEG_INFINITY), f64::NEG_INFINITY);
        assert!(f_trigamma(f64::NAN).is_nan());
    }
}