pxfm/gamma/
lgamma.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_integer, is_odd_integer};
30use crate::double_double::DoubleDouble;
31use crate::f_log;
32use crate::floor::FloorFinite;
33use crate::logs::{fast_log_d_to_dd, fast_log_dd, log_dd};
34use crate::polyeval::{f_polyeval4, f_polyeval5, f_polyeval6, f_polyeval10};
35use crate::sincospi::f_fast_sinpi_dd;
36
37#[inline]
38fn apply_sign_and_sum(z: DoubleDouble, parity: f64, s: DoubleDouble) -> DoubleDouble {
39    if parity >= 0. {
40        z
41    } else {
42        DoubleDouble::full_dd_sub(s, z)
43    }
44}
45
46#[inline]
47fn apply_sign_and_sum_quick(z: DoubleDouble, parity: f64, s: DoubleDouble) -> DoubleDouble {
48    if parity >= 0. {
49        z
50    } else {
51        DoubleDouble::quick_dd_sub(s, z)
52    }
53}
54
55#[cold]
56fn lgamma_0p5(dx: f64, v_log: DoubleDouble, f_res: DoubleDouble, sum_parity: f64) -> DoubleDouble {
57    // Log[Gamma[x+1]]
58    // <<FunctionApproximations`
59    // ClearAll["Global`*"]
60    // f[x_]:=LogGamma[x+1]
61    // {err0,approx}=MiniMaxApproximation[f[z],{z,{0.0000000000000001,0.5},8,7},WorkingPrecision->90]
62    // num=Numerator[approx][[1]];
63    // den=Denominator[approx][[1]];
64    // poly=den;
65    // coeffs=CoefficientList[poly,z];
66    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
67    const P: [(u64, u64); 9] = [
68        (0x349d90fba23c9118, 0xb7f552eee31fa8d2),
69        (0x3c56cb95ec26b8b0, 0xbfe2788cfc6fb619),
70        (0xbc98554b6e8ebe5c, 0xbff15a4f25fcc40b),
71        (0x3c51b37126e8f9ee, 0xbfc2d93ef9720645),
72        (0xbc85a532dd358f5d, 0x3fec506397ef590a),
73        (0x3c7dc535e77ac796, 0x3fe674f6812154ca),
74        (0x3c4ff02dbae30f4d, 0x3fc9aacc2b0173a0),
75        (0xbc3aa29e4d4e6c9d, 0x3f95d8cbe81572b6),
76        (0x3bed03960179ad28, 0x3f429e0ecfd47eb2),
77    ];
78
79    let mut p_num = DoubleDouble::mul_f64_add(
80        DoubleDouble::from_bit_pair(P[8]),
81        dx,
82        DoubleDouble::from_bit_pair(P[7]),
83    );
84    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[6]));
85    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[5]));
86    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[4]));
87    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[3]));
88    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[2]));
89    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[1]));
90    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[0]));
91
92    const Q: [(u64, u64); 8] = [
93        (0x0000000000000000, 0x3ff0000000000000),
94        (0xbc93d5d3e154eb7a, 0x400a6e37d475364e),
95        (0xbcb465441d96e6c2, 0x401112f3f3b083a7),
96        (0x3ca1b893e8188325, 0x4005cbfc6085847f),
97        (0xbc8608a5840fb86b, 0x3fec920358d35f3a),
98        (0x3c5dc5b89a3624bd, 0x3fc21011cbbc6923),
99        (0x3bfbe999bea0b965, 0x3f822ae49ffa14ce),
100        (0x3b90cb3bd523bf32, 0x3f20d6565fe86116),
101    ];
102
103    let mut p_den = DoubleDouble::mul_f64_add(
104        DoubleDouble::from_bit_pair(Q[7]),
105        dx,
106        DoubleDouble::from_bit_pair(Q[6]),
107    );
108    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[5]));
109    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[4]));
110    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[3]));
111    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[2]));
112    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[1]));
113    p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
114    let v0 = DoubleDouble::full_dd_sub(DoubleDouble::div(p_num, p_den), v_log);
115    apply_sign_and_sum(v0, sum_parity, f_res)
116}
117
118#[cold]
119fn lgamma_0p5_to_1(
120    dx: f64,
121    v_log: DoubleDouble,
122    f_res: DoubleDouble,
123    sum_parity: f64,
124) -> DoubleDouble {
125    // Log[Gamma[x+1]]
126    // <<FunctionApproximations`
127    // ClearAll["Global`*"]
128    // f[x_]:=LogGamma[x+1]
129    // {err0,approx}=MiniMaxApproximation[f[z],{z,{0.5,0.99999999999999999},9,9},WorkingPrecision->90]
130    // num=Numerator[approx][[1]];
131    // den=Denominator[approx][[1]];
132    // poly=den;
133    // coeffs=CoefficientList[poly,z];
134    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
135    const P: [(u64, u64); 10] = [
136        (0xb967217bfcc647c2, 0xbcc9b2c47481ff4f),
137        (0xbc7eff78623354d7, 0xbfe2788cfc6fb552),
138        (0xbc9886c6f42886e0, 0xbff3c6a7f676cf4c),
139        (0x3c77ed956ff8e661, 0xbfdaf57ee2a64253),
140        (0x3c8723f3a5de4fd5, 0x3feb961a3d8bbe89),
141        (0x3c8848ddf2e2726f, 0x3fedc9f11015d4ca),
142        (0xbc799f3b76da571b, 0x3fd7ac6e82c07787),
143        (0xbc5cb131b054a5f5, 0x3fb103e7e288f4da),
144        (0x3bfe93ab961d39a4, 0x3f74789410ab2cf5),
145        (0x3babcb1e8a573475, 0x3f1d74e78621d316),
146    ];
147
148    let mut p_num = DoubleDouble::mul_f64_add(
149        DoubleDouble::from_bit_pair(P[9]),
150        dx,
151        DoubleDouble::from_bit_pair(P[8]),
152    );
153    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[7]));
154    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[6]));
155    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[5]));
156    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[4]));
157    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[3]));
158    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[2]));
159    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[1]));
160    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[0]));
161
162    const Q: [(u64, u64); 10] = [
163        (0x0000000000000000, 0x3ff0000000000000),
164        (0xbca99871cf68ff41, 0x400c87945e972c56),
165        (0xbc8292764aa01c02, 0x401477d734273be4),
166        (0xbcaf0f1758e16cb3, 0x400e53c84c87f686),
167        (0xbc901825b170e576, 0x3ff8c99df9c66865),
168        (0x3c78af0564323160, 0x3fd629441413e902),
169        (0x3c3293dd176164f3, 0x3fa42e466fd1464e),
170        (0xbbf3fbc18666280b, 0x3f5f9cd4c60c4e7d),
171        (0xbb4036ba7be37458, 0x3efb2d3ed43aab6f),
172        (0xbaf7a3a7d5321f53, 0xbe665e3fadf143a6),
173    ];
174
175    let mut p_den = DoubleDouble::mul_f64_add(
176        DoubleDouble::from_bit_pair(Q[9]),
177        dx,
178        DoubleDouble::from_bit_pair(Q[8]),
179    );
180    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[7]));
181    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[6]));
182    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[5]));
183    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[4]));
184    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[3]));
185    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[2]));
186    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[1]));
187    p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
188    let v0 = DoubleDouble::full_dd_sub(DoubleDouble::div(p_num, p_den), v_log);
189    apply_sign_and_sum(v0, sum_parity, f_res)
190}
191
192#[cold]
193fn lgamma_1_to_4(
194    dx: f64,
195    v_log: DoubleDouble,
196    f_res: DoubleDouble,
197    sum_parity: f64,
198) -> DoubleDouble {
199    // Log[Gamma[x+1]]
200    // <<FunctionApproximations`
201    // ClearAll["Global`*"]
202    // f[x_]:=LogGamma[x+1]
203    // {err0,approx}=MiniMaxApproximation[f[z],{z,{1.0000000000000000001,4},13,13},WorkingPrecision->90]
204    // num=Numerator[approx][[1]];
205    // den=Denominator[approx][[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    const P: [(u64, u64); 14] = [
210        (0x398e8eb32d541d21, 0xbcf55335b06a5df1),
211        (0xbc60fc919ad95ffb, 0xbfe2788cfc6fb2c4),
212        (0x3c98d8e5e64ea307, 0xbffb5e5c82450af5),
213        (0x3c778a3bdabcdb1f, 0xbff824ecfcecbcd5),
214        (0x3c77768d59db11a0, 0x3fd6a097c5fa0036),
215        (0x3c9f4a74ad888f08, 0x3ff9297ba86cc14e),
216        (0xbc97e1e1819d78bd, 0x3ff3dd0025109756),
217        (0xbc6f1b00b6ce8a2a, 0x3fdfded97a03f2f3),
218        (0xbc55487a80e322fe, 0x3fbd5639f2258856),
219        (0x3bede86c7e0323ce, 0x3f8f7261ccd00da5),
220        (0x3bcde1ee8e81b7b5, 0x3f52fbfb5a8c7221),
221        (0x3ba8ed54c3db7fde, 0x3f07d48361a072b6),
222        (0xbb4c30e2cdaa48f2, 0x3eaa8661f0313183),
223        (0xbac575d303dd9b93, 0x3e3205d52df415e6),
224    ];
225
226    let mut p_num = DoubleDouble::mul_f64_add(
227        DoubleDouble::from_bit_pair(P[13]),
228        dx,
229        DoubleDouble::from_bit_pair(P[12]),
230    );
231    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[11]));
232    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[10]));
233    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[9]));
234    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[8]));
235    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[7]));
236    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[6]));
237    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[5]));
238    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[4]));
239    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[3]));
240    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[2]));
241    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[1]));
242    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[0]));
243
244    const Q: [(u64, u64); 14] = [
245        (0x0000000000000000, 0x3ff0000000000000),
246        (0x3cbde8d771be8aba, 0x40118da297250deb),
247        (0xbccdbba67547f5dc, 0x4020589156c4058c),
248        (0x3caca498a3ea822e, 0x4020e9442d441e22),
249        (0xbca6a7d3651d1b42, 0x40156480b1dc22fe),
250        (0x3ca79ba727e70ad6, 0x40013006d1e5d08c),
251        (0x3c8c2b824cfce390, 0x3fe1b0eb2f75de3f),
252        (0xbc5208ab1ad4d8e0, 0x3fb709e145993376),
253        (0xbc21d64934aab809, 0x3f825ee5a8799658),
254        (0x3bc29531f5a4518d, 0x3f40ed532b6bffdd),
255        (0x3b994fb11dace77e, 0x3ef052da92364c6d),
256        (0xbb1f0fb3869f715e, 0x3e8b6bbbe7aae5eb),
257        (0x3a81d8373548a8a8, 0x3e09b5304c6d77ba),
258        (0x3992e1e6b163e1f7, 0xbd50eee2d9d4e7b9),
259    ];
260
261    let mut p_den = DoubleDouble::mul_f64_add(
262        DoubleDouble::from_bit_pair(Q[13]),
263        dx,
264        DoubleDouble::from_bit_pair(Q[12]),
265    );
266    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[11]));
267    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[10]));
268    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[9]));
269    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[8]));
270    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[7]));
271    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[6]));
272    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[5]));
273    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[4]));
274    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[3]));
275    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[2]));
276    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[1]));
277    p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
278    let mut k = DoubleDouble::div(p_num, p_den);
279    k = DoubleDouble::from_exact_add(k.hi, k.lo);
280    let v0 = DoubleDouble::full_dd_sub(k, v_log);
281    apply_sign_and_sum(v0, sum_parity, f_res)
282}
283
284#[cold]
285fn lgamma_4_to_12(dx: f64, f_res: DoubleDouble, sum_parity: f64) -> DoubleDouble {
286    // Log[Gamma[x+1]]
287    // <<FunctionApproximations`
288    // ClearAll["Global`*"]
289    // f[x_]:=LogGamma[x]
290    // {err0,approx}=MiniMaxApproximation[f[z],{z,{4,12},10,10},WorkingPrecision->90]
291    // num=Numerator[approx][[1]];
292    // den=Denominator[approx][[1]];
293    // poly=num;
294    // coeffs=CoefficientList[poly,z];
295    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
296    const P: [(u64, u64); 11] = [
297        (0x3c9a767426b2d1e8, 0x400e45c3573057bb),
298        (0xbcccc22bbae40ac4, 0x403247d3a1b9b0e9),
299        (0xbc82c1d4989e7813, 0x3ffe955db10dd744),
300        (0x3cb4669cf9238f48, 0xc036a67a52498757),
301        (0x3cb4c4039d4da434, 0xc01a8d26ec4b2f7b),
302        (0x3ca43834ea54b437, 0x400c92b79f85b787),
303        (0x3c930746944a1bc1, 0x3ff8b3afe3e0525c),
304        (0xbc68499f7a8b0f87, 0x3fc8059ff9cb3e10),
305        (0x3c0e16d9f08b7e18, 0x3f81c282c77a3862),
306        (0xbbcec78600b42cd0, 0x3f22fe2577cf1017),
307        (0x3b1c222faf24d4a1, 0x3ea5511d126ad883),
308    ];
309
310    let mut p_num = DoubleDouble::mul_f64_add(
311        DoubleDouble::from_bit_pair(P[10]),
312        dx,
313        DoubleDouble::from_bit_pair(P[9]),
314    );
315    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[8]));
316    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[7]));
317    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[6]));
318    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[5]));
319    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[4]));
320    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[3]));
321    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[2]));
322    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[1]));
323    p_num = DoubleDouble::mul_f64_add(p_num, dx, DoubleDouble::from_bit_pair(P[0]));
324
325    const Q: [(u64, u64); 11] = [
326        (0x0000000000000000, 0x3ff0000000000000),
327        (0x3ccaa60b201a29f2, 0x40288b76f18d51d8),
328        (0xbcd831f5042551f9, 0x403d383173e1839d),
329        (0x3cb2b77730673e36, 0x4037e8a6dda469df),
330        (0x3cb8c229012ae276, 0x40207ddd53aa3098),
331        (0xbc8aba961299355d, 0x3ff4c90761849336),
332        (0xbc3e844b2b1f0edb, 0x3fb832c6fba5ff26),
333        (0xbbc70261cd94cb90, 0x3f68b185e16f05c1),
334        (0xbb6c1c2dfe90e592, 0x3f0303509bbb9cf8),
335        (0xbb0f160d6b147ae5, 0x3e7d1bd86b15f52e),
336        (0x3a5714448b9c17d2, 0xbdba82d4d2ccf533),
337    ];
338
339    let mut p_den = DoubleDouble::mul_f64_add(
340        DoubleDouble::from_bit_pair(Q[10]),
341        dx,
342        DoubleDouble::from_bit_pair(Q[9]),
343    );
344    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[8]));
345    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[7]));
346    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[6]));
347    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[5]));
348    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[4]));
349    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[3]));
350    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[2]));
351    p_den = DoubleDouble::mul_f64_add(p_den, dx, DoubleDouble::from_bit_pair(Q[1]));
352    p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
353    apply_sign_and_sum(DoubleDouble::div(p_num, p_den), sum_parity, f_res)
354}
355
356#[cold]
357fn stirling_accurate(dx: f64, parity: f64, f_res: DoubleDouble) -> DoubleDouble {
358    let y_recip = DoubleDouble::from_quick_recip(dx);
359    let y_sqr = DoubleDouble::mult(y_recip, y_recip);
360    // Bernoulli coefficients generated by SageMath:
361    // var('x')
362    // def bernoulli_terms(x, N):
363    //     S = 0
364    //     for k in range(1, N+1):
365    //         B = bernoulli(2*k)
366    //         term = (B / (2*k*(2*k-1))) * x**((2*k-1))
367    //         S += term
368    //     return S
369    //
370    // terms = bernoulli_terms(x, 10)
371    const BERNOULLI_C: [(u64, u64); 10] = [
372        (0x3c55555555555555, 0x3fb5555555555555),
373        (0x3bff49f49f49f49f, 0xbf66c16c16c16c17),
374        (0x3b8a01a01a01a01a, 0x3f4a01a01a01a01a),
375        (0x3befb1fb1fb1fb20, 0xbf43813813813814),
376        (0x3be5c3a9ce01b952, 0x3f4b951e2b18ff23),
377        (0x3bff82553c999b0e, 0xbf5f6ab0d9993c7d),
378        (0x3c10690690690690, 0x3f7a41a41a41a41a),
379        (0x3c21efcdab896745, 0xbf9e4286cb0f5398),
380        (0xbc279e2405a71f88, 0x3fc6fe96381e0680),
381        (0x3c724246319da678, 0xbff6476701181f3a),
382    ];
383    let bernoulli_poly = f_polyeval10(
384        y_sqr,
385        DoubleDouble::from_bit_pair(BERNOULLI_C[0]),
386        DoubleDouble::from_bit_pair(BERNOULLI_C[1]),
387        DoubleDouble::from_bit_pair(BERNOULLI_C[2]),
388        DoubleDouble::from_bit_pair(BERNOULLI_C[3]),
389        DoubleDouble::from_bit_pair(BERNOULLI_C[4]),
390        DoubleDouble::from_bit_pair(BERNOULLI_C[5]),
391        DoubleDouble::from_bit_pair(BERNOULLI_C[6]),
392        DoubleDouble::from_bit_pair(BERNOULLI_C[7]),
393        DoubleDouble::from_bit_pair(BERNOULLI_C[8]),
394        DoubleDouble::from_bit_pair(BERNOULLI_C[9]),
395    );
396
397    // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
398    const LOG2_PI_OVER_2: DoubleDouble =
399        DoubleDouble::from_bit_pair((0xbc865b5a1b7ff5df, 0x3fed67f1c864beb5));
400    let mut log_gamma = DoubleDouble::add(
401        DoubleDouble::mul_add_f64(bernoulli_poly, y_recip, -dx),
402        LOG2_PI_OVER_2,
403    );
404    let dy_log = log_dd(dx);
405    log_gamma = DoubleDouble::mul_add(
406        DoubleDouble::from_exact_add(dy_log.hi, dy_log.lo),
407        DoubleDouble::from_full_exact_add(dx, -0.5),
408        log_gamma,
409    );
410    apply_sign_and_sum(log_gamma, parity, f_res)
411}
412
413/// due to log(x) leading terms cancellation happens around 2,
414/// hence we're using different approximation around LogGamma(2).
415/// Coefficients are ill-conditioned here so Minimal Newton form is used
416#[cold]
417fn lgamma_around_2(x: f64, parity: f64, f_res: DoubleDouble) -> DoubleDouble {
418    let dx = DoubleDouble::from_full_exact_sub(x, 2.);
419
420    // Generated by Wolfram Mathematica:
421    // <<FunctionApproximations`
422    // ClearAll["Global`*"]
423    // f[x_]:=LogGamma[x]
424    // {err0,approx,err1}=MiniMaxApproximation[f[x],{x,{1.95,2.05},11,11},WorkingPrecision->90]
425    // num=Numerator[approx];
426    // den=Denominator[approx];
427    // poly=num;
428    // x0=SetPrecision[2,90];
429    // NumberForm[Series[num,{x,x0,9}],ExponentFunction->(Null&)]
430    // coeffs=Table[SeriesCoefficient[num,{x,x0,k}],{k,0,10}];
431    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
432    // poly=den;
433    // coeffs=Table[SeriesCoefficient[den,{x,x0,k}],{k,0,10}];
434    // NumberForm[Series[den,{x,x0,9}],ExponentFunction->(Null&)]
435    // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
436    const P: [(u64, u64); 10] = [
437        (0xbd8ef0f02676337a, 0x41000f48befb1b70),
438        (0x3dae9e186ab39649, 0x411aeb5de0ba157d),
439        (0xbdcfe4bca8f58298, 0x4122d673e1f69e2c),
440        (0x3db1e3de8e53cfce, 0x411d0fdd168e56a8),
441        (0x3d9114d199dacd01, 0x410b4a1ce996f910),
442        (0x3d9172663132c0b6, 0x40f02d95ceca2c8a),
443        (0x3d52988b7d30cfd3, 0x40c83a4c6e145abc),
444        (0xbd14327346a3db7b, 0x409632fe8dc80419),
445        (0x3cfe9f975e5e9984, 0x4057219509ba1d62),
446        (0xbca0fdd3506bf429, 0x4007988259d795c4),
447    ];
448
449    let dx2 = dx * dx;
450    let dx4 = dx2 * dx2;
451    let dx8 = dx4 * dx4;
452
453    let p0 = DoubleDouble::quick_mul_add(
454        dx,
455        DoubleDouble::from_bit_pair(P[1]),
456        DoubleDouble::from_bit_pair(P[0]),
457    );
458    let p1 = DoubleDouble::quick_mul_add(
459        dx,
460        DoubleDouble::from_bit_pair(P[3]),
461        DoubleDouble::from_bit_pair(P[2]),
462    );
463    let p2 = DoubleDouble::quick_mul_add(
464        dx,
465        DoubleDouble::from_bit_pair(P[5]),
466        DoubleDouble::from_bit_pair(P[4]),
467    );
468    let p3 = DoubleDouble::quick_mul_add(
469        dx,
470        DoubleDouble::from_bit_pair(P[7]),
471        DoubleDouble::from_bit_pair(P[6]),
472    );
473    let p4 = DoubleDouble::quick_mul_add(
474        dx,
475        DoubleDouble::from_bit_pair(P[9]),
476        DoubleDouble::from_bit_pair(P[8]),
477    );
478
479    let q0 = DoubleDouble::quick_mul_add(dx2, p1, p0);
480    let q1 = DoubleDouble::quick_mul_add(dx2, p3, p2);
481
482    let r0 = DoubleDouble::quick_mul_add(dx4, q1, q0);
483    let mut p_num = DoubleDouble::quick_mul_add(dx8, p4, r0);
484    p_num = DoubleDouble::quick_mult(p_num, dx);
485
486    const Q: [(u64, u64); 11] = [
487        (0x3da0f8e36723f959, 0x4112fe27249fc6f1),
488        (0xbdc289e4a571ddb8, 0x412897be1a39b97b),
489        (0xbdb8d18ab489860f, 0x412b4fcbc6d9cb44),
490        (0x3da0550c9f65a5ef, 0x4120fe765c0f6a79),
491        (0xbd90a121e792bf7f, 0x4109fa9eaa0f816a),
492        (0xbd7168de8e78812e, 0x40e9269d002372a5),
493        (0x3d4e4d052cd6982a, 0x40beab7f948d82f0),
494        (0xbd0cebe53b7e81bf, 0x4086a91c01d7241f),
495        (0xbcd2edf097020841, 0x4042a780e9d1b74b),
496        (0xbc73d1a40910845f, 0x3fecd007ff47224f),
497        (0xbc06e4de631047f3, 0x3f7ad97267872491),
498    ];
499
500    let q0 = DoubleDouble::quick_mul_add(
501        dx,
502        DoubleDouble::from_bit_pair(Q[1]),
503        DoubleDouble::from_bit_pair(Q[0]),
504    );
505    let q1 = DoubleDouble::quick_mul_add(
506        dx,
507        DoubleDouble::from_bit_pair(Q[3]),
508        DoubleDouble::from_bit_pair(Q[2]),
509    );
510    let q2 = DoubleDouble::quick_mul_add(
511        dx,
512        DoubleDouble::from_bit_pair(Q[5]),
513        DoubleDouble::from_bit_pair(Q[4]),
514    );
515    let q3 = DoubleDouble::quick_mul_add(
516        dx,
517        DoubleDouble::from_bit_pair(Q[7]),
518        DoubleDouble::from_bit_pair(Q[6]),
519    );
520    let q4 = DoubleDouble::quick_mul_add(
521        dx,
522        DoubleDouble::from_bit_pair(Q[9]),
523        DoubleDouble::from_bit_pair(Q[8]),
524    );
525
526    let r0 = DoubleDouble::quick_mul_add(dx2, q1, q0);
527    let r1 = DoubleDouble::quick_mul_add(dx2, q3, q2);
528
529    let s0 = DoubleDouble::quick_mul_add(dx4, r1, r0);
530    let s1 = DoubleDouble::quick_mul_add(dx2, DoubleDouble::from_bit_pair(Q[10]), q4);
531    let p_den = DoubleDouble::quick_mul_add(dx8, s1, s0);
532    apply_sign_and_sum_quick(DoubleDouble::div(p_num, p_den), parity, f_res)
533}
534
535#[inline]
536pub(crate) fn lgamma_core(x: f64) -> (DoubleDouble, i32) {
537    let ax = f64::from_bits(x.to_bits() & 0x7fff_ffff_ffff_ffff);
538    let dx = ax;
539
540    let is_positive = x.is_sign_positive();
541    let mut sum_parity = 1f64;
542
543    let mut signgam = 1i32;
544
545    if ax < f64::EPSILON {
546        if !is_positive {
547            signgam = -1i32;
548        }
549        return (DoubleDouble::new(0., -f_log(dx)), signgam);
550    }
551
552    let mut f_res = DoubleDouble::default();
553    // For negative x, since (G is gamma function)
554    // -x*G(-x)*G(x) = pi/sin(pi*x),
555    // we have
556    // G(x) = pi/(sin(pi*x)*(-x)*G(-x))
557    // since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
558    // Hence, for x<0, signgam = sign(sin(pi*x)) and
559    // lgamma(x) = log(|Gamma(x)|) = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
560    if !is_positive {
561        let y1 = ax.floor_finite();
562        let fraction = ax - y1; // excess over the boundary
563
564        let a = f_fast_sinpi_dd(fraction);
565
566        sum_parity = -1.;
567        const LOG_PI: DoubleDouble =
568            DoubleDouble::from_bit_pair((0x3c67abf2ad8d5088, 0x3ff250d048e7a1bd));
569        let mut den = DoubleDouble::quick_mult_f64(a, dx);
570        den = DoubleDouble::from_exact_add(den.hi, den.lo);
571        f_res = fast_log_dd(den);
572        f_res = DoubleDouble::from_exact_add(f_res.hi, f_res.lo);
573        f_res = DoubleDouble::quick_dd_sub(LOG_PI, f_res);
574
575        // gamma(x) is negative in (-2n-1,-2n), thus when fx is odd
576        let is_odd = (!is_odd_integer(y1)) as i32;
577        signgam = 1 - (is_odd << 1);
578    }
579
580    if ax <= 0.5 {
581        // Log[Gamma[x + 1]] poly generated by Wolfram
582        // <<FunctionApproximations`
583        // ClearAll["Global`*"]
584        // f[x_]:=LogGamma[x+1]
585        // {err0,approx}=MiniMaxApproximation[f[z],{z,{0.0000000000000001,0.5},7,7},WorkingPrecision->90]
586        // num=Numerator[approx][[1]];
587        // den=Denominator[approx][[1]];
588        // poly=den;
589        // coeffs=CoefficientList[poly,z];
590        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
591        let ps_num = f_polyeval4(
592            dx,
593            f64::from_bits(0x3fea8287cc94dc31),
594            f64::from_bits(0x3fe000cbc75a2ab7),
595            f64::from_bits(0x3fba69bac2495765),
596            f64::from_bits(0x3f78eb3984eb55ee),
597        );
598        let ps_den = f_polyeval5(
599            dx,
600            f64::from_bits(0x3ffee073402b349e),
601            f64::from_bits(0x3fe04c62232d12ec),
602            f64::from_bits(0x3fad09ff23ffb930),
603            f64::from_bits(0x3f5c253f6e8af7d2),
604            f64::from_bits(0xbee18a68b4ed9516),
605        );
606        let mut p_num = DoubleDouble::f64_mul_f64_add(
607            ps_num,
608            dx,
609            DoubleDouble::from_bit_pair((0x3c4d671927a50f5b, 0x3fb184cdffda39b0)),
610        );
611        p_num = DoubleDouble::mul_f64_add(
612            p_num,
613            dx,
614            DoubleDouble::from_bit_pair((0xbc8055e20f9945f5, 0xbfedba6e22cced8a)),
615        );
616        p_num = DoubleDouble::mul_f64_add(
617            p_num,
618            dx,
619            DoubleDouble::from_bit_pair((0x3c56cc8e006896be, 0xbfe2788cfc6fb619)),
620        );
621        p_num = DoubleDouble::mul_f64_add(
622            p_num,
623            dx,
624            DoubleDouble::from_bit_pair((0x34c1e175ecd02e7e, 0xb84ed528d8df1c88)),
625        );
626        let mut p_den = DoubleDouble::f64_mul_f64_add(
627            ps_den,
628            dx,
629            DoubleDouble::from_bit_pair((0xbc70ab0b3b299408, 0x400c16483bffaf57)),
630        );
631        p_den = DoubleDouble::mul_f64_add(
632            p_den,
633            dx,
634            DoubleDouble::from_bit_pair((0xbcac74fbe90fa7cb, 0x400846598b0e4750)),
635        );
636        p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
637        let d_log = fast_log_d_to_dd(dx);
638        let v0 = DoubleDouble::sub(DoubleDouble::div(p_num, p_den), d_log);
639        let l_res = f_res;
640        f_res = apply_sign_and_sum_quick(v0, sum_parity, f_res);
641        let err = f_fmla(
642            f_res.hi.abs(),
643            f64::from_bits(0x3c56a09e667f3bcd), // 2^-57.5
644            f64::from_bits(0x3c20000000000000), // 2^-61
645        );
646        let ub = f_res.hi + (f_res.lo + err);
647        let lb = f_res.hi + (f_res.lo - err);
648        if ub == lb {
649            return (f_res, signgam);
650        }
651        return (lgamma_0p5(dx, d_log, l_res, sum_parity), signgam);
652    } else if ax <= 1. {
653        let distance_to_2 = ax - 2.;
654        if distance_to_2.abs() < 0.05 {
655            return (lgamma_around_2(ax, sum_parity, f_res), signgam);
656        }
657
658        // Log[Gamma[x+1]]
659        // <<FunctionApproximations`
660        // ClearAll["Global`*"]
661        // f[x_]:=LogGamma[x+1]
662        // {err0,approx}=MiniMaxApproximation[f[z],{z,{0.0000000000000001,0.5},9,9},WorkingPrecision->90]
663        // num=Numerator[approx][[1]];
664        // den=Denominator[approx][[1]];
665        // poly=den;
666        // coeffs=CoefficientList[poly,z];
667        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
668
669        const P: [(u64, u64); 10] = [
670            (0xb967217bfcc647c2, 0xbcc9b2c47481ff4f),
671            (0xbc7eff78623354d7, 0xbfe2788cfc6fb552),
672            (0xbc9886c6f42886e0, 0xbff3c6a7f676cf4c),
673            (0x3c77ed956ff8e661, 0xbfdaf57ee2a64253),
674            (0x3c8723f3a5de4fd5, 0x3feb961a3d8bbe89),
675            (0x3c8848ddf2e2726f, 0x3fedc9f11015d4ca),
676            (0xbc799f3b76da571b, 0x3fd7ac6e82c07787),
677            (0xbc5cb131b054a5f5, 0x3fb103e7e288f4da),
678            (0x3bfe93ab961d39a4, 0x3f74789410ab2cf5),
679            (0x3babcb1e8a573475, 0x3f1d74e78621d316),
680        ];
681
682        let ps_den = f_polyeval4(
683            dx,
684            f64::from_bits(0x3fa42e466fd1464e),
685            f64::from_bits(0x3f5f9cd4c60c4e7d),
686            f64::from_bits(0x3efb2d3ed43aab6f),
687            f64::from_bits(0xbe665e3fadf143a6),
688        );
689
690        let dx2 = DoubleDouble::from_exact_mult(dx, dx);
691        let dx4 = DoubleDouble::quick_mult(dx2, dx2);
692        let dx8 = DoubleDouble::quick_mult(dx4, dx4);
693
694        let p0 = DoubleDouble::mul_f64_add(
695            DoubleDouble::from_bit_pair(P[1]),
696            dx,
697            DoubleDouble::from_bit_pair(P[0]),
698        );
699        let p1 = DoubleDouble::mul_f64_add(
700            DoubleDouble::from_bit_pair(P[3]),
701            dx,
702            DoubleDouble::from_bit_pair(P[2]),
703        );
704        let p2 = DoubleDouble::mul_f64_add(
705            DoubleDouble::from_bit_pair(P[5]),
706            dx,
707            DoubleDouble::from_bit_pair(P[4]),
708        );
709        let p3 = DoubleDouble::mul_f64_add(
710            DoubleDouble::from_bit_pair(P[7]),
711            dx,
712            DoubleDouble::from_bit_pair(P[6]),
713        );
714        let p4 = DoubleDouble::mul_f64_add(
715            DoubleDouble::from_bit_pair(P[9]),
716            dx,
717            DoubleDouble::from_bit_pair(P[8]),
718        );
719
720        let q0 = DoubleDouble::mul_add(dx2, p1, p0);
721        let q1 = DoubleDouble::mul_add(dx2, p3, p2);
722
723        let r0 = DoubleDouble::mul_add(dx4, q1, q0);
724        let p_num = DoubleDouble::mul_add(dx8, p4, r0);
725
726        let mut p_den = DoubleDouble::f64_mul_f64_add(
727            ps_den,
728            dx,
729            DoubleDouble::from_bit_pair((0x3c78af0564323160, 0x3fd629441413e902)),
730        );
731        p_den = DoubleDouble::mul_f64_add(
732            p_den,
733            dx,
734            DoubleDouble::from_bit_pair((0xbc901825b170e576, 0x3ff8c99df9c66865)),
735        );
736        p_den = DoubleDouble::mul_f64_add(
737            p_den,
738            dx,
739            DoubleDouble::from_bit_pair((0xbcaf0f1758e16cb3, 0x400e53c84c87f686)),
740        );
741        p_den = DoubleDouble::mul_f64_add(
742            p_den,
743            dx,
744            DoubleDouble::from_bit_pair((0xbc8292764aa01c02, 0x401477d734273be4)),
745        );
746        p_den = DoubleDouble::mul_f64_add(
747            p_den,
748            dx,
749            DoubleDouble::from_bit_pair((0xbca99871cf68ff41, 0x400c87945e972c56)),
750        );
751        p_den = DoubleDouble::mul_f64_add_f64(p_den, dx, f64::from_bits(0x3ff0000000000000));
752        let d_log = fast_log_d_to_dd(dx);
753        let v0 = DoubleDouble::sub(DoubleDouble::div(p_num, p_den), d_log);
754        let l_res = f_res;
755        f_res = apply_sign_and_sum_quick(v0, sum_parity, f_res);
756        let err = f_fmla(
757            f_res.hi.abs(),
758            f64::from_bits(0x3c40000000000000), // 2^-59
759            f64::from_bits(0x3c20000000000000), // 2^-61
760        );
761        let ub = f_res.hi + (f_res.lo + err);
762        let lb = f_res.hi + (f_res.lo - err);
763        if ub == lb {
764            return (f_res, signgam);
765        }
766        return (lgamma_0p5_to_1(dx, d_log, l_res, sum_parity), signgam);
767    } else if ax <= 4. {
768        let distance_to_2 = ax - 2.;
769        if distance_to_2.abs() < 0.05 {
770            return (lgamma_around_2(ax, sum_parity, f_res), signgam);
771        }
772        // <<FunctionApproximations`
773        // ClearAll["Global`*"]
774        // f[x_]:=LogGamma[x+1]
775        // {err0,approx}=MiniMaxApproximation[f[z],{z,{1.0000000000000000001,4},9,9},WorkingPrecision->90]
776        // num=Numerator[approx][[1]];
777        // den=Denominator[approx][[1]];
778        // poly=den;
779        // coeffs=CoefficientList[poly,z];
780        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
781
782        let x2 = DoubleDouble::from_exact_mult(x, x);
783        let x4 = DoubleDouble::quick_mult(x2, x2);
784        let x8 = DoubleDouble::quick_mult(x4, x4);
785
786        const P: [(u64, u64); 10] = [
787            (0x3a8ea8c71173ba1f, 0xbde8c6619bc06d43),
788            (0x3c8f2502d288b7e1, 0xbfe2788cfb0f13f7),
789            (0xbc873b33ddea3333, 0xbfebd290912f0200),
790            (0x3c223d47fd7b2e30, 0x3fb52786e934492b),
791            (0xbc8f4e91d7f48aa5, 0x3fe8204c68bc38f4),
792            (0x3c5356ff82d857c6, 0x3fde676d587374a4),
793            (0x3c5d8deef0e6c21f, 0x3fbef2e284faabe5),
794            (0xbc24ea363b4779fb, 0x3f8bb45183525b51),
795            (0xbbec808b7b332822, 0x3f43b11bd314773b),
796            (0xbb777d551025e6da, 0x3edf7931f7cb9cd1),
797        ];
798
799        let p0 = DoubleDouble::mul_f64_add(
800            DoubleDouble::from_bit_pair(P[1]),
801            dx,
802            DoubleDouble::from_bit_pair(P[0]),
803        );
804        let p1 = DoubleDouble::mul_f64_add(
805            DoubleDouble::from_bit_pair(P[3]),
806            dx,
807            DoubleDouble::from_bit_pair(P[2]),
808        );
809        let p2 = DoubleDouble::mul_f64_add(
810            DoubleDouble::from_bit_pair(P[5]),
811            dx,
812            DoubleDouble::from_bit_pair(P[4]),
813        );
814        let p3 = DoubleDouble::mul_f64_add(
815            DoubleDouble::from_bit_pair(P[7]),
816            dx,
817            DoubleDouble::from_bit_pair(P[6]),
818        );
819        let p4 = DoubleDouble::mul_f64_add(
820            DoubleDouble::from_bit_pair(P[9]),
821            dx,
822            DoubleDouble::from_bit_pair(P[8]),
823        );
824
825        let q0 = DoubleDouble::mul_add(x2, p1, p0);
826        let q1 = DoubleDouble::mul_add(x2, p3, p2);
827
828        let r0 = DoubleDouble::mul_add(x4, q1, q0);
829        let p_num = DoubleDouble::mul_add(x8, p4, r0);
830
831        const Q: [(u64, u64); 10] = [
832            (0x0000000000000000, 0x3ff0000000000000),
833            (0xbc9ad7b53d7da072, 0x4007730c69fb4a20),
834            (0xbc93d2a47740a995, 0x400ab6d03e2d6528),
835            (0x3c9cd643c37f1205, 0x3ffe2cccacb6740b),
836            (0xbc7b616646543538, 0x3fe1f36f793ad9c6),
837            (0xbc483b2cb83a34ba, 0x3fb630083527c66f),
838            (0x3c1089007cac404c, 0x3f7a6b85d9c297ea),
839            (0xbbcfc269fc4a2c55, 0x3f298d01a660c3d9),
840            (0xbb43342127aafe5a, 0x3eb9c8ba657b4b0a),
841            (0xbaac5e3aa213d878, 0xbe14186c41f01fd1),
842        ];
843
844        let p0 = DoubleDouble::mul_f64_add_f64(
845            DoubleDouble::from_bit_pair(Q[1]),
846            dx,
847            f64::from_bits(0x3ff0000000000000),
848        );
849        let p1 = DoubleDouble::mul_f64_add(
850            DoubleDouble::from_bit_pair(Q[3]),
851            dx,
852            DoubleDouble::from_bit_pair(Q[2]),
853        );
854        let p2 = DoubleDouble::mul_f64_add(
855            DoubleDouble::from_bit_pair(Q[5]),
856            dx,
857            DoubleDouble::from_bit_pair(Q[4]),
858        );
859        let p3 = DoubleDouble::mul_f64_add(
860            DoubleDouble::from_bit_pair(Q[7]),
861            dx,
862            DoubleDouble::from_bit_pair(Q[6]),
863        );
864        let p4 = DoubleDouble::f64_mul_f64_add(
865            f64::from_bits(0xbe14186c41f01fd1), // Q[9].hi
866            dx,
867            DoubleDouble::from_bit_pair(Q[8]),
868        );
869
870        let q0 = DoubleDouble::mul_add(x2, p1, p0);
871        let q1 = DoubleDouble::mul_add(x2, p3, p2);
872
873        let r0 = DoubleDouble::mul_add(x4, q1, q0);
874        let p_den = DoubleDouble::mul_add(x8, p4, r0);
875
876        let d_log = fast_log_d_to_dd(dx);
877        let prod = DoubleDouble::div(p_num, p_den);
878        let v0 = DoubleDouble::sub(prod, d_log);
879        let l_res = f_res;
880        f_res = apply_sign_and_sum_quick(v0, sum_parity, f_res);
881        let err = f_fmla(
882            f_res.hi.abs(),
883            f64::from_bits(0x3c40000000000000), // 2^-59
884            f64::from_bits(0x3c00000000000000), // 2^-63
885        );
886        let ub = f_res.hi + (f_res.lo + err);
887        let lb = f_res.hi + (f_res.lo - err);
888        if ub == lb {
889            return (f_res, signgam);
890        }
891        return (lgamma_1_to_4(dx, d_log, l_res, sum_parity), signgam);
892    } else if ax <= 12. {
893        // <<FunctionApproximations`
894        // ClearAll["Global`*"]
895        // f[x_]:=LogGamma[x]
896        // {err0,approx}=MiniMaxApproximation[f[z],{z,{4,12},9,9},WorkingPrecision->90]
897        // num=Numerator[approx][[1]];
898        // den=Denominator[approx][[1]];
899        // poly=den;
900        // coeffs=CoefficientList[poly,z];
901        // TableForm[Table[Row[{"'",NumberForm[coeffs[[i+1]],{50,50},ExponentFunction->(Null&)],"',"}],{i,0,Length[coeffs]-1}]]
902
903        let x2 = DoubleDouble::from_exact_mult(x, x);
904        let x4 = DoubleDouble::quick_mult(x2, x2);
905        let x8 = DoubleDouble::quick_mult(x4, x4);
906
907        const P: [(u64, u64); 10] = [
908            (0x3ca9a6c909e67304, 0x400c83f8e5e68934),
909            (0xbcc3a71a296d1f00, 0x40278d8a2abd6aec),
910            (0xbca623fa8857b35a, 0xc0144dd0190486d6),
911            (0x3cc1845532bca122, 0xc02920aaae63c5a7),
912            (0x3c57111721fd9df2, 0xbfc9a27952cac38f),
913            (0xbca78a77f8acae38, 0x400043078ac20503),
914            (0xbc721a88d770af7e, 0x3fdbeba4a1a95bfd),
915            (0xbc09c9e5917a665e, 0x3f9e18ff2504fd11),
916            (0xbbeb6e5c1cdf8c87, 0x3f45b0595f7eb903),
917            (0xbaf149d407d419d3, 0x3ecf5336ddb96b5f),
918        ];
919
920        let p0 = DoubleDouble::mul_f64_add(
921            DoubleDouble::from_bit_pair(P[1]),
922            dx,
923            DoubleDouble::from_bit_pair(P[0]),
924        );
925        let p1 = DoubleDouble::mul_f64_add(
926            DoubleDouble::from_bit_pair(P[3]),
927            dx,
928            DoubleDouble::from_bit_pair(P[2]),
929        );
930        let p2 = DoubleDouble::mul_f64_add(
931            DoubleDouble::from_bit_pair(P[5]),
932            dx,
933            DoubleDouble::from_bit_pair(P[4]),
934        );
935        let p3 = DoubleDouble::mul_f64_add(
936            DoubleDouble::from_bit_pair(P[7]),
937            dx,
938            DoubleDouble::from_bit_pair(P[6]),
939        );
940        let p4 = DoubleDouble::mul_f64_add(
941            DoubleDouble::from_bit_pair(P[9]),
942            dx,
943            DoubleDouble::from_bit_pair(P[8]),
944        );
945
946        let q0 = DoubleDouble::mul_add(x2, p1, p0);
947        let q1 = DoubleDouble::mul_add(x2, p3, p2);
948
949        let r0 = DoubleDouble::mul_add(x4, q1, q0);
950        let p_num = DoubleDouble::mul_add(x8, p4, r0);
951
952        const Q: [(u64, u64); 10] = [
953            (0x0000000000000000, 0x3ff0000000000000),
954            (0x3cc5e0cfc2a3ed2f, 0x4023561fd4efbbb2),
955            (0x3c829f67da778215, 0x403167ff0d04f99a),
956            (0x3ca3c8e7b81b165f, 0x4024f2d3c7d9439f),
957            (0xbca90199265d1bfc, 0x40045530db97bad5),
958            (0xbc3e89169d10977f, 0x3fd0d9f46084a388),
959            (0x3c1c2b7582566435, 0x3f872f7f248227bd),
960            (0x3ba8f3f0294e144f, 0x3f271e198971c58e),
961            (0x3b4d13ceca0a9bf7, 0x3ea62e85cd267c65),
962            (0xba896f9fc0c4f644, 0xbde99e5ee24ff6ba),
963        ];
964
965        let p0 = DoubleDouble::mul_f64_add(
966            DoubleDouble::from_bit_pair(Q[1]),
967            dx,
968            DoubleDouble::from_bit_pair(Q[0]),
969        );
970        let p1 = DoubleDouble::mul_f64_add(
971            DoubleDouble::from_bit_pair(Q[3]),
972            dx,
973            DoubleDouble::from_bit_pair(Q[2]),
974        );
975        let p2 = DoubleDouble::mul_f64_add(
976            DoubleDouble::from_bit_pair(Q[5]),
977            dx,
978            DoubleDouble::from_bit_pair(Q[4]),
979        );
980        let p3 = DoubleDouble::mul_f64_add(
981            DoubleDouble::from_bit_pair(Q[7]),
982            dx,
983            DoubleDouble::from_bit_pair(Q[6]),
984        );
985        let p4 = DoubleDouble::f64_mul_f64_add(
986            f64::from_bits(0xbde99e5ee24ff6ba), //Q[9].hi
987            dx,
988            DoubleDouble::from_bit_pair(Q[8]),
989        );
990
991        let q0 = DoubleDouble::mul_add(x2, p1, p0);
992        let q1 = DoubleDouble::mul_add(x2, p3, p2);
993
994        let r0 = DoubleDouble::mul_add(x4, q1, q0);
995        let p_den = DoubleDouble::mul_add(x8, p4, r0);
996
997        let l_res = f_res;
998        f_res = apply_sign_and_sum_quick(DoubleDouble::div(p_num, p_den), sum_parity, f_res);
999        let err = f_fmla(
1000            f_res.hi.abs(),
1001            f64::from_bits(0x3c50000000000000), // 2^-58
1002            f64::from_bits(0x3bd0000000000000), // 2^-66
1003        );
1004        let ub = f_res.hi + (f_res.lo + err);
1005        let lb = f_res.hi + (f_res.lo - err);
1006        if ub == lb {
1007            return (f_res, signgam);
1008        }
1009        return (lgamma_4_to_12(dx, l_res, sum_parity), signgam);
1010    }
1011    // Stirling's approximation of Log(Gamma) and then Exp[Log[Gamma]]
1012    let y_recip = DoubleDouble::from_quick_recip(dx);
1013    let y_sqr = DoubleDouble::mult(y_recip, y_recip);
1014    // Bernoulli coefficients generated by SageMath:
1015    // var('x')
1016    // def bernoulli_terms(x, N):
1017    //     S = 0
1018    //     for k in range(1, N+1):
1019    //         B = bernoulli(2*k)
1020    //         term = (B / (2*k*(2*k-1))) * x**((2*k-1))
1021    //         S += term
1022    //     return S
1023    //
1024    // terms = bernoulli_terms(x, 7)
1025    let bernoulli_poly_s = f_polyeval6(
1026        y_sqr.hi,
1027        f64::from_bits(0xbf66c16c16c16c17),
1028        f64::from_bits(0x3f4a01a01a01a01a),
1029        f64::from_bits(0xbf43813813813814),
1030        f64::from_bits(0x3f4b951e2b18ff23),
1031        f64::from_bits(0xbf5f6ab0d9993c7d),
1032        f64::from_bits(0x3f7a41a41a41a41a),
1033    );
1034    let bernoulli_poly = DoubleDouble::mul_f64_add(
1035        y_sqr,
1036        bernoulli_poly_s,
1037        DoubleDouble::from_bit_pair((0x3c55555555555555, 0x3fb5555555555555)),
1038    );
1039    // Log[Gamma(x)] = x*log(x) - x + 1/2*Log(2*PI/x) + bernoulli_terms
1040    const LOG2_PI_OVER_2: DoubleDouble =
1041        DoubleDouble::from_bit_pair((0xbc865b5a1b7ff5df, 0x3fed67f1c864beb5));
1042    let mut log_gamma = DoubleDouble::add(
1043        DoubleDouble::mul_add_f64(bernoulli_poly, y_recip, -dx),
1044        LOG2_PI_OVER_2,
1045    );
1046    let dy_log = fast_log_d_to_dd(dx);
1047    log_gamma = DoubleDouble::mul_add(
1048        DoubleDouble::from_exact_add(dy_log.hi, dy_log.lo),
1049        DoubleDouble::from_full_exact_add(dx, -0.5),
1050        log_gamma,
1051    );
1052    let l_res = f_res;
1053    f_res = apply_sign_and_sum_quick(log_gamma, sum_parity, f_res);
1054    let err = f_fmla(
1055        f_res.hi.abs(),
1056        f64::from_bits(0x3c30000000000000), // 2^-60
1057        f64::from_bits(0x3bc0000000000000), // 2^-67
1058    );
1059    let ub = f_res.hi + (f_res.lo + err);
1060    let lb = f_res.hi + (f_res.lo - err);
1061    if ub == lb {
1062        return (f_res, signgam);
1063    }
1064    (stirling_accurate(dx, sum_parity, l_res), signgam)
1065}
1066
1067/// Computes log(gamma(x))
1068///
1069/// ulp 0.52
1070pub fn f_lgamma(x: f64) -> f64 {
1071    let nx = x.to_bits().wrapping_shl(1);
1072    if nx >= 0xfeaea9b24f16a34cu64 || nx == 0 {
1073        // |x| >= 0x1.006df1bfac84ep+1015
1074        if x.is_infinite() {
1075            return f64::INFINITY;
1076        }
1077        if x.is_nan() {
1078            return f64::NAN;
1079        }
1080        return f64::INFINITY;
1081    }
1082
1083    if is_integer(x) {
1084        if x == 2. || x == 1. {
1085            return 0.;
1086        }
1087        if x.is_sign_negative() {
1088            return f64::INFINITY;
1089        }
1090    }
1091
1092    lgamma_core(x).0.to_f64()
1093}
1094
1095#[cfg(test)]
1096mod tests {
1097    use super::*;
1098    #[test]
1099    fn test_lgamma() {
1100        assert_eq!(f_lgamma(0.), f64::INFINITY);
1101        assert_eq!(f_lgamma(-0.), f64::INFINITY);
1102        assert_eq!(f_lgamma(-4.039410591125488), -0.001305303022594149);
1103        assert_eq!(f_lgamma(-2.000001907348633), 12.47664749001284);
1104        assert_eq!(
1105            f_lgamma(0.0000000000000006939032951805219),
1106            34.90419906721559
1107        );
1108        assert_eq!(f_lgamma(0.9743590354901843), 0.01534797880086699);
1109        assert_eq!(f_lgamma(1.9533844296518055), -0.019000687007583488);
1110        assert_eq!(f_lgamma(1.9614259600725743), -0.015824770893504085);
1111        assert_eq!(f_lgamma(1.961426168688831), -0.015824687947423532);
1112        assert_eq!(f_lgamma(2.0000000026484486), 0.0000000011197225706325235);
1113        assert_eq!(f_lgamma(f64::INFINITY), f64::INFINITY);
1114        assert!(f_lgamma(f64::NAN).is_nan());
1115        // assert_eq!(f_lgamma(1.2902249255008019e301),
1116        //            8932652024571557000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.);
1117        assert_eq!(
1118            f_lgamma(2.000000000000014),
1119            0.000000000000006008126761947661
1120        );
1121        assert_eq!(f_lgamma(-2.74999999558122), 0.004487888879321723);
1122        assert_eq!(
1123            f_lgamma(0.00000000000001872488985570349),
1124            31.608922747730112
1125        );
1126        assert_eq!(f_lgamma(-2.7484742253727745), 0.0015213685011468314);
1127        assert_eq!(f_lgamma(0.9759521409866919), 0.014362097695996162);
1128    }
1129}
pxfm/gamma/lgamma.rs

pxfm/gamma/
lgamma.rs
