1 /* mpfr_tlgamma -- test file for lgamma function
2
3 Copyright 2005-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23 #include "mpfr-test.h"
24
25 static int
26 mpfr_lgamma_nosign (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd)
27 {
28 int inex, sign;
29
30 inex = mpfr_lgamma (y, &sign, x, rnd);
31 if (!MPFR_IS_SINGULAR (y))
32 {
33 MPFR_ASSERTN (sign == 1 || sign == -1);
34 if (sign == -1 && (rnd == MPFR_RNDN || rnd == MPFR_RNDZ))
35 {
36 mpfr_neg (y, y, MPFR_RNDN);
37 inex = -inex;
38 /* This is a way to check with the generic tests, that the value
39 returned in the sign variable is consistent, but warning! The
40 tested function depends on the rounding mode: it is
41 * lgamma(x) = log(|Gamma(x)|) in MPFR_RNDD and MPFR_RNDU;
42 * lgamma(x) * sign(Gamma(x)) in MPFR_RNDN and MPFR_RNDZ. */
43 }
44 }
45
46 return inex;
47 }
48
49 #define TEST_FUNCTION mpfr_lgamma_nosign
50 #include "tgeneric.c"
51
52 static void
53 special (void)
54 {
55 mpfr_t x, y;
56 int inex;
57 int sign;
58 mpfr_exp_t emin, emax;
59
60 mpfr_init (x);
61 mpfr_init (y);
62
63 mpfr_set_nan (x);
64 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
65 if (!mpfr_nan_p (y))
66 {
67 printf ("Error for lgamma(NaN)\n");
68 exit (1);
69 }
70
71 mpfr_set_inf (x, -1);
72 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
73 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
74 {
75 printf ("Error for lgamma(-Inf)\n");
76 exit (1);
77 }
78
79 mpfr_set_inf (x, 1);
80 sign = -17;
81 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
82 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1)
83 {
84 printf ("Error for lgamma(+Inf)\n");
85 exit (1);
86 }
87
88 mpfr_set_ui (x, 0, MPFR_RNDN);
89 sign = -17;
90 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
91 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != 1)
92 {
93 printf ("Error for lgamma(+0)\n");
94 exit (1);
95 }
96
97 mpfr_set_ui (x, 0, MPFR_RNDN);
98 mpfr_neg (x, x, MPFR_RNDN);
99 sign = -17;
100 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
101 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || sign != -1)
102 {
103 printf ("Error for lgamma(-0)\n");
104 exit (1);
105 }
106
107 mpfr_set_ui (x, 1, MPFR_RNDN);
108 sign = -17;
109 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
110 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1)
111 {
112 printf ("Error for lgamma(1)\n");
113 exit (1);
114 }
115
116 mpfr_set_si (x, -1, MPFR_RNDN);
117 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
118 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
119 {
120 printf ("Error for lgamma(-1)\n");
121 exit (1);
122 }
123
124 mpfr_set_ui (x, 2, MPFR_RNDN);
125 sign = -17;
126 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
127 if (MPFR_IS_NAN (y) || mpfr_cmp_ui (y, 0) || MPFR_IS_NEG (y) || sign != 1)
128 {
129 printf ("Error for lgamma(2)\n");
130 exit (1);
131 }
132
133 mpfr_set_prec (x, 53);
134 mpfr_set_prec (y, 53);
135
136 #define CHECK_X1 "1.0762904832837976166"
137 #define CHECK_Y1 "-0.039418362817587634939"
138
139 mpfr_set_str (x, CHECK_X1, 10, MPFR_RNDN);
140 sign = -17;
141 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
142 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN);
143 if (mpfr_equal_p (y, x) == 0 || sign != 1)
144 {
145 printf ("mpfr_lgamma(" CHECK_X1 ") is wrong:\n"
146 "expected ");
147 mpfr_dump (x);
148 printf ("got ");
149 mpfr_dump (y);
150 exit (1);
151 }
152
153 #define CHECK_X2 "9.23709516716202383435e-01"
154 #define CHECK_Y2 "0.049010669407893718563"
155 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
156 sign = -17;
157 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
158 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
159 if (mpfr_equal_p (y, x) == 0 || sign != 1)
160 {
161 printf ("mpfr_lgamma(" CHECK_X2 ") is wrong:\n"
162 "expected ");
163 mpfr_dump (x);
164 printf ("got ");
165 mpfr_dump (y);
166 exit (1);
167 }
168
169 mpfr_set_prec (x, 8);
170 mpfr_set_prec (y, 175);
171 mpfr_set_ui (x, 33, MPFR_RNDN);
172 sign = -17;
173 mpfr_lgamma (y, &sign, x, MPFR_RNDU);
174 mpfr_set_prec (x, 175);
175 mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7");
176 if (mpfr_equal_p (x, y) == 0 || sign != 1)
177 {
178 printf ("Error in mpfr_lgamma (1)\n");
179 exit (1);
180 }
181
182 mpfr_set_prec (x, 21);
183 mpfr_set_prec (y, 8);
184 mpfr_set_ui (y, 120, MPFR_RNDN);
185 sign = -17;
186 mpfr_lgamma (x, &sign, y, MPFR_RNDZ);
187 mpfr_set_prec (y, 21);
188 mpfr_set_str_binary (y, "0.111000101000001100101E9");
189 if (mpfr_equal_p (x, y) == 0 || sign != 1)
190 {
191 printf ("Error in mpfr_lgamma (120)\n");
192 printf ("Expected "); mpfr_dump (y);
193 printf ("Got "); mpfr_dump (x);
194 exit (1);
195 }
196
197 mpfr_set_prec (x, 3);
198 mpfr_set_prec (y, 206);
199 mpfr_set_str_binary (x, "0.110e10");
200 sign = -17;
201 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN);
202 mpfr_set_prec (x, 206);
203 mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13");
204 if (mpfr_equal_p (x, y) == 0 || sign != 1)
205 {
206 printf ("Error in mpfr_lgamma (768)\n");
207 exit (1);
208 }
209 if (inex >= 0)
210 {
211 printf ("Wrong flag for mpfr_lgamma (768)\n");
212 exit (1);
213 }
214
215 mpfr_set_prec (x, 4);
216 mpfr_set_prec (y, 4);
217 mpfr_set_str_binary (x, "0.1100E-66");
218 sign = -17;
219 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
220 mpfr_set_str_binary (x, "0.1100E6");
221 if (mpfr_equal_p (x, y) == 0 || sign != 1)
222 {
223 printf ("Error for lgamma(0.1100E-66)\n");
224 printf ("Expected ");
225 mpfr_dump (x);
226 printf ("Got ");
227 mpfr_dump (y);
228 exit (1);
229 }
230
231 mpfr_set_prec (x, 256);
232 mpfr_set_prec (y, 32);
233 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN);
234 mpfr_add_ui (x, x, 1, MPFR_RNDN);
235 mpfr_div_2ui (x, x, 1, MPFR_RNDN);
236 sign = -17;
237 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
238 mpfr_set_prec (x, 32);
239 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207");
240 if (mpfr_equal_p (x, y) == 0 || sign != 1)
241 {
242 printf ("Error for lgamma(-2^199+0.5)\n");
243 printf ("Got ");
244 mpfr_dump (y);
245 printf ("instead of ");
246 mpfr_dump (x);
247 exit (1);
248 }
249
250 mpfr_set_prec (x, 256);
251 mpfr_set_prec (y, 32);
252 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN);
253 mpfr_sub_ui (x, x, 1, MPFR_RNDN);
254 mpfr_div_2ui (x, x, 1, MPFR_RNDN);
255 sign = -17;
256 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
257 mpfr_set_prec (x, 32);
258 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207");
259 if (mpfr_equal_p (x, y) == 0 || sign != -1)
260 {
261 printf ("Error for lgamma(-2^199-0.5)\n");
262 printf ("Got ");
263 mpfr_dump (y);
264 printf ("with sign %d instead of ", sign);
265 mpfr_dump (x);
266 printf ("with sign -1.\n");
267 exit (1);
268 }
269
270 mpfr_set_prec (x, 10);
271 mpfr_set_prec (y, 10);
272 mpfr_set_str_binary (x, "-0.1101111000E-3");
273 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN);
274 mpfr_set_str_binary (x, "10.01001011");
275 if (mpfr_equal_p (x, y) == 0 || sign != -1 || inex >= 0)
276 {
277 printf ("Error for lgamma(-0.1101111000E-3)\n");
278 printf ("Got ");
279 mpfr_dump (y);
280 printf ("instead of ");
281 mpfr_dump (x);
282 printf ("with sign %d instead of -1 (inex=%d).\n", sign, inex);
283 exit (1);
284 }
285
286 mpfr_set_prec (x, 18);
287 mpfr_set_prec (y, 28);
288 mpfr_set_str_binary (x, "-1.10001101010001101e-196");
289 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN);
290 mpfr_set_prec (x, 28);
291 mpfr_set_str_binary (x, "0.100001110110101011011010011E8");
292 MPFR_ASSERTN (mpfr_equal_p (x, y) && inex < 0);
293
294 /* values reported by Kaveh Ghazi on 14 Jul 2007, where mpfr_lgamma()
295 takes forever */
296 #define VAL1 "-0.11100001001010110111001010001001001011110100110000110E-55"
297 #define OUT1 "100110.01000000010111001110110101110101001001100110111"
298 #define VAL2 "-0.11100001001010110111001010001001001011110011111111100E-55"
299 #define OUT2 "100110.0100000001011100111011010111010100100110011111"
300 #define VAL3 "-0.11100001001010110111001010001001001001110101101010100E-55"
301 #define OUT3 "100110.01000000010111001110110101110101001011110111011"
302 #define VAL4 "-0.10001111110110110100100100000000001111110001001001011E-57"
303 #define OUT4 "101000.0001010111110011101101000101111111010001100011"
304 #define VAL5 "-0.10001111110110110100100100000000001111011111100001000E-57"
305 #define OUT5 "101000.00010101111100111011010001011111110100111000001"
306 #define VAL6 "-0.10001111110110110100100100000000001111011101100011001E-57"
307 #define OUT6 "101000.0001010111110011101101000101111111010011101111"
308
309 mpfr_set_prec (x, 53);
310 mpfr_set_prec (y, 53);
311
312 mpfr_set_str_binary (x, VAL1);
313 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
314 mpfr_set_str_binary (x, OUT1);
315 MPFR_ASSERTN(sign == -1 && mpfr_equal_p(x, y));
316
317 mpfr_set_str_binary (x, VAL2);
318 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
319 mpfr_set_str_binary (x, OUT2);
320 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
321
322 mpfr_set_str_binary (x, VAL3);
323 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
324 mpfr_set_str_binary (x, OUT3);
325 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
326
327 mpfr_set_str_binary (x, VAL4);
328 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
329 mpfr_set_str_binary (x, OUT4);
330 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
331
332 mpfr_set_str_binary (x, VAL5);
333 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
334 mpfr_set_str_binary (x, OUT5);
335 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
336
337 mpfr_set_str_binary (x, VAL6);
338 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
339 mpfr_set_str_binary (x, OUT6);
340 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
341
342 /* further test from Kaveh Ghazi */
343 mpfr_set_str_binary (x, "-0.10011010101001010010001110010111010111011101010111001E-53");
344 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
345 mpfr_set_str_binary (x, "100101.00111101101010000000101010111010001111001101111");
346 MPFR_ASSERTN(sign == -1 && mpfr_equal_p (x, y));
347
348 /* bug found by Kevin Rauch on 26 Oct 2007 */
349 emin = mpfr_get_emin ();
350 emax = mpfr_get_emax ();
351 set_emin (-1000000000);
352 set_emax (1000000000);
353 mpfr_set_ui (x, 1, MPFR_RNDN);
354 mpfr_lgamma (x, &sign, x, MPFR_RNDN);
355 MPFR_ASSERTN(mpfr_get_emin () == -1000000000);
356 MPFR_ASSERTN(mpfr_get_emax () == 1000000000);
357 set_emin (emin);
358 set_emax (emax);
359
360 /* two other bugs reported by Kevin Rauch on 27 Oct 2007 */
361 mpfr_set_prec (x, 128);
362 mpfr_set_prec (y, 128);
363 mpfr_set_str_binary (x, "0.11000110011110111111110010100110000000000000000000000000000000000000000000000000000000000000000001000011000110100100110111101010E-765689");
364 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN);
365 mpfr_set_str_binary (x, "10000001100100101111011011010000111010001001110000111010011000101001011111011111110011011010110100101111110111001001010100011101E-108");
366 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0);
367
368 mpfr_set_prec (x, 128);
369 mpfr_set_prec (y, 256);
370 mpfr_set_str_binary (x, "0.1011111111111111100000111011111E-31871");
371 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDN);
372 mpfr_set_prec (x, 256);
373 mpfr_set_str (x, "AC9729B83707E6797612D0D76DAF42B1240A677FF1B6E3783FD4E53037143B1P-237", 16, MPFR_RNDN);
374 MPFR_ASSERTN(inex < 0 && mpfr_cmp (y, x) == 0 && sign > 0);
375
376 mpfr_clear (x);
377 mpfr_clear (y);
378 }
379
380 static int
381 mpfr_lgamma1 (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r)
382 {
383 int sign;
384
385 return mpfr_lgamma (y, &sign, x, r);
386 }
387
388 /* Since r10377, the following test causes a "too much memory" error
389 when MPFR is built with Debian's tcc 0.9.27~git20151227.933c223-1
390 on x86_64. The problem came from __gmpfr_ceil_log2, now fixed in
391 r10443 (according to the integer promotion rules, this appeared to
392 be a bug in tcc, not in MPFR; however, relying on such an obscure
393 rule was not a good idea). */
394 static void
395 tcc_bug20160606 (void)
396 {
397 mpfr_t x, y;
398 int sign;
399
400 mpfr_init2 (x, 53);
401 mpfr_init2 (y, 53);
402 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN);
403 mpfr_lgamma (y, &sign, x, MPFR_RNDN);
404 mpfr_clear (x);
405 mpfr_clear (y);
406 }
407
408 /* With r12088, mpfr_lgamma is much too slow with a reduced emax that
409 yields an overflow, even though this case is easier. In practice,
410 this test will hang. */
411 static void
412 bug20180110 (void)
413 {
414 mpfr_exp_t emax, e;
415 mpfr_t x, y, z;
416 mpfr_flags_t flags, eflags;
417 int i, inex, sign;
418
419 emax = mpfr_get_emax ();
420
421 mpfr_init2 (x, 2);
422 mpfr_inits2 (64, y, z, (mpfr_ptr) 0);
423 eflags = MPFR_FLAGS_INEXACT | MPFR_FLAGS_OVERFLOW;
424
425 for (i = 10; i <= 30; i++)
426 {
427 mpfr_set_si_2exp (x, -1, -(1L << i), MPFR_RNDN); /* -2^(-2^i) */
428 mpfr_lgamma (y, &sign, x, MPFR_RNDZ);
429 e = mpfr_get_exp (y);
430 set_emax (e - 1);
431 mpfr_clear_flags ();
432 inex = mpfr_lgamma (y, &sign, x, MPFR_RNDZ);
433 flags = __gmpfr_flags;
434 mpfr_set_inf (z, 1);
435 mpfr_nextbelow (z);
436 set_emax (emax);
437 if (! (mpfr_equal_p (y, z) && SAME_SIGN (inex, -1) && flags == eflags))
438 {
439 printf ("Error in bug20180110 for i = %d:\n", i);
440 printf ("Expected ");
441 mpfr_dump (z);
442 printf ("with inex = %d and flags =", -1);
443 flags_out (eflags);
444 printf ("Got ");
445 mpfr_dump (y);
446 printf ("with inex = %d and flags =", inex);
447 flags_out (flags);
448 exit (1);
449 }
450 }
451
452 mpfr_clears (x, y, z, (mpfr_ptr) 0);
453 }
454
455 int
456 main (void)
457 {
458 tests_start_mpfr ();
459
460 tcc_bug20160606 ();
461 bug20180110 ();
462
463 special ();
464 test_generic (MPFR_PREC_MIN, 100, 2);
465
466 data_check ("data/lgamma", mpfr_lgamma1, "mpfr_lgamma");
467
468 tests_end_mpfr ();
469 return 0;
470 }