pxfm/gamma/digammaf.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::logs::simple_fast_log;
31use crate::polyeval::{
32 f_estrin_polyeval7, f_estrin_polyeval8, f_estrin_polyeval9, f_polyeval4, f_polyeval6,
33 f_polyeval10,
34};
35use crate::tangent::cotpif_core;
36
37#[inline]
38fn approx_digamma(x: f64) -> f64 {
39 if x <= 1. {
40 // Generated in Wolfram Mathematica:
41 // <<FunctionApproximations`
42 // ClearAll["Global`*"]
43 // f[x_]:=PolyGamma[x + 1]
44 // {err0,approx, err1}=MiniMaxApproximation[f[z],{z,{-0.99999999,0},9,8},WorkingPrecision->75,MaxIterations->100]
45 // num=Numerator[approx];
46 // den=Denominator[approx];
47 // poly=num;
48 // coeffs=CoefficientList[poly,z];
49 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
50 // poly=den;
51 // coeffs=CoefficientList[poly,z];
52 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
53 let p_num = f_polyeval10(
54 x,
55 f64::from_bits(0xbfe2788cfc6fb619),
56 f64::from_bits(0x3fca347925788707),
57 f64::from_bits(0x3ff887e0f068df69),
58 f64::from_bits(0x3ff543446198b4d2),
59 f64::from_bits(0x3fe03e4455fbad95),
60 f64::from_bits(0x3fb994be8389e4f6),
61 f64::from_bits(0x3f84eb98b830c9b1),
62 f64::from_bits(0x3f4025193ac4ad97),
63 f64::from_bits(0x3ee18c1d683d866a),
64 f64::from_bits(0x3e457cb5b4a07c95),
65 );
66 let p_den = f_estrin_polyeval9(
67 x,
68 f64::from_bits(0x3ff0000000000000),
69 f64::from_bits(0x4003f5f42d95aca8),
70 f64::from_bits(0x4002f96e541d0513),
71 f64::from_bits(0x3ff22c34843313fa),
72 f64::from_bits(0x3fd33574180109bf),
73 f64::from_bits(0x3fa6c07b99ebb277),
74 f64::from_bits(0x3f6cdd7b8fa68926),
75 f64::from_bits(0x3f212b74d39e287f),
76 f64::from_bits(0x3ebabd247f366583),
77 );
78 return p_num / p_den - 1. / x;
79 } else if x < 1.461632144 {
80 // exception
81 if x == 1.4616321325302124 {
82 return -1.2036052e-8;
83 }
84 // Generated in Wolfram Mathematica:
85 // <<FunctionApproximations`
86 // ClearAll["Global`*"]
87 // f[x_]:=PolyGamma[x+1]
88 // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{0,0.461632144},8,8},WorkingPrecision->75,MaxIterations->100]
89 // num=Numerator[approx];
90 // den=Denominator[approx];
91 // poly=num;
92 // coeffs=CoefficientList[poly,z];
93 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
94 // poly=den;
95 // coeffs=CoefficientList[poly,z];
96 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
97 let p_num = f_estrin_polyeval9(
98 x,
99 f64::from_bits(0xbfe2788cfc6fb619),
100 f64::from_bits(0x3fd0ad221e8c3b8b),
101 f64::from_bits(0x3ff813be4dee2e90),
102 f64::from_bits(0x3ff2f64cbfa7d1a4),
103 f64::from_bits(0x3fd9a8c4798f426c),
104 f64::from_bits(0x3fb111a34898f6bf),
105 f64::from_bits(0x3f75dd3efac1e579),
106 f64::from_bits(0x3f272596b2582f0d),
107 f64::from_bits(0x3eb9b074f4ca6263),
108 );
109 let p_den = f_estrin_polyeval9(
110 x,
111 f64::from_bits(0x3ff0000000000000),
112 f64::from_bits(0x40032fd3a1fe3a25),
113 f64::from_bits(0x40012969bcd7fef3),
114 f64::from_bits(0x3fee1a267ee7a97a),
115 f64::from_bits(0x3fcc1522178a69a6),
116 f64::from_bits(0x3f9bd89421334af0),
117 f64::from_bits(0x3f5b40bc3203df4c),
118 f64::from_bits(0x3f05ac6be0b79fac),
119 f64::from_bits(0x3e9047c3d8071f18),
120 );
121 return p_num / p_den - 1. / x;
122 } else if x <= 2. {
123 // Generated in Wolfram Mathematica:
124 // <<FunctionApproximations`
125 // ClearAll["Global`*"]
126 // f[x_]:=PolyGamma[x+1]
127 // {err0,approx,err1}=MiniMaxApproximation[f[z],{z,{0.461632144,1},7,6},WorkingPrecision->75,MaxIterations->100]
128 // num=Numerator[approx];
129 // den=Denominator[approx];
130 // poly=num;
131 // coeffs=CoefficientList[poly,z];
132 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
133 // poly=den;
134 // coeffs=CoefficientList[poly,z];
135 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
136 let p_num = f_estrin_polyeval8(
137 x,
138 f64::from_bits(0xbfe2788cfc6fb613),
139 f64::from_bits(0x3fd690553caa1c6b),
140 f64::from_bits(0x3ff721cf4d9e008f),
141 f64::from_bits(0x3fee9f4096f34b09),
142 f64::from_bits(0x3fd055e88830fc71),
143 f64::from_bits(0x3f9e66bceee16960),
144 f64::from_bits(0x3f55d436778b3403),
145 f64::from_bits(0x3eeef6bc214819b3),
146 );
147 let p_den = f_estrin_polyeval8(
148 x,
149 f64::from_bits(0x3ff0000000000000),
150 f64::from_bits(0x4001e96eaab05729),
151 f64::from_bits(0x3ffcb1aa289077da),
152 f64::from_bits(0x3fe5499e89b757b6),
153 f64::from_bits(0x3fbee531a912bca9),
154 f64::from_bits(0x3f84d46f2121ceb7),
155 f64::from_bits(0x3f35abd7eb7440e6),
156 f64::from_bits(0x3ec43bf7c110aad1),
157 );
158 return p_num / p_den - 1. / x;
159 } else if x <= 3. {
160 // Generated in Wolfram Mathematica:
161 // <<FunctionApproximations`
162 // ClearAll["Global`*"]
163 // f[x_]:=PolyGamma[x+1]
164 // {err0,approx}=MiniMaxApproximation[f[z],{z,{1,2},7,6},WorkingPrecision->75,MaxIterations->100]
165 // num=Numerator[approx][[1]];
166 // den=Denominator[approx][[1]];
167 // poly=num;
168 // coeffs=CoefficientList[poly,z];
169 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
170 // poly=den;
171 // coeffs=CoefficientList[poly,z];
172 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
173 let p_num = f_estrin_polyeval8(
174 x,
175 f64::from_bits(0xbfe2788cfc63695f),
176 f64::from_bits(0x3fdb63791eb688ea),
177 f64::from_bits(0x3ff625ed84968583),
178 f64::from_bits(0x3fe900ea36e59e02),
179 f64::from_bits(0x3fc5319409f4fec6),
180 f64::from_bits(0x3f8b6e7cacff2a59),
181 f64::from_bits(0x3f34a7e591bf2af3),
182 f64::from_bits(0x3e9c323866d138db),
183 );
184 let p_den = f_estrin_polyeval7(
185 x,
186 f64::from_bits(0x3ff0000000000000),
187 f64::from_bits(0x4000ddf448e34181),
188 f64::from_bits(0x3ff87188f2414f79),
189 f64::from_bits(0x3fdeff74d18f811a),
190 f64::from_bits(0x3fb1a0cddeb3a320),
191 f64::from_bits(0x3f701050c1344800),
192 f64::from_bits(0x3f10480a4ea8cf57),
193 );
194 return p_num / p_den - 1. / x;
195 } else if x <= 8. {
196 // Generated in Wolfram Mathematica:
197 // <<FunctionApproximations`
198 // ClearAll["Global`*"]
199 // f[x_]:=PolyGamma[x + 1]
200 // {err0,approx}=MiniMaxApproximation[f[z],{z,{2,7},7,7},WorkingPrecision->75,MaxIterations->100]
201 // num=Numerator[approx][[1]];
202 // den=Denominator[approx][[1]];
203 // poly=num;
204 // coeffs=CoefficientList[poly,z];
205 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
206 // poly=den;
207 // coeffs=CoefficientList[poly,z];
208 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
209 let p_num = f_estrin_polyeval8(
210 x,
211 f64::from_bits(0xbfe2788c3725fd5e),
212 f64::from_bits(0x3fde39f54e5a651a),
213 f64::from_bits(0x3ff56983b839b94f),
214 f64::from_bits(0x3fe60d6118d8fc08),
215 f64::from_bits(0x3fc0e889ace69a30),
216 f64::from_bits(0x3f844e10e399bd93),
217 f64::from_bits(0x3f3099741afda7cb),
218 f64::from_bits(0x3eb74a15144af8e9),
219 );
220 let p_den = f_estrin_polyeval8(
221 x,
222 f64::from_bits(0x3ff0000000000000),
223 f64::from_bits(0x4000409ed08c0553),
224 f64::from_bits(0x3ff63746cb6183e3),
225 f64::from_bits(0x3fda1196b1a75351),
226 f64::from_bits(0x3fab4ba9fad2d187),
227 f64::from_bits(0x3f67de6e6987e3a3),
228 f64::from_bits(0x3f0c9d85ca31182e),
229 f64::from_bits(0x3e8b269f154c8f12),
230 );
231 return p_num / p_den - 1. / x;
232 } else if x <= 15. {
233 // Generated in Wolfram Mathematica:
234 // <<FunctionApproximations`
235 // ClearAll["Global`*"]
236 // f[x_]:=PolyGamma[x + 1]
237 // {err0,approx}=MiniMaxApproximation[f[z],{z,{7,14},7,7},WorkingPrecision->75,MaxIterations->100]
238 // num=Numerator[approx][[1]];
239 // den=Denominator[approx][[1]];
240 // poly=num;
241 // coeffs=CoefficientList[poly,z];
242 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
243 // poly=den;
244 // coeffs=CoefficientList[poly,z];
245 // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50}, ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
246 let p_num = f_estrin_polyeval8(
247 x,
248 f64::from_bits(0xbfe272e75f131710),
249 f64::from_bits(0x3fe53fce507de081),
250 f64::from_bits(0x3ff182957d4f961a),
251 f64::from_bits(0x3fd6d1652dea00e9),
252 f64::from_bits(0x3fa45e16488abe0f),
253 f64::from_bits(0x3f5a52a8f3f3663f),
254 f64::from_bits(0x3ef5b767554d208e),
255 f64::from_bits(0x3e6d2393b100353d),
256 );
257 let p_den = f_estrin_polyeval8(
258 x,
259 f64::from_bits(0x3ff0000000000000),
260 f64::from_bits(0x3ffb1f295c6e5fc5),
261 f64::from_bits(0x3feb88eb913eb117),
262 f64::from_bits(0x3fc570f3aed83ff7),
263 f64::from_bits(0x3f8afe819fdfa5a5),
264 f64::from_bits(0x3f3a2cec9041f361),
265 f64::from_bits(0x3ed0549335964bb9),
266 f64::from_bits(0x3e3ebdcb0002d63e),
267 );
268 return p_num / p_den - 1. / x;
269 }
270 // digamma asymptotic expansion
271 // digamma(x) ~ ln(z)+1/(2z)-sum_(n=1)^(infty)(Bernoulli(2n))/(2nz^(2n))
272 // Generated in SageMath:
273 // var('x')
274 //
275 // def bernoulli_terms(x, N):
276 // S = 0
277 // S += QQ(1)/QQ(2)/x
278 // for k in range(1, N+1):
279 // B = bernoulli(2*k)
280 // term = (B / QQ(2*k))*x**(-2*k)
281 // S += term
282 // return S
283 //
284 // terms = bernoulli_terms(x, 5)
285 //
286 // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
287 // for k in range(1, 13):
288 // c = terms.coefficient(x, -k) # coefficient of x^(-k)
289 // if c == 0:
290 // continue
291 // print("f64::from_bits(" + double_to_hex(c) + "),")
292 let rcp = 1. / x;
293 let rcp_sqr = rcp * rcp;
294 let p = f_polyeval6(
295 rcp_sqr,
296 f64::from_bits(0x3fb5555555555555),
297 f64::from_bits(0xbf81111111111111),
298 f64::from_bits(0x3f70410410410410),
299 f64::from_bits(0xbf71111111111111),
300 f64::from_bits(0x3f7f07c1f07c1f08),
301 f64::from_bits(0xbf95995995995996),
302 );
303 let v_log = simple_fast_log(x);
304 v_log - f_fmla(rcp, f64::from_bits(0x3fe0000000000000), p * rcp_sqr)
305}
306
307/// Computes digamma(x)
308pub fn f_digammaf(x: f32) -> f32 {
309 let xb = x.to_bits();
310 // filter out exceptional cases
311 if xb >= 0xffu32 << 23 || xb == 0 {
312 if x.is_infinite() {
313 return if x.is_sign_negative() {
314 f32::NAN
315 } else {
316 f32::INFINITY
317 };
318 }
319 if x.is_nan() {
320 return f32::NAN;
321 }
322 if xb == 0 {
323 return f32::INFINITY;
324 }
325 }
326
327 let ax = x.to_bits() & 0x7fff_ffff;
328
329 if ax <= 0x32abcc77u32 {
330 // |x| < 2e-8
331 // digamma near where x -> 1 ~ Digamma[x] = -euler + O(x)
332 // considering identity Digamma[x+1] = Digamma[x] + 1/x,
333 // hence x < ulp(1) then x+1 ~= 1 that gives
334 // Digamma[x] = Digamma[x+1] - 1/x = -euler - 1/x
335 const EULER: f64 = f64::from_bits(0x3fe2788cfc6fb619);
336 return (-EULER - 1. / x as f64) as f32;
337 } else if ax <= 0x377ba882u32 {
338 // |x| <= 0.000015
339 // Laurent series of digamma(x)
340 // Generated by SageMath:
341 // from mpmath import mp
342 // mp.prec = 150
343 // R = RealField(150)
344 // var('x')
345 // def laurent_terms(x, N):
346 // S = 0
347 // S += -1/x - R(mp.euler)
348 // S1 = 0
349 // for k in range(1, N+1):
350 // zet = R(mp.zeta(k + 1))
351 // term = zet*(-x)**k
352 // S1 += term
353 // return S - S1
354 //
355 // terms = laurent_terms(x, 4)
356 //
357 // coeffs = [RealField(150)(terms.coefficient(x, n)) for n in range(0, terms.degree(x)+1, 1)]
358 // for k in range(1, 13):
359 // c = terms.coefficient(x, k) # coefficient of x^(-k)
360 // if c == 0:
361 // continue
362 // print("f64::from_bits(" + double_to_hex(c) + "),")
363 const EULER: f64 = f64::from_bits(0x3fe2788cfc6fb619);
364 let dx = x as f64;
365 let start = -1. / dx;
366 let neg_dx = -dx;
367 let z = f_polyeval4(
368 neg_dx,
369 f64::from_bits(0x3ffa51a6625307d3),
370 f64::from_bits(0xbff33ba004f00621),
371 f64::from_bits(0x3ff151322ac7d848),
372 f64::from_bits(0xbff097418eca7cce),
373 );
374 let r = f_fmla(z, dx, -EULER) + start;
375 return r as f32;
376 }
377
378 let mut dx = x as f64;
379 let mut result: f64 = 0.;
380 if x < 0. {
381 // at negative integers it's inf
382 if is_integerf(x) {
383 return f32::NAN;
384 }
385
386 // reflection Gamma(1-x) + Gamma(x) = Pi/tan(PI*x)
387 const PI: f64 = f64::from_bits(0x400921fb54442d18);
388 let cot_x_angle = -dx;
389 dx = 1. - dx;
390 result = PI * cotpif_core(cot_x_angle);
391 }
392 let approx = approx_digamma(dx);
393 result += approx;
394 result as f32
395}
396
397#[cfg(test)]
398mod tests {
399 use super::*;
400
401 #[test]
402 fn test_digamma() {
403 assert_eq!(f_digammaf(-13.999000000012591), -996.9182);
404 assert_eq!(f_digammaf(15.3796425), 2.700182);
405 assert_eq!(f_digammaf(0.0005187988), -1928.1058);
406 assert_eq!(f_digammaf(0.0019531252), -512.574);
407 assert_eq!(f_digammaf(-96.353516), 6.1304626);
408 assert_eq!(f_digammaf(-31.06964), 17.582127);
409 assert_eq!(f_digammaf(-0.000000000000001191123), 839543830000000.);
410 assert_eq!(f_digammaf(f32::INFINITY), f32::INFINITY);
411 assert_eq!(f_digammaf(0.), f32::INFINITY);
412 assert_eq!(f_digammaf(-0.), f32::INFINITY);
413 assert!(f_digammaf(f32::NEG_INFINITY).is_nan());
414 assert!(f_digammaf(f32::NAN).is_nan());
415 }
416}