pxfm/gamma/
digamma.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::gamma::digamma_coeffs::digamma_coeefs;
32use crate::logs::fast_log_d_to_dd;
33use crate::tangent::cotpi_core;
34
35#[inline]
36fn approx_digamma_hard(x: f64) -> DoubleDouble {
37    if x <= 12. {
38        let x2 = DoubleDouble::from_exact_mult(x, x);
39        let x4 = DoubleDouble::quick_mult(x2, x2);
40        let x8 = DoubleDouble::quick_mult(x4, x4);
41        let (p, q) = digamma_coeefs(x);
42
43        let e0 = DoubleDouble::mul_f64_add(
44            DoubleDouble::from_bit_pair(p[1]),
45            x,
46            DoubleDouble::from_bit_pair(p[0]),
47        );
48        let e1 = DoubleDouble::mul_f64_add(
49            DoubleDouble::from_bit_pair(p[3]),
50            x,
51            DoubleDouble::from_bit_pair(p[2]),
52        );
53        let e2 = DoubleDouble::mul_f64_add(
54            DoubleDouble::from_bit_pair(p[5]),
55            x,
56            DoubleDouble::from_bit_pair(p[4]),
57        );
58        let e3 = DoubleDouble::mul_f64_add(
59            DoubleDouble::from_bit_pair(p[7]),
60            x,
61            DoubleDouble::from_bit_pair(p[6]),
62        );
63        let e4 = DoubleDouble::mul_f64_add(
64            DoubleDouble::from_bit_pair(p[9]),
65            x,
66            DoubleDouble::from_bit_pair(p[8]),
67        );
68        let e5 = DoubleDouble::mul_f64_add(
69            DoubleDouble::from_bit_pair(p[11]),
70            x,
71            DoubleDouble::from_bit_pair(p[10]),
72        );
73
74        let f0 = DoubleDouble::mul_add(x2, e1, e0);
75        let f1 = DoubleDouble::mul_add(x2, e3, e2);
76        let f2 = DoubleDouble::mul_add(x2, e5, e4);
77
78        let g0 = DoubleDouble::mul_add(x4, f1, f0);
79
80        let p_num = DoubleDouble::mul_add(x8, f2, g0);
81
82        let rcp = DoubleDouble::from_quick_recip(x);
83
84        let e0 = DoubleDouble::mul_f64_add_f64(
85            DoubleDouble::from_bit_pair(q[1]),
86            x,
87            f64::from_bits(0x3ff0000000000000),
88        );
89        let e1 = DoubleDouble::mul_f64_add(
90            DoubleDouble::from_bit_pair(q[3]),
91            x,
92            DoubleDouble::from_bit_pair(q[2]),
93        );
94        let e2 = DoubleDouble::mul_f64_add(
95            DoubleDouble::from_bit_pair(q[5]),
96            x,
97            DoubleDouble::from_bit_pair(q[4]),
98        );
99        let e3 = DoubleDouble::mul_f64_add(
100            DoubleDouble::from_bit_pair(q[7]),
101            x,
102            DoubleDouble::from_bit_pair(q[6]),
103        );
104        let e4 = DoubleDouble::mul_f64_add(
105            DoubleDouble::from_bit_pair(q[9]),
106            x,
107            DoubleDouble::from_bit_pair(q[8]),
108        );
109        let e5 = DoubleDouble::mul_f64_add(
110            DoubleDouble::from_bit_pair(q[11]),
111            x,
112            DoubleDouble::from_bit_pair(q[10]),
113        );
114
115        let f0 = DoubleDouble::mul_add(x2, e1, e0);
116        let f1 = DoubleDouble::mul_add(x2, e3, e2);
117        let f2 = DoubleDouble::mul_add(x2, e5, e4);
118
119        let g0 = DoubleDouble::mul_add(x4, f1, f0);
120
121        let p_den = DoubleDouble::mul_add(x8, f2, g0);
122
123        let q = DoubleDouble::div(p_num, p_den);
124        let r = DoubleDouble::quick_dd_sub(q, rcp);
125        return r;
126    }
127    // digamma asymptotic expansion
128    // digamma(x) ~ ln(z)+1/(2z)-sum_(n=1)^(infty)(Bernoulli(2n))/(2nz^(2n))
129    // Generated in SageMath:
130    // var('x')
131    //
132    // def bernoulli_terms(x, N):
133    //     S = 0
134    //     S += QQ(1)/QQ(2)/x
135    //     for k in range(1, N+1):
136    //         B = bernoulli(2*k)
137    //         term = (B / QQ(2*k))*x**(-2*k)
138    //         S += term
139    //     return S
140    //
141    // terms = bernoulli_terms(x, 9)
142    //
143    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
144    // for k in range(1, 13):
145    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
146    //     if c == 0:
147    //         continue
148    //     print("f64::from_bits(" + double_to_hex(c) + "),")
149    const C: [(u64, u64); 10] = [
150        (0x3c55555555555555, 0x3fb5555555555555),
151        (0xbc01111111111111, 0xbf81111111111111),
152        (0x3c10410410410410, 0x3f70410410410410),
153        (0xbbf1111111111111, 0xbf71111111111111),
154        (0xbc0f07c1f07c1f08, 0x3f7f07c1f07c1f08),
155        (0x3c39a99a99a99a9a, 0xbf95995995995996),
156        (0x3c55555555555555, 0x3fb5555555555555),
157        (0xbc77979797979798, 0xbfdc5e5e5e5e5e5e),
158        (0xbc69180646019180, 0x40086e7f9b9fe6e8),
159        (0x3ccad759ad759ad7, 0xc03a74ca514ca515),
160    ];
161
162    let rcp = DoubleDouble::from_quick_recip(x);
163    let rcp_sqr = DoubleDouble::quick_mult(rcp, rcp);
164
165    let rx = rcp_sqr;
166    let x2 = DoubleDouble::quick_mult(rx, rx);
167    let x4 = DoubleDouble::quick_mult(x2, x2);
168    let x8 = DoubleDouble::quick_mult(x4, x4);
169
170    let q0 = DoubleDouble::quick_mul_add(
171        DoubleDouble::from_bit_pair(C[1]),
172        rx,
173        DoubleDouble::from_bit_pair(C[0]),
174    );
175    let q1 = DoubleDouble::quick_mul_add(
176        DoubleDouble::from_bit_pair(C[3]),
177        rx,
178        DoubleDouble::from_bit_pair(C[2]),
179    );
180    let q2 = DoubleDouble::quick_mul_add(
181        DoubleDouble::from_bit_pair(C[5]),
182        rx,
183        DoubleDouble::from_bit_pair(C[4]),
184    );
185    let q3 = DoubleDouble::quick_mul_add(
186        DoubleDouble::from_bit_pair(C[7]),
187        rx,
188        DoubleDouble::from_bit_pair(C[6]),
189    );
190    let q4 = DoubleDouble::quick_mul_add(
191        DoubleDouble::from_bit_pair(C[9]),
192        rx,
193        DoubleDouble::from_bit_pair(C[8]),
194    );
195
196    let q0 = DoubleDouble::quick_mul_add(x2, q1, q0);
197    let q1 = DoubleDouble::quick_mul_add(x2, q3, q2);
198
199    let r0 = DoubleDouble::quick_mul_add(x4, q1, q0);
200    let mut p = DoubleDouble::quick_mul_add(x8, q4, r0);
201    p = DoubleDouble::quick_mul_f64_add(
202        rcp,
203        f64::from_bits(0x3fe0000000000000),
204        DoubleDouble::quick_mult(p, rcp_sqr),
205    );
206
207    let v_log = fast_log_d_to_dd(x);
208    DoubleDouble::quick_dd_sub(v_log, p)
209}
210
211#[inline]
212fn approx_digamma_hard_dd(x: DoubleDouble) -> DoubleDouble {
213    if x.hi <= 12. {
214        let x2 = DoubleDouble::quick_mult(x, x);
215        let x4 = DoubleDouble::quick_mult(x2, x2);
216        let x8 = DoubleDouble::quick_mult(x4, x4);
217        let (p, q) = digamma_coeefs(x.hi);
218
219        let e0 = DoubleDouble::mul_add(
220            DoubleDouble::from_bit_pair(p[1]),
221            x,
222            DoubleDouble::from_bit_pair(p[0]),
223        );
224        let e1 = DoubleDouble::mul_add(
225            DoubleDouble::from_bit_pair(p[3]),
226            x,
227            DoubleDouble::from_bit_pair(p[2]),
228        );
229        let e2 = DoubleDouble::mul_add(
230            DoubleDouble::from_bit_pair(p[5]),
231            x,
232            DoubleDouble::from_bit_pair(p[4]),
233        );
234        let e3 = DoubleDouble::mul_add(
235            DoubleDouble::from_bit_pair(p[7]),
236            x,
237            DoubleDouble::from_bit_pair(p[6]),
238        );
239        let e4 = DoubleDouble::mul_add(
240            DoubleDouble::from_bit_pair(p[9]),
241            x,
242            DoubleDouble::from_bit_pair(p[8]),
243        );
244        let e5 = DoubleDouble::mul_add(
245            DoubleDouble::from_bit_pair(p[11]),
246            x,
247            DoubleDouble::from_bit_pair(p[10]),
248        );
249
250        let f0 = DoubleDouble::mul_add(x2, e1, e0);
251        let f1 = DoubleDouble::mul_add(x2, e3, e2);
252        let f2 = DoubleDouble::mul_add(x2, e5, e4);
253
254        let g0 = DoubleDouble::mul_add(x4, f1, f0);
255
256        let p_num = DoubleDouble::mul_add(x8, f2, g0);
257
258        let rcp = x.recip();
259
260        let e0 = DoubleDouble::mul_add_f64(
261            DoubleDouble::from_bit_pair(q[1]),
262            x,
263            f64::from_bits(0x3ff0000000000000),
264        );
265        let e1 = DoubleDouble::mul_add(
266            DoubleDouble::from_bit_pair(q[3]),
267            x,
268            DoubleDouble::from_bit_pair(q[2]),
269        );
270        let e2 = DoubleDouble::mul_add(
271            DoubleDouble::from_bit_pair(q[5]),
272            x,
273            DoubleDouble::from_bit_pair(q[4]),
274        );
275        let e3 = DoubleDouble::mul_add(
276            DoubleDouble::from_bit_pair(q[7]),
277            x,
278            DoubleDouble::from_bit_pair(q[6]),
279        );
280        let e4 = DoubleDouble::mul_add(
281            DoubleDouble::from_bit_pair(q[9]),
282            x,
283            DoubleDouble::from_bit_pair(q[8]),
284        );
285        let e5 = DoubleDouble::mul_add(
286            DoubleDouble::from_bit_pair(q[11]),
287            x,
288            DoubleDouble::from_bit_pair(q[10]),
289        );
290
291        let f0 = DoubleDouble::mul_add(x2, e1, e0);
292        let f1 = DoubleDouble::mul_add(x2, e3, e2);
293        let f2 = DoubleDouble::mul_add(x2, e5, e4);
294
295        let g0 = DoubleDouble::mul_add(x4, f1, f0);
296
297        let p_den = DoubleDouble::mul_add(x8, f2, g0);
298
299        let q = DoubleDouble::div(p_num, p_den);
300        let r = DoubleDouble::quick_dd_sub(q, rcp);
301        return r;
302    }
303    // digamma asymptotic expansion
304    // digamma(x) ~ ln(z)+1/(2z)-sum_(n=1)^(infty)(Bernoulli(2n))/(2nz^(2n))
305    // Generated in SageMath:
306    // var('x')
307    //
308    // def bernoulli_terms(x, N):
309    //     S = 0
310    //     S += QQ(1)/QQ(2)/x
311    //     for k in range(1, N+1):
312    //         B = bernoulli(2*k)
313    //         term = (B / QQ(2*k))*x**(-2*k)
314    //         S += term
315    //     return S
316    //
317    // terms = bernoulli_terms(x, 9)
318    //
319    // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
320    // for k in range(1, 13):
321    //     c = terms.coefficient(x, -k)  # coefficient of x^(-k)
322    //     if c == 0:
323    //         continue
324    //     print("f64::from_bits(" + double_to_hex(c) + "),")
325    const C: [(u64, u64); 10] = [
326        (0x3c55555555555555, 0x3fb5555555555555),
327        (0xbc01111111111111, 0xbf81111111111111),
328        (0x3c10410410410410, 0x3f70410410410410),
329        (0xbbf1111111111111, 0xbf71111111111111),
330        (0xbc0f07c1f07c1f08, 0x3f7f07c1f07c1f08),
331        (0x3c39a99a99a99a9a, 0xbf95995995995996),
332        (0x3c55555555555555, 0x3fb5555555555555),
333        (0xbc77979797979798, 0xbfdc5e5e5e5e5e5e),
334        (0xbc69180646019180, 0x40086e7f9b9fe6e8),
335        (0x3ccad759ad759ad7, 0xc03a74ca514ca515),
336    ];
337
338    let rcp = x.recip();
339    let rcp_sqr = DoubleDouble::quick_mult(rcp, rcp);
340
341    let rx = rcp_sqr;
342    let x2 = DoubleDouble::quick_mult(rx, rx);
343    let x4 = DoubleDouble::quick_mult(x2, x2);
344    let x8 = DoubleDouble::quick_mult(x4, x4);
345
346    let q0 = DoubleDouble::quick_mul_add(
347        DoubleDouble::from_bit_pair(C[1]),
348        rx,
349        DoubleDouble::from_bit_pair(C[0]),
350    );
351    let q1 = DoubleDouble::quick_mul_add(
352        DoubleDouble::from_bit_pair(C[3]),
353        rx,
354        DoubleDouble::from_bit_pair(C[2]),
355    );
356    let q2 = DoubleDouble::quick_mul_add(
357        DoubleDouble::from_bit_pair(C[5]),
358        rx,
359        DoubleDouble::from_bit_pair(C[4]),
360    );
361    let q3 = DoubleDouble::quick_mul_add(
362        DoubleDouble::from_bit_pair(C[7]),
363        rx,
364        DoubleDouble::from_bit_pair(C[6]),
365    );
366    let q4 = DoubleDouble::quick_mul_add(
367        DoubleDouble::from_bit_pair(C[9]),
368        rx,
369        DoubleDouble::from_bit_pair(C[8]),
370    );
371
372    let q0 = DoubleDouble::quick_mul_add(x2, q1, q0);
373    let q1 = DoubleDouble::quick_mul_add(x2, q3, q2);
374
375    let r0 = DoubleDouble::quick_mul_add(x4, q1, q0);
376    let mut p = DoubleDouble::quick_mul_add(x8, q4, r0);
377    p = DoubleDouble::quick_mul_f64_add(
378        rcp,
379        f64::from_bits(0x3fe0000000000000),
380        DoubleDouble::quick_mult(p, rcp_sqr),
381    );
382
383    let mut v_log = fast_log_d_to_dd(x.hi);
384    v_log.lo += x.lo / x.hi;
385    DoubleDouble::quick_dd_sub(v_log, p)
386}
387
388/// Computes digamma(x)
389///
390pub fn f_digamma(x: f64) -> f64 {
391    let xb = x.to_bits();
392    if !x.is_normal() {
393        if x.is_infinite() {
394            return if x.is_sign_negative() {
395                f64::NAN
396            } else {
397                f64::INFINITY
398            };
399        }
400        if x.is_nan() {
401            return f64::NAN;
402        }
403        if xb.wrapping_shl(1) == 0 {
404            // |x| == 0
405            return f64::INFINITY;
406        }
407    }
408
409    let dx = x;
410
411    if x.abs() <= f64::EPSILON {
412        // |x| < ulp(1)
413        // digamma near where x -> 1 ~ Digamma[x] = -euler + O(x)
414        // considering identity Digamma[x+1] = Digamma[x] + 1/x,
415        // hence x < ulp(1) then x+1 ~= 1 that gives
416        // Digamma[x] = Digamma[x+1] - 1/x = -euler - 1/x
417        // euler is dropped since 1/x >> euler
418        // that gives:
419        // Digamma[x] = Digamma[x+1] - 1/x = -1/x
420        return -1. / x;
421    }
422
423    if x < 0. {
424        // at negative integers it's inf
425        if is_integer(x) {
426            return f64::NAN;
427        }
428
429        // reflection Gamma(1-x) + Gamma(x) = Pi/tan(PI*x)
430        const PI: DoubleDouble =
431            DoubleDouble::from_bit_pair((0x3ca1a62633145c07, 0x400921fb54442d18));
432        let r = DoubleDouble::from_full_exact_sub(1., x);
433        let mut result = PI * cotpi_core(-x);
434        let app = approx_digamma_hard_dd(r);
435        result = DoubleDouble::quick_dd_add(result, app);
436        result.to_f64()
437    } else {
438        let app = approx_digamma_hard(dx);
439        app.to_f64()
440    }
441}
442
443#[cfg(test)]
444mod tests {
445    use super::*;
446    #[test]
447    fn test_digamma() {
448        assert_eq!(f_digamma(0.0019531200000040207), -512.5753182109892);
449        assert_eq!(f_digamma(-13.999000000012591), -997.3224450000563);
450        assert_eq!(f_digamma(13.999000000453323), 2.602844047257257);
451        assert_eq!(f_digamma(0.), f64::INFINITY);
452        assert_eq!(f_digamma(-0.), f64::INFINITY);
453        assert!(f_digamma(f64::NAN).is_nan());
454    }
455}
pxfm/gamma/digamma.rs

pxfm/gamma/
digamma.rs
