1 /* Test file for mpfr_compound_si.
2
3 Copyright 2021-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_compound_si
26 #define INTEGER_TYPE long
27 #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS)
28 #define test_generic_ui test_generic_si
29 #include "tgeneric_ui.c"
30
31 /* Special cases from IEEE 754-2019 */
32 static void
33 check_ieee754 (void)
34 {
35 mpfr_t x, y;
36 long i;
37 long t[] = { 0, 1, 2, 3, 17, LONG_MAX-1, LONG_MAX };
38 int j;
39 mpfr_prec_t prec = 2; /* we need at least 2 so that 3/4 is exact */
40
41 mpfr_init2 (x, prec);
42 mpfr_init2 (y, prec);
43
44 /* compound(x,n) = NaN for x < -1, and set invalid exception */
45 for (i = 0; i < numberof(t); i++)
46 for (j = 0; j < 2; j++)
47 {
48 const char *s;
49
50 mpfr_clear_nanflag ();
51 if (j == 0)
52 {
53 mpfr_set_si (x, -2, MPFR_RNDN);
54 s = "-2";
55 }
56 else
57 {
58 mpfr_set_inf (x, -1);
59 s = "-Inf";
60 }
61 mpfr_compound_si (y, x, t[i], MPFR_RNDN);
62 if (!mpfr_nan_p (y))
63 {
64 printf ("Error, compound(%s,%ld) should give NaN\n", s, t[i]);
65 printf ("got "); mpfr_dump (y);
66 exit (1);
67 }
68 if (!mpfr_nanflag_p ())
69 {
70 printf ("Error, compound(%s,%ld) should raise invalid flag\n",
71 s, t[i]);
72 exit (1);
73 }
74 }
75
76 /* compound(x,0) = 1 for x >= -1 or x = NaN */
77 for (i = -2; i <= 2; i++)
78 {
79 if (i == -2)
80 mpfr_set_nan (x);
81 else if (i == 2)
82 mpfr_set_inf (x, 1);
83 else
84 mpfr_set_si (x, i, MPFR_RNDN);
85 mpfr_compound_si (y, x, 0, MPFR_RNDN);
86 if (mpfr_cmp_ui (y, 1) != 0)
87 {
88 printf ("Error, compound(x,0) should give 1 on\nx = ");
89 mpfr_dump (x);
90 printf ("got "); mpfr_dump (y);
91 exit (1);
92 }
93 }
94
95 /* compound(-1,n) = +Inf for n < 0, and raise divide-by-zero flag */
96 mpfr_clear_divby0 ();
97 mpfr_set_si (x, -1, MPFR_RNDN);
98 mpfr_compound_si (y, x, -1, MPFR_RNDN);
99 if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0)
100 {
101 printf ("Error, compound(-1,-1) should give +Inf\n");
102 printf ("got "); mpfr_dump (y);
103 exit (1);
104 }
105 if (!mpfr_divby0_p ())
106 {
107 printf ("Error, compound(-1,-1) should raise divide-by-zero flag\n");
108 exit (1);
109 }
110
111 /* compound(-1,n) = +0 for n > 0 */
112 mpfr_set_si (x, -1, MPFR_RNDN);
113 mpfr_compound_si (y, x, 1, MPFR_RNDN);
114 if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0)
115 {
116 printf ("Error, compound(-1,1) should give +0\n");
117 printf ("got "); mpfr_dump (y);
118 exit (1);
119 }
120
121 /* compound(+/-0,n) = 1 */
122 for (i = -1; i <= 1; i++)
123 {
124 mpfr_set_zero (x, -1);
125 mpfr_compound_si (y, x, i, MPFR_RNDN);
126 if (mpfr_cmp_ui (y, 1) != 0)
127 {
128 printf ("Error1, compound(x,%ld) should give 1\non x = ", i);
129 mpfr_dump (x);
130 printf ("got "); mpfr_dump (y);
131 exit (1);
132 }
133 mpfr_set_zero (x, +1);
134 mpfr_compound_si (y, x, i, MPFR_RNDN);
135 if (mpfr_cmp_ui (y, 1) != 0)
136 {
137 printf ("Error, compound(x,%ld) should give 1\non x = ", i);
138 mpfr_dump (x);
139 printf ("got "); mpfr_dump (y);
140 exit (1);
141 }
142 }
143
144 /* compound(+Inf,n) = +Inf for n > 0 */
145 mpfr_set_inf (x, 1);
146 mpfr_compound_si (y, x, 1, MPFR_RNDN);
147 if (!mpfr_inf_p (y) || MPFR_SIGN(y) < 0)
148 {
149 printf ("Error, compound(+Inf,1) should give +Inf\n");
150 printf ("got "); mpfr_dump (y);
151 exit (1);
152 }
153
154 /* compound(+Inf,n) = +0 for n < 0 */
155 mpfr_set_inf (x, 1);
156 mpfr_compound_si (y, x, -1, MPFR_RNDN);
157 if (!mpfr_zero_p (y) || MPFR_SIGN(y) < 0)
158 {
159 printf ("Error, compound(+Inf,-1) should give +0\n");
160 printf ("got "); mpfr_dump (y);
161 exit (1);
162 }
163
164 /* compound(NaN,n) = NaN for n <> 0 */
165 mpfr_set_nan (x);
166 mpfr_compound_si (y, x, -1, MPFR_RNDN);
167 if (!mpfr_nan_p (y))
168 {
169 printf ("Error, compound(NaN,-1) should give NaN\n");
170 printf ("got "); mpfr_dump (y);
171 exit (1);
172 }
173 mpfr_compound_si (y, x, +1, MPFR_RNDN);
174 if (!mpfr_nan_p (y))
175 {
176 printf ("Error, compound(NaN,+1) should give NaN\n");
177 printf ("got "); mpfr_dump (y);
178 exit (1);
179 }
180
181 /* hard-coded test: x is the 32-bit nearest approximation of 17/42 */
182 mpfr_set_prec (x, 32);
183 mpfr_set_prec (y, 32);
184 mpfr_set_ui_2exp (x, 3476878287UL, -33, MPFR_RNDN);
185 mpfr_compound_si (y, x, 12, MPFR_RNDN);
186 mpfr_set_ui_2exp (x, 1981447393UL, -25, MPFR_RNDN);
187 if (!mpfr_equal_p (y, x))
188 {
189 printf ("Error for compound(3476878287/2^33,12)\n");
190 printf ("expected "); mpfr_dump (x);
191 printf ("got "); mpfr_dump (y);
192 exit (1);
193 }
194
195 /* test for negative n */
196 i = -1;
197 while (1)
198 {
199 /* i has the form -(2^k-1) */
200 mpfr_set_si_2exp (x, -1, -1, MPFR_RNDN); /* x = -0.5 */
201 mpfr_compound_si (y, x, i, MPFR_RNDN);
202 mpfr_set_ui_2exp (x, 1, -i, MPFR_RNDN);
203 if (!mpfr_equal_p (y, x))
204 {
205 printf ("Error for compound(-0.5,%ld)\n", i);
206 printf ("expected "); mpfr_dump (x);
207 printf ("got "); mpfr_dump (y);
208 exit (1);
209 }
210 if (i == -2147483647) /* largest possible value on 32-bit machine */
211 break;
212 i = 2 * i - 1;
213 }
214
215 /* The "#if" makes sure that 64-bit constants are supported, avoiding
216 a compilation failure. The "if" makes sure that the constant is
217 representable in a long (this would not be the case with 32-bit
218 unsigned long and 64-bit limb). */
219 #if GMP_NUMB_BITS >= 64 || MPFR_PREC_BITS >= 64
220 if (4994322635099777669 <= LONG_MAX)
221 {
222 i = -4994322635099777669;
223 mpfr_set_ui (x, 1, MPFR_RNDN);
224 mpfr_compound_si (y, x, i, MPFR_RNDN);
225 mpfr_set_si (x, 1, MPFR_RNDN);
226 mpfr_mul_2si (x, x, i, MPFR_RNDN);
227 if (!mpfr_equal_p (y, x))
228 {
229 printf ("Error for compound(1,%ld)\n", i);
230 printf ("expected "); mpfr_dump (x);
231 printf ("got "); mpfr_dump (y);
232 exit (1);
233 }
234 }
235 #endif
236
237 mpfr_clear (x);
238 mpfr_clear (y);
239 }
240
241 /* Failure with mpfr_compound_si from 2021-02-15 due to
242 incorrect underflow detection. */
243 static void
244 bug_20230206 (void)
245 {
246 if (MPFR_PREC_MIN == 1)
247 {
248 mpfr_t x, y1, y2;
249 int inex1, inex2;
250 mpfr_flags_t flags1, flags2;
251 #if MPFR_PREC_BITS >= 64
252 mpfr_exp_t emin;
253 #endif
254
255 mpfr_inits2 (1, x, y1, y2, (mpfr_ptr) 0);
256 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); /* x = 1/2 */
257
258 /* This first test is useful mainly for a 32-bit mpfr_exp_t type
259 (no failure with a 64-bit mpfr_exp_t type since the underflow
260 threshold in the extended exponent range is much lower). */
261
262 mpfr_set_ui_2exp (y1, 1, -1072124363, MPFR_RNDN);
263 inex1 = -1;
264 flags1 = MPFR_FLAGS_INEXACT;
265 mpfr_clear_flags ();
266 /* -1832808704 ~= -2^30 / log2(3/2) */
267 inex2 = mpfr_compound_si (y2, x, -1832808704, MPFR_RNDN);
268 flags2 = __gmpfr_flags;
269 if (!(mpfr_equal_p (y1, y2) &&
270 SAME_SIGN (inex1, inex2) &&
271 flags1 == flags2))
272 {
273 printf ("Error in bug_20230206 (1):\n");
274 printf ("Expected ");
275 mpfr_dump (y1);
276 printf (" with inex = %d, flags =", inex1);
277 flags_out (flags1);
278 printf ("Got ");
279 mpfr_dump (y2);
280 printf (" with inex = %d, flags =", inex2);
281 flags_out (flags2);
282 exit (1);
283 }
284
285 /* This second test is for a 64-bit mpfr_exp_t type
286 (it is disabled with a 32-bit mpfr_exp_t type). */
287
288 /* The "#if" makes sure that 64-bit constants are supported, avoiding
289 a compilation failure. The "if" makes sure that the constant is
290 representable in a long (this would not be the case with 32-bit
291 unsigned long and 64-bit limb). It also ensures that mpfr_exp_t
292 has at least 64 bits. */
293 #if MPFR_PREC_BITS >= 64
294 emin = mpfr_get_emin ();
295 set_emin (MPFR_EMIN_MIN);
296 mpfr_set_ui_2exp (y1, 1, -4611686018427366846, MPFR_RNDN);
297 inex1 = 1;
298 flags1 = MPFR_FLAGS_INEXACT;
299 mpfr_clear_flags ();
300 /* -7883729320669216768 ~= -2^62 / log2(3/2) */
301 inex2 = mpfr_compound_si (y2, x, -7883729320669216768, MPFR_RNDN);
302 flags2 = __gmpfr_flags;
303 if (!(mpfr_equal_p (y1, y2) &&
304 SAME_SIGN (inex1, inex2) &&
305 flags1 == flags2))
306 {
307 printf ("Error in bug_20230206 (2):\n");
308 printf ("Expected ");
309 mpfr_dump (y1);
310 printf (" with inex = %d, flags =", inex1);
311 flags_out (flags1);
312 printf ("Got ");
313 mpfr_dump (y2);
314 printf (" with inex = %d, flags =", inex2);
315 flags_out (flags2);
316 exit (1);
317 }
318 set_emin (emin);
319 #endif
320
321 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
322 }
323 }
324
325 /* Reported by Patrick Pelissier on 2023-02-11 for the master branch
326 (tgeneric_ui.c with GMP_CHECK_RANDOMIZE=1412991715).
327 On a 32-bit host, one gets Inf (overflow) instead of 0.1E1071805703.
328 */
329 static void
330 bug_20230211 (void)
331 {
332 mpfr_t x, y1, y2;
333 int inex1, inex2;
334 mpfr_flags_t flags1, flags2;
335
336 mpfr_inits2 (1, x, y1, y2, (mpfr_ptr) 0);
337 mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); /* x = 1/2 */
338 mpfr_set_ui_2exp (y1, 1, 1071805702, MPFR_RNDN);
339 inex1 = 1;
340 flags1 = MPFR_FLAGS_INEXACT;
341 mpfr_clear_flags ();
342 inex2 = mpfr_compound_si (y2, x, 1832263949, MPFR_RNDN);
343 flags2 = __gmpfr_flags;
344 if (!(mpfr_equal_p (y1, y2) &&
345 SAME_SIGN (inex1, inex2) &&
346 flags1 == flags2))
347 {
348 printf ("Error in bug_20230211:\n");
349 printf ("Expected ");
350 mpfr_dump (y1);
351 printf (" with inex = %d, flags =", inex1);
352 flags_out (flags1);
353 printf ("Got ");
354 mpfr_dump (y2);
355 printf (" with inex = %d, flags =", inex2);
356 flags_out (flags2);
357 exit (1);
358 }
359 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
360 }
361
362 /* Integer overflow with compound.c d04caeae04c6a83276916c4fbac1fe9b0cec3c8b
363 (2023-02-23) or 952fb0f5cc2df1fffde3eb54c462fdae5f123ea6 in the 4.2 branch
364 on "n * (kx - 1) + 1". Note: if the only effect is just a random value,
365 this probably doesn't affect the result (one might enter the "if" while
366 one shouldn't, but the real check is done inside the "if"). This test
367 fails if -fsanitize=undefined -fno-sanitize-recover is used or if the
368 processor emits a signal in case of integer overflow.
369 This test has been made obsolete by the "kx < ex" condition
370 in 2cb3123891dd46fe0258d4aec7f8655b8ec69aaf (master branch)
371 or f5cb40571bc3d1559f05b230cf4ffecaf0952852 (4.2 branch). */
372 static void
373 bug_20230517 (void)
374 {
375 mpfr_exp_t old_emax;
376 mpfr_t x;
377
378 old_emax = mpfr_get_emax ();
379 set_emax (MPFR_EMAX_MAX);
380
381 mpfr_init2 (x, 123456);
382 mpfr_set_ui (x, 65536, MPFR_RNDN);
383 mpfr_nextabove (x);
384 mpfr_compound_si (x, x, LONG_MAX >> 16, MPFR_RNDN);
385 mpfr_clear (x);
386
387 set_emax (old_emax);
388 }
389
390 /* Inverse function on non-special cases...
391 One has x = (1+y)^n with y > -1 and x > 0. Thus y = x^(1/n) - 1.
392 The inverse function is useful
393 - to build and check hard-to-round cases (see bad_cases() in tests.c);
394 - to test the behavior close to the overflow and underflow thresholds.
395 The case x = 0 actually needs to be handled as it may occur with
396 bad_cases() due to rounding.
397 */
398 static int
399 inv_compound (mpfr_ptr y, mpfr_srcptr x, long n, mpfr_rnd_t rnd_mode)
400 {
401 mpfr_t t;
402 int inexact;
403 mpfr_prec_t precy, prect;
404 MPFR_ZIV_DECL (loop);
405 MPFR_SAVE_EXPO_DECL (expo);
406
407 MPFR_ASSERTN (n != 0);
408
409 if (MPFR_UNLIKELY (MPFR_IS_ZERO (x)))
410 {
411 if (n > 0)
412 return mpfr_set_si (y, -1, rnd_mode);
413 else
414 {
415 MPFR_SET_INF (y);
416 MPFR_SET_POS (y);
417 MPFR_RET (0);
418 }
419 }
420
421 MPFR_SAVE_EXPO_MARK (expo);
422
423 if (mpfr_equal_p (x, __gmpfr_one))
424 {
425 MPFR_SAVE_EXPO_FREE (expo);
426 mpfr_set_zero (y, 1);
427 MPFR_RET (0);
428 }
429
430 precy = MPFR_GET_PREC (y);
431 prect = precy + 20;
432 mpfr_init2 (t, prect);
433
434 MPFR_ZIV_INIT (loop, prect);
435 for (;;)
436 {
437 mpfr_exp_t expt1, expt2, err;
438 unsigned int inext;
439
440 if (mpfr_rootn_si (t, x, n, MPFR_RNDN) == 0)
441 {
442 /* With a huge t, this case would yield inext != 0 and a
443 MPFR_CAN_ROUND failure until a huge precision is reached
444 (as the result is very close to an exact point). Fortunately,
445 since t is exact, we can obtain the correctly rounded result
446 by doing the second operation to the target precision directly.
447 */
448 inexact = mpfr_sub_ui (y, t, 1, rnd_mode);
449 goto end;
450 }
451 expt1 = MPFR_GET_EXP (t);
452 /* |error| <= 2^(expt1-prect-1) */
453 inext = mpfr_sub_ui (t, t, 1, MPFR_RNDN);
454 if (MPFR_UNLIKELY (MPFR_IS_ZERO (t)))
455 goto cont; /* cannot round yet */
456 expt2 = MPFR_GET_EXP (t);
457 err = 1;
458 if (expt2 < expt1)
459 err += expt1 - expt2;
460 /* |error(rootn)| <= 2^(err+expt2-prect-2)
461 and if mpfr_sub_ui is inexact:
462 |error| <= 2^(err+expt2-prect-2) + 2^(expt2-prect-1)
463 <= (2^(err-1) + 1) * 2^(expt2-prect-1)
464 <= 2^((err+1)+expt2-prect-2) */
465 if (inext)
466 err++;
467 /* |error| <= 2^(err+expt2-prect-2) */
468 if (MPFR_CAN_ROUND (t, prect + 2 - err, precy, rnd_mode))
469 break;
470
471 cont:
472 MPFR_ZIV_NEXT (loop, prect);
473 mpfr_set_prec (t, prect);
474 }
475
476 inexact = mpfr_set (y, t, rnd_mode);
477
478 end:
479 MPFR_ZIV_FREE (loop);
480 mpfr_clear (t);
481 MPFR_SAVE_EXPO_FREE (expo);
482 return mpfr_check_range (y, inexact, rnd_mode);
483 }
484
485 #define DEFN(N) \
486 static int mpfr_compound##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) \
487 { return mpfr_compound_si (y, x, N, r); } \
488 static int inv_compound##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r) \
489 { return inv_compound (y, x, N, r); }
490
491 DEFN(2)
492 DEFN(3)
493 DEFN(4)
494 DEFN(5)
495 DEFN(17)
496 DEFN(120)
497
498 #define TEST_FUNCTION mpfr_compound2
499 #define test_generic test_generic_compound2
500 #include "tgeneric.c"
501
502 #define TEST_FUNCTION mpfr_compound3
503 #define test_generic test_generic_compound3
504 #include "tgeneric.c"
505
506 #define TEST_FUNCTION mpfr_compound4
507 #define test_generic test_generic_compound4
508 #include "tgeneric.c"
509
510 #define TEST_FUNCTION mpfr_compound5
511 #define test_generic test_generic_compound5
512 #include "tgeneric.c"
513
514 #define TEST_FUNCTION mpfr_compound17
515 #define test_generic test_generic_compound17
516 #include "tgeneric.c"
517
518 #define TEST_FUNCTION mpfr_compound120
519 #define test_generic test_generic_compound120
520 #include "tgeneric.c"
521
522 int
523 main (void)
524 {
525 tests_start_mpfr ();
526
527 check_ieee754 ();
528 bug_20230206 ();
529 bug_20230211 ();
530 bug_20230517 ();
531
532 test_generic_si (MPFR_PREC_MIN, 100, 100);
533
534 test_generic_compound2 (MPFR_PREC_MIN, 100, 100);
535 test_generic_compound3 (MPFR_PREC_MIN, 100, 100);
536 test_generic_compound4 (MPFR_PREC_MIN, 100, 100);
537 test_generic_compound5 (MPFR_PREC_MIN, 100, 100);
538 test_generic_compound17 (MPFR_PREC_MIN, 100, 100);
539 test_generic_compound120 (MPFR_PREC_MIN, 100, 100);
540
541 /* Note: For small n, we need a psup high enough to avoid too many
542 "f exact while f^(-1) inexact" occurrences in bad_cases(). */
543 bad_cases (mpfr_compound2, inv_compound2, "mpfr_compound2",
544 0, -256, 255, 4, 128, 240, 40);
545 bad_cases (mpfr_compound3, inv_compound3, "mpfr_compound3",
546 0, -256, 255, 4, 128, 120, 40);
547 bad_cases (mpfr_compound4, inv_compound4, "mpfr_compound4",
548 0, -256, 255, 4, 128, 80, 40);
549 bad_cases (mpfr_compound5, inv_compound5, "mpfr_compound5",
550 0, -256, 255, 4, 128, 80, 40);
551 bad_cases (mpfr_compound17, inv_compound17, "mpfr_compound17",
552 0, -256, 255, 4, 128, 80, 40);
553 bad_cases (mpfr_compound120, inv_compound120, "mpfr_compound120",
554 0, -256, 255, 4, 128, 80, 40);
555
556 tests_end_mpfr ();
557 return 0;
558 }