pxfm/gamma/
trigamma.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::is_integer;
30use crate::double_double::DoubleDouble;
31use crate::sincospi::f_fast_sinpi_dd;
32
33// Generated in Wolfram Mathematica:
34// <<FunctionApproximations`
35// ClearAll["Global`*"]
36// f[x_]:=PolyGamma[1, x+1]
37// {err0,approx}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},11,11},WorkingPrecision->75,MaxIterations->100]
38// num=Numerator[approx][[1]];
39// poly=num;
40// coeffs=CoefficientList[poly,z];
41// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
42static P_1: [(u64, u64); 12] = [
43    (0x3c81873d89121fec, 0x3ffa51a6625307d3),
44    (0x3cb78bf7af507504, 0x4016a65cbac476ca),
45    (0xbc94179c65e021c2, 0x40218234a0a79582),
46    (0x3cb842a8ab5e0994, 0x401fde32175f8515),
47    (0xbc8768b33f5776b7, 0x4012de6bde49abff),
48    (0x3c8e06354a27f081, 0x3ffe5d6ef3a2eac6),
49    (0xbc8b391d09fed17a, 0x3fe0d186688252cf),
50    (0xbc59a5c46bb8b8cc, 0x3fb958f9a0f156b7),
51    (0x3c2ac44c6a197244, 0x3f88e605f24e1a89),
52    (0x3bd2d05fa8be27f2, 0x3f4cd369f5d68104),
53    (0x3b84ad0a748fdd22, 0x3efde955ebb17874),
54    (0xb96aa8b9a65e0899, 0xbce053d04459ead7),
55];
56
57// Generated by Wolfram Mathematica:
58// <<FunctionApproximations`
59// ClearAll["Global`*"]
60// f[x_]:=PolyGamma[1, x+1]
61// {err0,approx}=MiniMaxApproximation[f[z],{z,{0,3},11,11},WorkingPrecision->75,MaxIterations->100]
62// num=Numerator[approx][[1]];
63// poly=num;
64// coeffs=CoefficientList[poly,z];
65// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
66static P_2: [(u64, u64); 12] = [
67    (0x3c81873d8912236c, 0x3ffa51a6625307d3),
68    (0xbcbc731a62342288, 0x40176c23f7ea51e6),
69    (0xbcc45a6fd00e67a8, 0x4022cb5ae67657ef),
70    (0x3cc0876fde7fe4e6, 0x4021d766062b9550),
71    (0xbcaec4a4859cba1d, 0x401629f91cd4f291),
72    (0x3c76184014e4d7e3, 0x4002d43da3352004),
73    (0x3c812c7609483e0e, 0x3fe62e3266eef8c7),
74    (0xbc5f991047f52d2b, 0x3fc1eacb910b951c),
75    (0x3c28b9f38d603f2f, 0x3f930960a301df34),
76    (0x3bf9b620eb930504, 0x3f5814f8e057b14b),
77    (0xbb990860b88b54e4, 0x3f0b9f67c71aa3bf),
78    (0x38e5cb6acfbaab77, 0xbc4194b8c01afe9a),
79];
80
81// Generated by Wolfram Mathematica:
82// <<FunctionApproximations`
83// ClearAll["Global`*"]
84// f[x_]:=PolyGamma[1, x+1]
85// {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{3,9},11,11},WorkingPrecision->75,MaxIterations->100]
86// num=Numerator[approx];
87// poly=num;
88// coeffs=CoefficientList[poly,z];
89// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
90static P_3: [(u64, u64); 12] = [
91    (0x3c9cd56dbb295efc, 0x3ffa51a662556db9),
92    (0x3c9f4ee74f5f9daf, 0x4018ff913088cb34),
93    (0x3ccf08737350609c, 0x402593d55686b8b1),
94    (0xbcc6cd4ed33afebb, 0x402641d10de4def5),
95    (0xbcb24d1957c1303c, 0x401e682c37e8e2cf),
96    (0x3ca30ac79162ceb2, 0x400ccfc7c4566f55),
97    (0x3c9efea5ff293dc9, 0x3ff33eb2c6e89d0b),
98    (0x3c74670a11068abc, 0x3fd1fbf456e5c6f0),
99    (0x3c47b5dcdea19c36, 0x3fa6a6a2148c482c),
100    (0xbc14642012a1cc1e, 0x3f71851e927f52e7),
101    (0x3bc7db88a4ec5478, 0x3f29a45059a43475),
102    (0xb7bc31e55271eab0, 0xbb375518529c52fb),
103];
104
105// Generated in Wolfram Mathematica:
106// <<FunctionApproximations`
107// ClearAll["Global`*"]
108// f[x_]:=PolyGamma[1, x+1]
109// {err0,approx}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},11,11},WorkingPrecision->75,MaxIterations->100]
110// den=Denominator[approx][[1]];
111// poly=den;
112// coeffs=CoefficientList[poly,z];
113// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
114static Q_1: [(u64, u64); 12] = [
115    (0x0000000000000000, 0x3ff0000000000000),
116    (0xbcb84c43a11fc28a, 0x40139d9587da0fb5),
117    (0x3ca1cf3dbcddbb57, 0x402507cb6225f0f0),
118    (0x3cb01aa6ddcc3cfd, 0x402a1b416d0ed4e6),
119    (0xbcbc31c216b5ff66, 0x4024ec8829e535d4),
120    (0xbcb335c23022f43e, 0x4016d2ba6d1a18e6),
121    (0x3cafbfffc03ad28a, 0x400158c4611ed51f),
122    (0xbc8d1fb10a031a27, 0x3fe26f15bb52f89b),
123    (0x3c56a9fea160eecb, 0x3fbaec13f663049d),
124    (0xbc2f4ee869ba9364, 0x3f89cde0500cd68f),
125    (0x3be3f23afc9398b6, 0x3f4d4b0f4dcf3eb8),
126    (0x3b67a3ed4795d33e, 0x3efde955eafac9c2),
127];
128
129// Generated by Wolfram Mathematica:
130// <<FunctionApproximations`
131// ClearAll["Global`*"]
132// f[x_]:=PolyGamma[1, x+1]
133// {err0,approx}=MiniMaxApproximation[f[z],{z,{0,3},11,11},WorkingPrecision->75,MaxIterations->100]
134// den=Denominator[approx][[1]];
135// poly=den;
136// coeffs=CoefficientList[poly,z];
137// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
138static Q_2: [(u64, u64); 12] = [
139    (0x0000000000000000, 0x3ff0000000000000),
140    (0xbc81a3e4e026b7b1, 0x401415d1a20a9339),
141    (0x3cc279576dfe3ec9, 0x402627c1a95d33d2),
142    (0x3c9a94b5cf0cae88, 0x402c724dc5cf4577),
143    (0xbc8aa1fa0c3820a8, 0x4027b7a332bb07f4),
144    (0xbc96968367088d66, 0x401b14376177bdd7),
145    (0x3ca2d3dfa5847f4d, 0x4005b0511cd98f2c),
146    (0xbc8cfad394d41dd1, 0x3fe877bc2d02c7f3),
147    (0xbc51592b8ec81a92, 0x3fc31f52afc72b95),
148    (0x3c2cbef277d587e9, 0x3f93cb2f0e574376),
149    (0xbbfbb670fd94f6ba, 0x3f5883767f745a92),
150    (0xbb931b04d74e5893, 0x3f0b9f67c71a60f3),
151];
152
153// Generated by Wolfram Mathematica:
154// <<FunctionApproximations`
155// ClearAll["Global`*"]
156// f[x_]:=PolyGamma[1, x+1]
157// {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{3,9},11,11},WorkingPrecision->75,MaxIterations->100]
158// den=Denominator[approx];
159// poly=den;
160// coeffs=CoefficientList[poly,z];
161// TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
162static Q_3: [(u64, u64); 12] = [
163    (0x0000000000000000, 0x3ff0000000000000),
164    (0x3cbafcb4b6d646d9, 0x40150b12a79fc9cf),
165    (0x3c989ef814b8dd2a, 0x40288c1d26ffdca5),
166    (0x3cb0282cfea9c473, 0x4030d737c893cd5f),
167    (0x3cc955b8aaadb37d, 0x402e5c289b6de3e0),
168    (0x3cb377161f8861d2, 0x4022fb66d87bd522),
169    (0xbcb4b0e4cff46ad6, 0x4010e5d13c2a5907),
170    (0xbc8824539e4b1bd6, 0x3ff58c8fe8f26fca),
171    (0xbc7d34220d810ea0, 0x3fd36c1351f43e66),
172    (0xbc4cbdbe85570017, 0x3fa7c1170466605e),
173    (0xbc0c3afb98775c53, 0x3f71ebafd3e5e3b9),
174    (0x3bc0b0b7f16afd0a, 0x3f29a45059a43475),
175];
176
177#[inline]
178fn approx_trigamma(x: f64) -> DoubleDouble {
179    if x <= 10. {
180        let (p, q) = if x <= 1. {
181            (&P_1, &Q_1)
182        } else if x <= 4. {
183            (&P_2, &Q_2)
184        } else {
185            (&P_3, &Q_3)
186        };
187        let x2 = DoubleDouble::from_exact_mult(x, x);
188        let x4 = DoubleDouble::quick_mult(x2, x2);
189        let x8 = DoubleDouble::quick_mult(x4, x4);
190
191        let e0 = DoubleDouble::mul_f64_add(
192            DoubleDouble::from_bit_pair(p[1]),
193            x,
194            DoubleDouble::from_bit_pair(p[0]),
195        );
196        let e1 = DoubleDouble::mul_f64_add(
197            DoubleDouble::from_bit_pair(p[3]),
198            x,
199            DoubleDouble::from_bit_pair(p[2]),
200        );
201        let e2 = DoubleDouble::mul_f64_add(
202            DoubleDouble::from_bit_pair(p[5]),
203            x,
204            DoubleDouble::from_bit_pair(p[4]),
205        );
206        let e3 = DoubleDouble::mul_f64_add(
207            DoubleDouble::from_bit_pair(p[7]),
208            x,
209            DoubleDouble::from_bit_pair(p[6]),
210        );
211        let e4 = DoubleDouble::mul_f64_add(
212            DoubleDouble::from_bit_pair(p[9]),
213            x,
214            DoubleDouble::from_bit_pair(p[8]),
215        );
216        let e5 = DoubleDouble::mul_f64_add(
217            DoubleDouble::from_bit_pair(p[11]),
218            x,
219            DoubleDouble::from_bit_pair(p[10]),
220        );
221
222        let f0 = DoubleDouble::mul_add(x2, e1, e0);
223        let f1 = DoubleDouble::mul_add(x2, e3, e2);
224        let f2 = DoubleDouble::mul_add(x2, e5, e4);
225
226        let g0 = DoubleDouble::mul_add(x4, f1, f0);
227
228        let p_num = DoubleDouble::mul_add(x8, f2, g0);
229
230        let rcp = DoubleDouble::from_quick_recip(x);
231        let rcp2 = DoubleDouble::quick_mult(rcp, rcp);
232
233        let e0 = DoubleDouble::mul_f64_add_f64(
234            DoubleDouble::from_bit_pair(q[1]),
235            x,
236            f64::from_bits(0x3ff0000000000000),
237        );
238        let e1 = DoubleDouble::mul_f64_add(
239            DoubleDouble::from_bit_pair(q[3]),
240            x,
241            DoubleDouble::from_bit_pair(q[2]),
242        );
243        let e2 = DoubleDouble::mul_f64_add(
244            DoubleDouble::from_bit_pair(q[5]),
245            x,
246            DoubleDouble::from_bit_pair(q[4]),
247        );
248        let e3 = DoubleDouble::mul_f64_add(
249            DoubleDouble::from_bit_pair(q[7]),
250            x,
251            DoubleDouble::from_bit_pair(q[6]),
252        );
253        let e4 = DoubleDouble::mul_f64_add(
254            DoubleDouble::from_bit_pair(q[9]),
255            x,
256            DoubleDouble::from_bit_pair(q[8]),
257        );
258        let e5 = DoubleDouble::mul_f64_add(
259            DoubleDouble::from_bit_pair(q[11]),
260            x,
261            DoubleDouble::from_bit_pair(q[10]),
262        );
263
264        let f0 = DoubleDouble::mul_add(x2, e1, e0);
265        let f1 = DoubleDouble::mul_add(x2, e3, e2);
266        let f2 = DoubleDouble::mul_add(x2, e5, e4);
267
268        let g0 = DoubleDouble::mul_add(x4, f1, f0);
269
270        let p_den = DoubleDouble::mul_add(x8, f2, g0);
271
272        let q = DoubleDouble::div(p_num, p_den);
273        let r = DoubleDouble::quick_dd_add(q, rcp2);
274        return r;
275    }
276    // asymptotic expansion Trigamma[x] = 1/x + 1/x^2 + sum(Bernoulli(2*k)/x^(2*k + 1))
277    // Generated in SageMath:
278    // var('x')
279    // def bernoulli_terms(x, N):
280    //     S = 0
281    //     for k in range(1, N+1):
282    //         B = bernoulli(2*k)
283    //         term = B*x**(-(2*k+1))
284    //         S += term
285    //     return S
286    //
287    // terms = bernoulli_terms(x, 10)
288    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
289    // for k in range(0, 14):
290    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
291    //     if c == 0:
292    //         continue
293    //     print("f64::from_bits(" + double_to_hex(c) + "),")
294    const C: [(u64, u64); 10] = [
295        (0x3c65555555555555, 0x3fc5555555555555),
296        (0xbc21111111111111, 0xbfa1111111111111),
297        (0x3c38618618618618, 0x3f98618618618618),
298        (0xbc21111111111111, 0xbfa1111111111111),
299        (0xbc4364d9364d9365, 0x3fb364d9364d9365),
300        (0xbc6981981981981a, 0xbfd0330330330330),
301        (0xbc95555555555555, 0x3ff2aaaaaaaaaaab),
302        (0xbcb7979797979798, 0xc01c5e5e5e5e5e5e),
303        (0xbcac3b070ec1c3b0, 0x404b7c4f8f13e3c5),
304        (0x3cc8d3018d3018d3, 0xc08088fe72cfe72d),
305    ];
306
307    let rcp = DoubleDouble::from_quick_recip(x);
308
309    let q = DoubleDouble::quick_mult(rcp, rcp);
310
311    let q2 = DoubleDouble::quick_mult(q, q);
312    let q4 = q2 * q2;
313    let q8 = q4 * q4;
314
315    let e0 = DoubleDouble::quick_mul_add(
316        DoubleDouble::from_bit_pair(C[1]),
317        q,
318        DoubleDouble::from_bit_pair(C[0]),
319    );
320    let e1 = DoubleDouble::quick_mul_add(
321        DoubleDouble::from_bit_pair(C[3]),
322        q,
323        DoubleDouble::from_bit_pair(C[2]),
324    );
325    let e2 = DoubleDouble::quick_mul_add(
326        DoubleDouble::from_bit_pair(C[5]),
327        q,
328        DoubleDouble::from_bit_pair(C[4]),
329    );
330    let e3 = DoubleDouble::quick_mul_add(
331        DoubleDouble::from_bit_pair(C[7]),
332        q,
333        DoubleDouble::from_bit_pair(C[6]),
334    );
335    let e4 = DoubleDouble::quick_mul_add(
336        DoubleDouble::from_bit_pair(C[9]),
337        q,
338        DoubleDouble::from_bit_pair(C[8]),
339    );
340
341    let q0 = DoubleDouble::quick_mul_add(q2, e1, e0);
342    let q1 = DoubleDouble::quick_mul_add(q2, e3, e2);
343
344    let r0 = DoubleDouble::quick_mul_add(q4, q1, q0);
345    let mut p = DoubleDouble::quick_mul_add(q8, e4, r0);
346
347    let q_over_2 = DoubleDouble::quick_mult_f64(q, 0.5);
348    p = DoubleDouble::quick_mult(p, q);
349    p = DoubleDouble::quick_mult(p, rcp);
350    p = DoubleDouble::quick_dd_add(q_over_2, p);
351    p = DoubleDouble::quick_dd_add(p, rcp);
352    p
353}
354
355#[inline]
356fn approx_trigamma_dd(x: DoubleDouble) -> DoubleDouble {
357    if x.hi <= 10. {
358        let (p, q) = if x.hi <= 1. {
359            (&P_1, &Q_1)
360        } else if x.hi <= 4. {
361            (&P_2, &Q_2)
362        } else {
363            (&P_3, &Q_3)
364        };
365        let x2 = DoubleDouble::quick_mult(x, x);
366        let x4 = DoubleDouble::quick_mult(x2, x2);
367        let x8 = DoubleDouble::quick_mult(x4, x4);
368
369        let e0 = DoubleDouble::mul_add(
370            DoubleDouble::from_bit_pair(p[1]),
371            x,
372            DoubleDouble::from_bit_pair(p[0]),
373        );
374        let e1 = DoubleDouble::mul_add(
375            DoubleDouble::from_bit_pair(p[3]),
376            x,
377            DoubleDouble::from_bit_pair(p[2]),
378        );
379        let e2 = DoubleDouble::mul_add(
380            DoubleDouble::from_bit_pair(p[5]),
381            x,
382            DoubleDouble::from_bit_pair(p[4]),
383        );
384        let e3 = DoubleDouble::mul_add(
385            DoubleDouble::from_bit_pair(p[7]),
386            x,
387            DoubleDouble::from_bit_pair(p[6]),
388        );
389        let e4 = DoubleDouble::mul_add(
390            DoubleDouble::from_bit_pair(p[9]),
391            x,
392            DoubleDouble::from_bit_pair(p[8]),
393        );
394        let e5 = DoubleDouble::mul_add(
395            DoubleDouble::from_bit_pair(p[11]),
396            x,
397            DoubleDouble::from_bit_pair(p[10]),
398        );
399
400        let f0 = DoubleDouble::mul_add(x2, e1, e0);
401        let f1 = DoubleDouble::mul_add(x2, e3, e2);
402        let f2 = DoubleDouble::mul_add(x2, e5, e4);
403
404        let g0 = DoubleDouble::mul_add(x4, f1, f0);
405
406        let p_num = DoubleDouble::mul_add(x8, f2, g0);
407
408        let rcp = x.recip();
409        let rcp2 = DoubleDouble::quick_mult(rcp, rcp);
410
411        let e0 = DoubleDouble::mul_add_f64(
412            DoubleDouble::from_bit_pair(q[1]),
413            x,
414            f64::from_bits(0x3ff0000000000000),
415        );
416        let e1 = DoubleDouble::mul_add(
417            DoubleDouble::from_bit_pair(q[3]),
418            x,
419            DoubleDouble::from_bit_pair(q[2]),
420        );
421        let e2 = DoubleDouble::mul_add(
422            DoubleDouble::from_bit_pair(q[5]),
423            x,
424            DoubleDouble::from_bit_pair(q[4]),
425        );
426        let e3 = DoubleDouble::mul_add(
427            DoubleDouble::from_bit_pair(q[7]),
428            x,
429            DoubleDouble::from_bit_pair(q[6]),
430        );
431        let e4 = DoubleDouble::mul_add(
432            DoubleDouble::from_bit_pair(q[9]),
433            x,
434            DoubleDouble::from_bit_pair(q[8]),
435        );
436        let e5 = DoubleDouble::mul_add(
437            DoubleDouble::from_bit_pair(q[11]),
438            x,
439            DoubleDouble::from_bit_pair(q[10]),
440        );
441
442        let f0 = DoubleDouble::mul_add(x2, e1, e0);
443        let f1 = DoubleDouble::mul_add(x2, e3, e2);
444        let f2 = DoubleDouble::mul_add(x2, e5, e4);
445
446        let g0 = DoubleDouble::mul_add(x4, f1, f0);
447
448        let p_den = DoubleDouble::mul_add(x8, f2, g0);
449
450        let q = DoubleDouble::div(p_num, p_den);
451        let r = DoubleDouble::quick_dd_add(q, rcp2);
452        return r;
453    }
454    // asymptotic expansion Trigamma[x] = 1/x + 1/x^2 + sum(Bernoulli(2*k)/x^(2*k + 1))
455    // Generated in SageMath:
456    // var('x')
457    // def bernoulli_terms(x, N):
458    //     S = 0
459    //     for k in range(1, N+1):
460    //         B = bernoulli(2*k)
461    //         term = B*x**(-(2*k+1))
462    //         S += term
463    //     return S
464    //
465    // terms = bernoulli_terms(x, 10)
466    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
467    // for k in range(0, 14):
468    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
469    //     if c == 0:
470    //         continue
471    //     print("f64::from_bits(" + double_to_hex(c) + "),")
472    const C: [(u64, u64); 10] = [
473        (0x3c65555555555555, 0x3fc5555555555555),
474        (0xbc21111111111111, 0xbfa1111111111111),
475        (0x3c38618618618618, 0x3f98618618618618),
476        (0xbc21111111111111, 0xbfa1111111111111),
477        (0xbc4364d9364d9365, 0x3fb364d9364d9365),
478        (0xbc6981981981981a, 0xbfd0330330330330),
479        (0xbc95555555555555, 0x3ff2aaaaaaaaaaab),
480        (0xbcb7979797979798, 0xc01c5e5e5e5e5e5e),
481        (0xbcac3b070ec1c3b0, 0x404b7c4f8f13e3c5),
482        (0x3cc8d3018d3018d3, 0xc08088fe72cfe72d),
483    ];
484
485    let rcp = x.recip();
486
487    let q = DoubleDouble::quick_mult(rcp, rcp);
488
489    let q2 = DoubleDouble::quick_mult(q, q);
490    let q4 = q2 * q2;
491    let q8 = q4 * q4;
492
493    let e0 = DoubleDouble::quick_mul_add(
494        DoubleDouble::from_bit_pair(C[1]),
495        q,
496        DoubleDouble::from_bit_pair(C[0]),
497    );
498    let e1 = DoubleDouble::quick_mul_add(
499        DoubleDouble::from_bit_pair(C[3]),
500        q,
501        DoubleDouble::from_bit_pair(C[2]),
502    );
503    let e2 = DoubleDouble::quick_mul_add(
504        DoubleDouble::from_bit_pair(C[5]),
505        q,
506        DoubleDouble::from_bit_pair(C[4]),
507    );
508    let e3 = DoubleDouble::quick_mul_add(
509        DoubleDouble::from_bit_pair(C[7]),
510        q,
511        DoubleDouble::from_bit_pair(C[6]),
512    );
513    let e4 = DoubleDouble::quick_mul_add(
514        DoubleDouble::from_bit_pair(C[9]),
515        q,
516        DoubleDouble::from_bit_pair(C[8]),
517    );
518
519    let q0 = DoubleDouble::quick_mul_add(q2, e1, e0);
520    let q1 = DoubleDouble::quick_mul_add(q2, e3, e2);
521
522    let r0 = DoubleDouble::quick_mul_add(q4, q1, q0);
523    let mut p = DoubleDouble::quick_mul_add(q8, e4, r0);
524
525    let q_over_2 = DoubleDouble::quick_mult_f64(q, 0.5);
526    p = DoubleDouble::quick_mult(p, q);
527    p = DoubleDouble::quick_mult(p, rcp);
528    p = DoubleDouble::quick_dd_add(q_over_2, p);
529    p = DoubleDouble::quick_dd_add(p, rcp);
530    p
531}
532
533/// Computes the trigamma function ψ₁(x).
534///
535/// The trigamma function is the second derivative of the logarithm of the gamma function.
536pub fn f_trigamma(x: f64) -> f64 {
537    let xb = x.to_bits();
538    if !x.is_normal() {
539        if x.is_infinite() {
540            return if x.is_sign_negative() {
541                f64::NEG_INFINITY
542            } else {
543                0.
544            };
545        }
546        if x.is_nan() {
547            return f64::NAN;
548        }
549        if xb == 0 {
550            return f64::INFINITY;
551        }
552    }
553
554    let x_e = (x.to_bits() >> 52) & 0x7ff;
555
556    const E_BIAS: u64 = (1u64 << (11 - 1u64)) - 1u64;
557
558    if x_e < E_BIAS - 52 {
559        // |x| < 2^-52
560        let dx = x;
561        return 1. / (dx * dx);
562    }
563
564    if x < 0. {
565        if is_integer(x) {
566            return f64::INFINITY;
567        }
568        // reflection formula
569        // Trigamma[1-x] + Trigamma[x] = PI^2 / sinpi^2(x)
570        const SQR_PI: DoubleDouble =
571            DoubleDouble::from_bit_pair((0x3cc692b71366cc05, 0x4023bd3cc9be45de)); // pi^2
572        let sinpi_ax = f_fast_sinpi_dd(-x);
573        let dx = DoubleDouble::from_full_exact_sub(1., x);
574        let result = DoubleDouble::div(SQR_PI, DoubleDouble::quick_mult(sinpi_ax, sinpi_ax));
575        let trigamma_x = approx_trigamma_dd(dx);
576        return DoubleDouble::quick_dd_sub(result, trigamma_x).to_f64();
577    }
578
579    approx_trigamma(x).to_f64()
580}
581
582#[cfg(test)]
583mod tests {
584
585    use super::*;
586
587    #[test]
588    fn test_trigamma() {
589        assert_eq!(f_trigamma(-27.058018), 300.35629698636757);
590        assert_eq!(f_trigamma(27.058018), 0.037648965757704725);
591        assert_eq!(f_trigamma(8.058018), 0.13211796975281037);
592        assert_eq!(f_trigamma(-8.058018), 300.2758629255111);
593        assert_eq!(f_trigamma(2.23432), 0.5621320243666134);
594        assert_eq!(f_trigamma(-2.4653), 9.653674003034206);
595        assert_eq!(f_trigamma(0.123541), 66.91128231455282);
596        assert_eq!(f_trigamma(-0.54331), 9.154415950366596);
597        assert_eq!(f_trigamma(-5.), f64::INFINITY);
598        assert_eq!(f_trigamma(f64::INFINITY), 0.0);
599        assert_eq!(f_trigamma(f64::NEG_INFINITY), f64::NEG_INFINITY);
600        assert!(f_trigamma(f64::NAN).is_nan());
601    }
602}
pxfm/gamma/trigamma.rs

pxfm/gamma/
trigamma.rs
