1 /* mpfr_tlngamma -- test file for lngamma 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 #define TEST_FUNCTION mpfr_lngamma
26 #define TEST_RANDOM_POS 16
27 #include "tgeneric.c"
28
29 static void
30 special (void)
31 {
32 mpfr_t x, y;
33 int i, inex;
34
35 mpfr_init (x);
36 mpfr_init (y);
37
38 mpfr_set_nan (x);
39 mpfr_lngamma (y, x, MPFR_RNDN);
40 if (!mpfr_nan_p (y))
41 {
42 printf ("Error for lngamma(NaN)\n");
43 exit (1);
44 }
45
46 mpfr_set_inf (x, 1);
47 mpfr_clear_flags ();
48 mpfr_lngamma (y, x, MPFR_RNDN);
49 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || __gmpfr_flags != 0)
50 {
51 printf ("Error for lngamma(+Inf)\n");
52 exit (1);
53 }
54
55 mpfr_set_inf (x, -1);
56 mpfr_clear_flags ();
57 mpfr_lngamma (y, x, MPFR_RNDN);
58 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 || __gmpfr_flags != 0)
59 {
60 printf ("Error for lngamma(-Inf)\n");
61 exit (1);
62 }
63
64 mpfr_set_ui (x, 0, MPFR_RNDN);
65 mpfr_clear_flags ();
66 mpfr_lngamma (y, x, MPFR_RNDN);
67 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
68 __gmpfr_flags != MPFR_FLAGS_DIVBY0)
69 {
70 printf ("Error for lngamma(+0)\n");
71 exit (1);
72 }
73
74 mpfr_set_ui (x, 0, MPFR_RNDN);
75 mpfr_neg (x, x, MPFR_RNDN);
76 mpfr_clear_flags ();
77 mpfr_lngamma (y, x, MPFR_RNDN);
78 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
79 __gmpfr_flags != MPFR_FLAGS_DIVBY0)
80 {
81 printf ("Error for lngamma(-0)\n");
82 exit (1);
83 }
84
85 mpfr_set_ui (x, 1, MPFR_RNDN);
86 mpfr_clear_flags ();
87 mpfr_lngamma (y, x, MPFR_RNDN);
88 if (mpfr_cmp_ui0 (y, 0) || MPFR_IS_NEG (y))
89 {
90 printf ("Error for lngamma(1)\n");
91 exit (1);
92 }
93
94 for (i = 1; i <= 5; i++)
95 {
96 int c;
97
98 mpfr_set_si (x, -i, MPFR_RNDN);
99 mpfr_clear_flags ();
100 mpfr_lngamma (y, x, MPFR_RNDN);
101 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0 ||
102 __gmpfr_flags != MPFR_FLAGS_DIVBY0)
103 {
104 printf ("Error for lngamma(-%d)\n", i);
105 exit (1);
106 }
107 if (i & 1)
108 {
109 mpfr_nextabove (x);
110 c = '+';
111 }
112 else
113 {
114 mpfr_nextbelow (x);
115 c = '-';
116 }
117 mpfr_lngamma (y, x, MPFR_RNDN);
118 if (!mpfr_nan_p (y))
119 {
120 printf ("Error for lngamma(-%d%cepsilon)\n", i, c);
121 exit (1);
122 }
123 }
124
125 mpfr_set_ui (x, 2, MPFR_RNDN);
126 mpfr_lngamma (y, x, MPFR_RNDN);
127 if (mpfr_cmp_ui0 (y, 0) || MPFR_IS_NEG (y))
128 {
129 printf ("Error for lngamma(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 mpfr_lngamma (y, x, MPFR_RNDN);
141 mpfr_set_str (x, CHECK_Y1, 10, MPFR_RNDN);
142 if (MPFR_IS_NAN (y) || mpfr_cmp (y, x))
143 {
144 printf ("mpfr_lngamma(" CHECK_X1 ") is wrong:\n"
145 "expected ");
146 mpfr_dump (x);
147 printf ("got ");
148 mpfr_dump (y);
149 exit (1);
150 }
151
152 #define CHECK_X2 "9.23709516716202383435e-01"
153 #define CHECK_Y2 "0.049010669407893718563"
154 mpfr_set_str (x, CHECK_X2, 10, MPFR_RNDN);
155 mpfr_lngamma (y, x, MPFR_RNDN);
156 mpfr_set_str (x, CHECK_Y2, 10, MPFR_RNDN);
157 if (mpfr_cmp0 (y, x))
158 {
159 printf ("mpfr_lngamma(" CHECK_X2 ") is wrong:\n"
160 "expected ");
161 mpfr_dump (x);
162 printf ("got ");
163 mpfr_dump (y);
164 exit (1);
165 }
166
167 mpfr_set_prec (x, 8);
168 mpfr_set_prec (y, 175);
169 mpfr_set_ui (x, 33, MPFR_RNDN);
170 mpfr_lngamma (y, x, MPFR_RNDU);
171 mpfr_set_prec (x, 175);
172 mpfr_set_str_binary (x, "0.1010001100011101101011001101110010100001000001000001110011000001101100001111001001000101011011100100010101011110100111110101010100010011010010000101010111001100011000101111E7");
173 if (mpfr_cmp0 (x, y))
174 {
175 printf ("Error in mpfr_lngamma (1)\n");
176 exit (1);
177 }
178
179 mpfr_set_prec (x, 21);
180 mpfr_set_prec (y, 8);
181 mpfr_set_ui (y, 120, MPFR_RNDN);
182 mpfr_lngamma (x, y, MPFR_RNDZ);
183 mpfr_set_prec (y, 21);
184 mpfr_set_str_binary (y, "0.111000101000001100101E9");
185 if (mpfr_cmp0 (x, y))
186 {
187 printf ("Error in mpfr_lngamma (120)\n");
188 printf ("Expected "); mpfr_dump (y);
189 printf ("Got "); mpfr_dump (x);
190 exit (1);
191 }
192
193 mpfr_set_prec (x, 3);
194 mpfr_set_prec (y, 206);
195 mpfr_set_str_binary (x, "0.110e10");
196 inex = mpfr_lngamma (y, x, MPFR_RNDN);
197 mpfr_set_prec (x, 206);
198 mpfr_set_str_binary (x, "0.10000111011000000011100010101001100110001110000111100011000100100110110010001011011110101001111011110110000001010100111011010000000011100110110101100111000111010011110010000100010111101010001101000110101001E13");
199 if (mpfr_cmp0 (x, y))
200 {
201 printf ("Error in mpfr_lngamma (768)\n");
202 exit (1);
203 }
204 if (inex >= 0)
205 {
206 printf ("Wrong flag for mpfr_lngamma (768)\n");
207 exit (1);
208 }
209
210 mpfr_set_prec (x, 4);
211 mpfr_set_prec (y, 4);
212 mpfr_set_str_binary (x, "0.1100E-66");
213 mpfr_lngamma (y, x, MPFR_RNDN);
214 mpfr_set_str_binary (x, "0.1100E6");
215 if (mpfr_cmp0 (x, y))
216 {
217 printf ("Error for lngamma(0.1100E-66)\n");
218 exit (1);
219 }
220
221 mpfr_set_prec (x, 256);
222 mpfr_set_prec (y, 32);
223 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN);
224 mpfr_add_ui (x, x, 1, MPFR_RNDN);
225 mpfr_div_2ui (x, x, 1, MPFR_RNDN);
226 mpfr_lngamma (y, x, MPFR_RNDN);
227 mpfr_set_prec (x, 32);
228 mpfr_set_str_binary (x, "-0.10001000111011111011000010100010E207");
229 if (mpfr_cmp0 (x, y))
230 {
231 printf ("Error for lngamma(-2^199+0.5)\n");
232 printf ("Got ");
233 mpfr_dump (y);
234 printf ("instead of ");
235 mpfr_dump (x);
236 exit (1);
237 }
238
239 mpfr_set_prec (x, 256);
240 mpfr_set_prec (y, 32);
241 mpfr_set_si_2exp (x, -1, 200, MPFR_RNDN);
242 mpfr_sub_ui (x, x, 1, MPFR_RNDN);
243 mpfr_div_2ui (x, x, 1, MPFR_RNDN);
244 mpfr_lngamma (y, x, MPFR_RNDN);
245 if (!mpfr_nan_p (y))
246 {
247 printf ("Error for lngamma(-2^199-0.5)\n");
248 exit (1);
249 }
250
251 mpfr_clear (x);
252 mpfr_clear (y);
253 }
254
255 /* test failing with GMP_CHECK_RANDOMIZE=1513869588 */
256 static void
257 bug20171220 (void)
258 {
259 mpfr_t x, y, z;
260 int inex;
261
262 mpfr_init2 (x, 15);
263 mpfr_init2 (y, 15);
264 mpfr_init2 (z, 15);
265
266 mpfr_set_str (x, "1.01111e+00", 10, MPFR_RNDN); /* x = 8283/8192 */
267 inex = mpfr_lngamma (y, x, MPFR_RNDN);
268 mpfr_set_str (z, "-0.0063109971733698154140545190234", 10, MPFR_RNDN);
269 MPFR_ASSERTN(mpfr_equal_p (y, z));
270 MPFR_ASSERTN(inex > 0);
271
272 mpfr_set_prec (x, 43);
273 mpfr_set_prec (y, 1);
274 mpfr_set_prec (z, 1);
275 mpfr_set_str (x, "9.8888652212918e-01", 10, MPFR_RNDN);
276 /* lngamma(x) = 0.00651701007222520, should be rounded up to 0.0078125 */
277 inex = mpfr_lngamma (y, x, MPFR_RNDU);
278 mpfr_set_ui_2exp (z, 1, -7, MPFR_RNDN);
279 MPFR_ASSERTN(mpfr_equal_p (y, z));
280 MPFR_ASSERTN(inex > 0);
281
282 mpfr_clear (x);
283 mpfr_clear (y);
284 mpfr_clear (z);
285 }
286
287 /* taway failing with GMP_CHECK_RANDOMIZE=1513888822 */
288 static void
289 bug20171220a (void)
290 {
291 mpfr_t x, y, z;
292 int inex;
293
294 mpfr_init2 (x, 198);
295 mpfr_init2 (y, 1);
296 mpfr_init2 (z, 1);
297
298 mpfr_set_str (x, "9.999962351340362288585900348170984233205352566408878552154832e-01", 10, MPFR_RNDN);
299 inex = mpfr_lngamma (y, x, MPFR_RNDA);
300 /* lngamma(x) ~ 2.1731512683e0-6 ~ 2^-18.81, should be rounded to 2^-18 */
301 mpfr_set_si_2exp (z, 1, -18, MPFR_RNDN);
302 MPFR_ASSERTN(mpfr_equal_p (y, z));
303 MPFR_ASSERTN(inex > 0);
304
305 mpfr_clear (x);
306 mpfr_clear (y);
307 mpfr_clear (z);
308 }
309
310 int
311 main (void)
312 {
313 tests_start_mpfr ();
314
315 bug20171220 ();
316 bug20171220a ();
317
318 special ();
319 test_generic (MPFR_PREC_MIN, 100, 2);
320
321 tests_end_mpfr ();
322 return 0;
323 }