pxfm/gamma/
trigammaf.rs

1/*
2 * // Copyright (c) Radzivon Bartoshyk 8/2025. All rights reserved.
3 * //
4 * // Redistribution and use in source and binary forms, with or without modification,
5 * // are permitted provided that the following conditions are met:
6 * //
7 * // 1.  Redistributions of source code must retain the above copyright notice, this
8 * // list of conditions and the following disclaimer.
9 * //
10 * // 2.  Redistributions in binary form must reproduce the above copyright notice,
11 * // this list of conditions and the following disclaimer in the documentation
12 * // and/or other materials provided with the distribution.
13 * //
14 * // 3.  Neither the name of the copyright holder nor the names of its
15 * // contributors may be used to endorse or promote products derived from
16 * // this software without specific prior written permission.
17 * //
18 * // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 * // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 * // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
21 * // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
22 * // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
23 * // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
24 * // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
25 * // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
26 * // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27 * // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29use crate::common::{f_fmla, is_integerf};
30use crate::polyeval::{f_polyeval6, f_polyeval10};
31use crate::sin_cosf::fast_sinpif;
32
33#[inline]
34fn approx_trigamma(x: f64) -> f64 {
35    if x <= 1. {
36        // Polynomial for Trigamma[x+1]
37        // <<FunctionApproximations`
38        // ClearAll["Global`*"]
39        // f[x_]:=PolyGamma[1, x+1]
40        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},9,9},WorkingPrecision->75,MaxIterations->100]
41        // num=Numerator[approx];
42        // den=Denominator[approx];
43        // poly=num;
44        // coeffs=CoefficientList[poly,z];
45        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
46        // poly=den;
47        // coeffs=CoefficientList[poly,z];
48        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
49        let p_num = f_polyeval10(
50            x,
51            f64::from_bits(0x3ffa51a6625307d3),
52            f64::from_bits(0x40142fc9c4f02b3f),
53            f64::from_bits(0x401b30a762805b44),
54            f64::from_bits(0x4014dc84c95656bd),
55            f64::from_bits(0x4003e44f4c820b4c),
56            f64::from_bits(0x3fe81f37523197d3),
57            f64::from_bits(0x3fc22bffe2490221),
58            f64::from_bits(0x3f8f221a6329ea36),
59            f64::from_bits(0x3f47406930b9563c),
60            f64::from_bits(0xbd99cd44c6ad497a),
61        );
62        let p_den = f_polyeval10(
63            x,
64            f64::from_bits(0x3ff0000000000000),
65            f64::from_bits(0x40121e3db4e0a2f3),
66            f64::from_bits(0x40218e97a5430c4f),
67            f64::from_bits(0x402329897737b159),
68            f64::from_bits(0x401a0fdc27807c2d),
69            f64::from_bits(0x4006ff242e1f3a51),
70            f64::from_bits(0x3fea6eda129c4e85),
71            f64::from_bits(0x3fc32700b2ae2e88),
72            f64::from_bits(0x3f8fdc1dc6116d41),
73            f64::from_bits(0x3f4740690261cfbc),
74        );
75        return p_num / p_den + 1. / (x * x);
76    } else if x <= 4. {
77        // Polynomial for Trigamma[x+1]
78        // <<FunctionApproximations`
79        // ClearAll["Global`*"]
80        // f[x_]:=PolyGamma[1, x+1]
81        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{0,3},9,9},WorkingPrecision->75,MaxIterations->100]
82        // num=Numerator[approx];
83        // den=Denominator[approx];
84        // poly=num;
85        // coeffs=CoefficientList[poly,z];
86        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
87        // poly=den;
88        // coeffs=CoefficientList[poly,z];
89        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
90        let p_num = f_polyeval10(
91            x,
92            f64::from_bits(0x3ffa51a6625307d3),
93            f64::from_bits(0x4015167fa2d2b5a8),
94            f64::from_bits(0x401dd40865d5985e),
95            f64::from_bits(0x4018353d9425fb58),
96            f64::from_bits(0x4008a12aa45851fa),
97            f64::from_bits(0x3ff018736d0c5dbe),
98            f64::from_bits(0x3fca715702bdb519),
99            f64::from_bits(0x3f9908a9d73d983c),
100            f64::from_bits(0x3f54fd9cbbb46314),
101            f64::from_bits(0xbd00b8a28c578ab5),
102        );
103        let p_den = f_polyeval10(
104            x,
105            f64::from_bits(0x3ff0000000000000),
106            f64::from_bits(0x4012aa7f041a768b),
107            f64::from_bits(0x4022c2604e5f9c7a),
108            f64::from_bits(0x4025655b63c2db22),
109            f64::from_bits(0x401eaa8e59c8295d),
110            f64::from_bits(0x400cc8724a58809c),
111            f64::from_bits(0x3ff1c7a91c8e3c40),
112            f64::from_bits(0x3fcc05613a11183e),
113            f64::from_bits(0x3f99b096bd3ce542),
114            f64::from_bits(0x3f54fd9cbb9c6167),
115        );
116        return p_num / p_den + 1. / (x * x);
117    } else if x <= 10. {
118        // Polynomial for Trigamma[x+1]
119        // <<FunctionApproximations`
120        // ClearAll["Global`*"]
121        // f[x_]:=PolyGamma[1, x+1]
122        // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{3,9},9,9},WorkingPrecision->75,MaxIterations->100]
123        // num=Numerator[approx];
124        // den=Denominator[approx];
125        // poly=num;
126        // coeffs=CoefficientList[poly,z];
127        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
128        // poly=den;
129        // coeffs=CoefficientList[poly,z];
130        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
131        let p_num = f_polyeval10(
132            x,
133            f64::from_bits(0x3ffa51a664b1b211),
134            f64::from_bits(0x4016d7f75881312a),
135            f64::from_bits(0x4021b10defb47bcc),
136            f64::from_bits(0x401fe9633665e1bf),
137            f64::from_bits(0x4012601cce6766d7),
138            f64::from_bits(0x3ffbd0ece1c435f1),
139            f64::from_bits(0x3fdb3fd0e233c485),
140            f64::from_bits(0x3faffdedea90b870),
141            f64::from_bits(0x3f71a4bbf0d00147),
142            f64::from_bits(0xbc0aae286498a357),
143        );
144        let p_den = f_polyeval10(
145            x,
146            f64::from_bits(0x3ff0000000000000),
147            f64::from_bits(0x4013bbbd604d685d),
148            f64::from_bits(0x40253a4fca05438e),
149            f64::from_bits(0x402a4aba4634f14f),
150            f64::from_bits(0x4024cdd23bd6284a),
151            f64::from_bits(0x4015fbe371275c3f),
152            f64::from_bits(0x3fff4d7ebf7d1ed0),
153            f64::from_bits(0x3fdd459154d7bc72),
154            f64::from_bits(0x3fb08c1cd4cedca3),
155            f64::from_bits(0x3f71a4bbf0d0012d),
156        );
157        return p_num / p_den + 1. / (x * x);
158    }
159    // asymptotic expansion Trigamma[x] = 1/x + 1/x^2 + sum(Bernoulli(2*k)/x^(2*k + 1))
160    // Generated in SageMath:
161    // var('x')
162    // def bernoulli_terms(x, N):
163    //     S = 0
164    //     for k in range(1, N+1):
165    //         B = bernoulli(2*k)
166    //         term = B*x**(-(2*k+1))
167    //         S += term
168    //     return S
169    //
170    // terms = bernoulli_terms(x, 7)
171    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
172    // for k in range(0, 14):
173    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
174    //     if c == 0:
175    //         continue
176    //     print("f64::from_bits(" + double_to_hex(c) + "),")
177    let r = 1. / x;
178    let r2 = r * r;
179    let p = f_polyeval6(
180        r2,
181        f64::from_bits(0x3fc5555555555555),
182        f64::from_bits(0xbfa1111111111111),
183        f64::from_bits(0x3f98618618618618),
184        f64::from_bits(0xbfa1111111111111),
185        f64::from_bits(0x3fb364d9364d9365),
186        f64::from_bits(0xbfd0330330330330),
187    );
188    f_fmla(p, r2 * r, f_fmla(r2, 0.5, r))
189}
190
191/// Computes the trigamma function ψ₁(x).
192///
193/// The trigamma function is the second derivative of the logarithm of the gamma function.
194pub fn f_trigammaf(x: f32) -> f32 {
195    let xb = x.to_bits();
196    // filter out exceptional cases
197    if xb >= 0xffu32 << 23 || xb == 0 {
198        if x.is_infinite() {
199            return if x.is_sign_negative() {
200                f32::NEG_INFINITY
201            } else {
202                0.
203            };
204        }
205        if x.is_nan() {
206            return f32::NAN;
207        }
208        if xb.wrapping_shl(1) == 0 {
209            return f32::INFINITY;
210        }
211    }
212
213    let ax = x.to_bits() & 0x7fff_ffff;
214
215    if ax <= 0x34000000u32 {
216        // |x| < f32::EPSILON
217        let dx = x as f64;
218        return (1. / (dx * dx)) as f32;
219    }
220
221    let mut dx = x as f64;
222
223    let mut result = 0.;
224    let mut sum_parity: f64 = 1.0;
225
226    if x < 0. {
227        // singularity at negative integers
228        if is_integerf(x) {
229            return f32::INFINITY;
230        }
231        // reflection formula
232        // Trigamma[1-x] + Trigamma[x] = PI^2 / sinpi^2(x)
233        const SQR_PI: f64 = f64::from_bits(0x4023bd3cc9be45de); // pi^2
234        let sinpi_ax = fast_sinpif(-x);
235        dx = 1. - dx;
236        result = SQR_PI / (sinpi_ax * sinpi_ax);
237        sum_parity = -1.;
238    }
239
240    let r = approx_trigamma(dx) * sum_parity;
241    result += r;
242    result as f32
243}
244
245#[cfg(test)]
246mod tests {
247    use crate::f_trigammaf;
248
249    #[test]
250    fn test_trigamma() {
251        assert_eq!(f_trigammaf(-27.058018), 300.35904);
252        assert_eq!(f_trigammaf(27.058018), 0.037648965);
253        assert_eq!(f_trigammaf(8.058018), 0.13211797);
254        assert_eq!(f_trigammaf(-8.058018), 300.27863);
255        assert_eq!(f_trigammaf(2.23432), 0.56213206);
256        assert_eq!(f_trigammaf(-2.4653), 9.653673);
257        assert_eq!(f_trigammaf(0.123541), 66.911285);
258        assert_eq!(f_trigammaf(-0.54331), 9.154416);
259        assert_eq!(f_trigammaf(-5.), f32::INFINITY);
260        assert_eq!(f_trigammaf(-0.), f32::INFINITY);
261        assert_eq!(f_trigammaf(0.), f32::INFINITY);
262        assert_eq!(f_trigammaf(f32::INFINITY), 0.0);
263        assert_eq!(f_trigammaf(f32::NEG_INFINITY), f32::NEG_INFINITY);
264        assert!(f_trigammaf(f32::NAN).is_nan());
265    }
266}
pxfm/gamma/trigammaf.rs

pxfm/gamma/
trigammaf.rs
