1 /* Test file for mpfr_root (also used for mpfr_rootn_ui).
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 /* Note: troot.c is included by trootn_ui.c with TF = mpfr_rootn_ui */
24 #ifdef TF
25 # define TF_IS_MPFR_ROOT 0
26 #else
27 # define TF_IS_MPFR_ROOT 1
28 # define TF mpfr_root
29 # define _MPFR_NO_DEPRECATED_ROOT
30 #endif
31
32 #include "mpfr-test.h"
33
34 #include <time.h>
35
36 /* return the cpu time in seconds */
37 static double
38 cputime (void)
39 {
40 return (double) clock () / (double) CLOCKS_PER_SEC;
41 }
42
43 #define DEFN(N) \
44 static int root##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
45 { return TF (y, x, N, rnd); } \
46 static int pow##N (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t rnd) \
47 { return mpfr_pow_ui (y, x, N, rnd); }
48
49 DEFN(2)
50 DEFN(3)
51 DEFN(4)
52 DEFN(5)
53 DEFN(17)
54 DEFN(120)
55
56 static void
57 special (void)
58 {
59 mpfr_t x, y;
60 int i;
61
62 mpfr_init (x);
63 mpfr_init (y);
64
65 /* root(NaN) = NaN */
66 mpfr_set_nan (x);
67 i = TF (y, x, 17, MPFR_RNDN);
68 if (!mpfr_nan_p (y))
69 {
70 printf ("Error: root(NaN,17) <> NaN\n");
71 exit (1);
72 }
73 MPFR_ASSERTN (i == 0);
74
75 /* root(+Inf) = +Inf */
76 mpfr_set_inf (x, 1);
77 i = TF (y, x, 42, MPFR_RNDN);
78 if (!mpfr_inf_p (y) || mpfr_sgn (y) < 0)
79 {
80 printf ("Error: root(+Inf,42) <> +Inf\n");
81 exit (1);
82 }
83 MPFR_ASSERTN (i == 0);
84
85 /* root(-Inf, 17) = -Inf */
86 mpfr_set_inf (x, -1);
87 i = TF (y, x, 17, MPFR_RNDN);
88 if (!mpfr_inf_p (y) || mpfr_sgn (y) > 0)
89 {
90 printf ("Error: root(-Inf,17) <> -Inf\n");
91 exit (1);
92 }
93 MPFR_ASSERTN (i == 0);
94
95 /* root(-Inf, 42) = NaN */
96 mpfr_set_inf (x, -1);
97 i = TF (y, x, 42, MPFR_RNDN);
98 if (!mpfr_nan_p (y))
99 {
100 printf ("Error: root(-Inf,42) <> -Inf\n");
101 exit (1);
102 }
103 MPFR_ASSERTN (i == 0);
104
105 /* root(+/-0, k) = +/-0, with the sign depending on TF.
106 * Before calling the function, we set y to NaN with the wrong sign,
107 * so that if the code of the function forgets to do something, this
108 * will be detected.
109 */
110 mpfr_set_zero (x, 1); /* x is +0 */
111 MPFR_SET_NAN (y);
112 MPFR_SET_NEG (y);
113 i = TF (y, x, 17, MPFR_RNDN);
114 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
115 {
116 printf ("Error: root(+0,17) <> +0\n");
117 exit (1);
118 }
119 MPFR_ASSERTN (i == 0);
120 MPFR_SET_NAN (y);
121 MPFR_SET_NEG (y);
122 i = TF (y, x, 42, MPFR_RNDN);
123 if (MPFR_NOTZERO (y) || MPFR_IS_NEG (y))
124 {
125 printf ("Error: root(+0,42) <> +0\n");
126 exit (1);
127 }
128 MPFR_ASSERTN (i == 0);
129 mpfr_set_zero (x, -1); /* x is -0 */
130 MPFR_SET_NAN (y);
131 MPFR_SET_POS (y);
132 i = TF (y, x, 17, MPFR_RNDN);
133 if (MPFR_NOTZERO (y) || MPFR_IS_POS (y))
134 {
135 printf ("Error: root(-0,17) <> -0\n");
136 exit (1);
137 }
138 MPFR_ASSERTN (i == 0);
139 MPFR_SET_NAN (y);
140 if (TF_IS_MPFR_ROOT)
141 MPFR_SET_POS (y);
142 else
143 MPFR_SET_NEG (y);
144 i = TF (y, x, 42, MPFR_RNDN);
145 if (MPFR_NOTZERO (y) ||
146 (TF_IS_MPFR_ROOT ? MPFR_IS_POS (y) : MPFR_IS_NEG (y)))
147 {
148 printf ("Error: root(-0,42) <> %c0\n",
149 TF_IS_MPFR_ROOT ? '-' : '+');
150 exit (1);
151 }
152 MPFR_ASSERTN (i == 0);
153
154 mpfr_set_prec (x, 53);
155 mpfr_set_str (x, "8.39005285514734966412e-01", 10, MPFR_RNDN);
156 TF (x, x, 3, MPFR_RNDN);
157 if (mpfr_cmp_str1 (x, "9.43166207799662426048e-01"))
158 {
159 printf ("Error in root3 (1)\n");
160 printf ("expected 9.43166207799662426048e-01\n");
161 printf ("got ");
162 mpfr_dump (x);
163 exit (1);
164 }
165
166 mpfr_set_prec (x, 32);
167 mpfr_set_prec (y, 32);
168 mpfr_set_str_binary (x, "0.10000100001100101001001001011001");
169 TF (x, x, 3, MPFR_RNDN);
170 mpfr_set_str_binary (y, "0.11001101011000100111000111111001");
171 if (mpfr_cmp (x, y))
172 {
173 printf ("Error in root3 (2)\n");
174 exit (1);
175 }
176
177 mpfr_set_prec (x, 32);
178 mpfr_set_prec (y, 32);
179 mpfr_set_str_binary (x, "-0.1100001110110000010101011001011");
180 TF (x, x, 3, MPFR_RNDD);
181 mpfr_set_str_binary (y, "-0.11101010000100100101000101011001");
182 if (mpfr_cmp (x, y))
183 {
184 printf ("Error in root3 (3)\n");
185 exit (1);
186 }
187
188 mpfr_set_prec (x, 82);
189 mpfr_set_prec (y, 27);
190 mpfr_set_str_binary (x, "0.1010001111011101011011000111001011001101100011110110010011011011011010011001100101e-7");
191 TF (y, x, 3, MPFR_RNDD);
192 mpfr_set_str_binary (x, "0.101011110001110001000100011E-2");
193 if (mpfr_cmp (x, y))
194 {
195 printf ("Error in root3 (4)\n");
196 exit (1);
197 }
198
199 mpfr_set_prec (x, 204);
200 mpfr_set_prec (y, 38);
201 mpfr_set_str_binary (x, "0.101000000001101000000001100111111011111001110110100001111000100110100111001101100111110001110001011011010110010011100101111001111100001010010100111011101100000011011000101100010000000011000101001010001001E-5");
202 TF (y, x, 3, MPFR_RNDD);
203 mpfr_set_str_binary (x, "0.10001001111010011011101000010110110010E-1");
204 if (mpfr_cmp (x, y))
205 {
206 printf ("Error in root3 (5)\n");
207 exit (1);
208 }
209
210 /* Worst case found on 2006-11-25 */
211 mpfr_set_prec (x, 53);
212 mpfr_set_prec (y, 53);
213 mpfr_set_str_binary (x, "1.0100001101101101001100110001001000000101001101100011E28");
214 TF (y, x, 35, MPFR_RNDN);
215 mpfr_set_str_binary (x, "1.1100000010110101100011101011000010100001101100100011E0");
216 if (mpfr_cmp (x, y))
217 {
218 printf ("Error in rootn (y, x, 35, MPFR_RNDN) for\n"
219 "x = 1.0100001101101101001100110001001000000101001101100011E28\n"
220 "Expected ");
221 mpfr_dump (x);
222 printf ("Got ");
223 mpfr_dump (y);
224 exit (1);
225 }
226 /* Worst cases found on 2006-11-26 */
227 mpfr_set_str_binary (x, "1.1111010011101110001111010110000101110000110110101100E17");
228 TF (y, x, 36, MPFR_RNDD);
229 mpfr_set_str_binary (x, "1.0110100111010001101001010111001110010100111111000010E0");
230 if (mpfr_cmp (x, y))
231 {
232 printf ("Error in rootn (y, x, 36, MPFR_RNDD) for\n"
233 "x = 1.1111010011101110001111010110000101110000110110101100E17\n"
234 "Expected ");
235 mpfr_dump (x);
236 printf ("Got ");
237 mpfr_dump (y);
238 exit (1);
239 }
240 mpfr_set_str_binary (x, "1.1100011101101101100010110001000001110001111110010000E23");
241 TF (y, x, 36, MPFR_RNDU);
242 mpfr_set_str_binary (x, "1.1001010100001110000110111111100011011101110011000100E0");
243 if (mpfr_cmp (x, y))
244 {
245 printf ("Error in rootn (y, x, 36, MPFR_RNDU) for\n"
246 "x = 1.1100011101101101100010110001000001110001111110010000E23\n"
247 "Expected ");
248 mpfr_dump (x);
249 printf ("Got ");
250 mpfr_dump (y);
251 exit (1);
252 }
253
254 /* Check for k = 1 */
255 mpfr_set_ui (x, 17, MPFR_RNDN);
256 i = TF (y, x, 1, MPFR_RNDN);
257 if (mpfr_cmp_ui (x, 17) || i != 0)
258 {
259 printf ("Error in root for 17^(1/1)\n");
260 exit (1);
261 }
262
263 mpfr_set_ui (x, 0, MPFR_RNDN);
264 i = TF (y, x, 0, MPFR_RNDN);
265 if (!MPFR_IS_NAN (y) || i != 0)
266 {
267 printf ("Error in root for (+0)^(1/0)\n");
268 exit (1);
269 }
270 mpfr_neg (x, x, MPFR_RNDN);
271 i = TF (y, x, 0, MPFR_RNDN);
272 if (!MPFR_IS_NAN (y) || i != 0)
273 {
274 printf ("Error in root for (-0)^(1/0)\n");
275 exit (1);
276 }
277
278 mpfr_set_ui (x, 1, MPFR_RNDN);
279 i = TF (y, x, 0, MPFR_RNDN);
280 if (!MPFR_IS_NAN (y) || i != 0)
281 {
282 printf ("Error in root for 1^(1/0)\n");
283 exit (1);
284 }
285
286 /* Check for k==2 */
287 mpfr_set_si (x, -17, MPFR_RNDD);
288 i = TF (y, x, 2, MPFR_RNDN);
289 if (!MPFR_IS_NAN (y) || i != 0)
290 {
291 printf ("Error in root for (-17)^(1/2)\n");
292 exit (1);
293 }
294
295 mpfr_clear (x);
296 mpfr_clear (y);
297 }
298
299 /* https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=812779
300 * https://bugzilla.gnome.org/show_bug.cgi?id=756960
301 * is a GNOME Calculator bug (mpfr_root applied on a negative integer,
302 * which is converted to an unsigned integer), but the strange result
303 * is also due to a bug in MPFR.
304 */
305 static void
306 bigint (void)
307 {
308 mpfr_t x, y;
309
310 mpfr_inits2 (64, x, y, (mpfr_ptr) 0);
311
312 mpfr_set_ui (x, 10, MPFR_RNDN);
313 if (sizeof (unsigned long) * CHAR_BIT == 64)
314 {
315 TF (x, x, ULONG_MAX, MPFR_RNDN);
316 mpfr_set_ui_2exp (y, 1, -63, MPFR_RNDN);
317 mpfr_add_ui (y, y, 1, MPFR_RNDN);
318 if (! mpfr_equal_p (x, y))
319 {
320 printf ("Error in bigint for ULONG_MAX\n");
321 printf ("Expected ");
322 mpfr_dump (y);
323 printf ("Got ");
324 mpfr_dump (x);
325 exit (1);
326 }
327 }
328
329 mpfr_set_ui (x, 10, MPFR_RNDN);
330 TF (x, x, 1234567890, MPFR_RNDN);
331 mpfr_set_str_binary (y,
332 "1.00000000000000000000000000001000000000101011000101000110010001");
333 if (! mpfr_equal_p (x, y))
334 {
335 printf ("Error in bigint for 1234567890\n");
336 printf ("Expected ");
337 mpfr_dump (y);
338 printf ("Got ");
339 mpfr_dump (x);
340 exit (1);
341 }
342
343 mpfr_clears (x, y, (mpfr_ptr) 0);
344 }
345
346 #define TEST_FUNCTION TF
347 #define INTEGER_TYPE unsigned long
348 #define INT_RAND_FUNCTION() \
349 (INTEGER_TYPE) (RAND_BOOL () ? randlimb () : randlimb () % 3 + 2)
350 #include "tgeneric_ui.c"
351
352 static void
353 exact_powers (unsigned long bmax, unsigned long kmax)
354 {
355 long b, k;
356 mpz_t z;
357 mpfr_t x, y;
358 int inex, neg;
359
360 mpz_init (z);
361 for (b = 2; b <= bmax; b++)
362 for (k = 1; k <= kmax; k++)
363 {
364 mpz_ui_pow_ui (z, b, k);
365 mpfr_init2 (x, mpz_sizeinbase (z, 2));
366 mpfr_set_ui (x, b, MPFR_RNDN);
367 mpfr_pow_ui (x, x, k, MPFR_RNDN);
368 mpz_set_ui (z, b);
369 mpfr_init2 (y, mpz_sizeinbase (z, 2));
370 for (neg = 0; neg <= 1; neg++)
371 {
372 inex = TF (y, x, k, MPFR_RNDN);
373 if (inex != 0)
374 {
375 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
376 printf ("Expected inex=0, got %d\n", inex);
377 exit (1);
378 }
379 if (neg && (k & 1) == 0)
380 {
381 if (!MPFR_IS_NAN (y))
382 {
383 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
384 printf ("Expected y=NaN\n");
385 printf ("Got ");
386 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
387 printf ("\n");
388 exit (1);
389 }
390 }
391 else if (MPFR_IS_NAN (y) || mpfr_cmp_si (y, b) != 0)
392 {
393 printf ("Error in exact_powers, b=%ld, k=%ld\n", b, k);
394 printf ("Expected y=%ld\n", b);
395 printf ("Got ");
396 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDN);
397 printf ("\n");
398 exit (1);
399 }
400 mpfr_neg (x, x, MPFR_RNDN);
401 b = -b;
402 }
403 mpfr_clear (x);
404 mpfr_clear (y);
405 }
406 mpz_clear (z);
407 }
408
409 /* Compare root(x,2^h) with pow(x,2^(-h)). */
410 static void
411 cmp_pow (void)
412 {
413 mpfr_t x, y1, y2;
414 int h;
415
416 mpfr_inits2 (128, x, y1, y2, (mpfr_ptr) 0);
417
418 for (h = 1; h < sizeof (unsigned long) * CHAR_BIT; h++)
419 {
420 unsigned long k = (unsigned long) 1 << h;
421 int i;
422
423 for (i = 0; i < 10; i++)
424 {
425 mpfr_rnd_t rnd;
426 mpfr_flags_t flags1, flags2;
427 int inex1, inex2;
428
429 tests_default_random (x, 0, __gmpfr_emin, __gmpfr_emax, 1);
430 rnd = RND_RAND_NO_RNDF ();
431 mpfr_set_ui_2exp (y1, 1, -h, MPFR_RNDN);
432 mpfr_clear_flags ();
433 inex1 = mpfr_pow (y1, x, y1, rnd);
434 flags1 = __gmpfr_flags;
435 mpfr_clear_flags ();
436 inex2 = TF (y2, x, k, rnd);
437 flags2 = __gmpfr_flags;
438 if (!(mpfr_equal_p (y1, y2) && SAME_SIGN (inex1, inex2) &&
439 flags1 == flags2))
440 {
441 printf ("Error in cmp_pow on h=%d, i=%d, rnd=%s\n",
442 h, i, mpfr_print_rnd_mode ((mpfr_rnd_t) rnd));
443 printf ("x = ");
444 mpfr_dump (x);
445 printf ("pow = ");
446 mpfr_dump (y1);
447 printf ("with inex = %d, flags =", inex1);
448 flags_out (flags1);
449 printf ("root = ");
450 mpfr_dump (y2);
451 printf ("with inex = %d, flags =", inex2);
452 flags_out (flags2);
453 exit (1);
454 }
455 }
456 }
457
458 mpfr_clears (x, y1, y2, (mpfr_ptr) 0);
459 }
460
461 static void
462 bug20171214 (void)
463 {
464 mpfr_t x, y;
465 int inex;
466
467 mpfr_init2 (x, 805);
468 mpfr_init2 (y, 837);
469 mpfr_set_ui (x, 1, MPFR_RNDN);
470 inex = TF (y, x, 120, MPFR_RNDN);
471 MPFR_ASSERTN (inex == 0);
472 MPFR_ASSERTN (mpfr_cmp_ui (y, 1) == 0);
473 mpfr_set_si (x, -1, MPFR_RNDN);
474 inex = TF (y, x, 121, MPFR_RNDN);
475 MPFR_ASSERTN (inex == 0);
476 MPFR_ASSERTN (mpfr_cmp_si (y, -1) == 0);
477 mpfr_clear (x);
478 mpfr_clear (y);
479 }
480
481 int
482 main (int argc, char *argv[])
483 {
484 mpfr_t x;
485 int r;
486 mpfr_prec_t p;
487 unsigned long k;
488
489 if (argc == 3) /* troot prec k */
490 {
491 double st1, st2;
492 unsigned long k;
493 int l;
494 mpfr_t y;
495 p = strtoul (argv[1], NULL, 10);
496 k = strtoul (argv[2], NULL, 10);
497 mpfr_init2 (x, p);
498 mpfr_init2 (y, p);
499 mpfr_const_pi (y, MPFR_RNDN);
500 TF (x, y, k, MPFR_RNDN); /* to warm up cache */
501 st1 = cputime ();
502 for (l = 0; cputime () - st1 < 1.0; l++)
503 TF (x, y, k, MPFR_RNDN);
504 st1 = (cputime () - st1) / l;
505 printf ("%-15s took %.2es\n", MAKE_STR(TF), st1);
506
507 /* compare with x^(1/k) = exp(1/k*log(x)) */
508 /* first warm up cache */
509 mpfr_swap (x, y);
510 mpfr_log (y, x, MPFR_RNDN);
511 mpfr_div_ui (y, y, k, MPFR_RNDN);
512 mpfr_exp (y, y, MPFR_RNDN);
513
514 st2 = cputime ();
515 for (l = 0; cputime () - st2 < 1.0; l++)
516 {
517 mpfr_log (y, x, MPFR_RNDN);
518 mpfr_div_ui (y, y, k, MPFR_RNDN);
519 mpfr_exp (y, y, MPFR_RNDN);
520 }
521 st2 = (cputime () - st2) / l;
522 printf ("exp(1/k*log(x)) took %.2es\n", st2);
523
524 mpfr_clear (x);
525 mpfr_clear (y);
526 return 0;
527 }
528
529 tests_start_mpfr ();
530
531 bug20171214 ();
532 exact_powers (3, 1000);
533 special ();
534 bigint ();
535 cmp_pow ();
536
537 mpfr_init (x);
538
539 for (p = MPFR_PREC_MIN; p < 100; p++)
540 {
541 mpfr_set_prec (x, p);
542 RND_LOOP (r)
543 {
544 mpfr_set_ui (x, 1, MPFR_RNDN);
545 k = 2 + randlimb () % 4; /* 2 <= k <= 5 */
546 TF (x, x, k, (mpfr_rnd_t) r);
547 if (mpfr_cmp_ui (x, 1))
548 {
549 printf ("Error in rootn[%lu] for x=1, rnd=%s\ngot ",
550 k, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
551 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
552 printf ("\n");
553 exit (1);
554 }
555 mpfr_set_si (x, -1, MPFR_RNDN);
556 if (k % 2)
557 {
558 TF (x, x, k, (mpfr_rnd_t) r);
559 if (mpfr_cmp_si (x, -1))
560 {
561 printf ("Error in rootn[%lu] for x=-1, rnd=%s\ngot ",
562 k, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
563 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
564 printf ("\n");
565 exit (1);
566 }
567 }
568
569 if (p >= 5)
570 {
571 int i;
572 for (i = -12; i <= 12; i++)
573 {
574 mpfr_set_ui (x, 27, MPFR_RNDN);
575 mpfr_mul_2si (x, x, 3*i, MPFR_RNDN);
576 TF (x, x, 3, MPFR_RNDN);
577 if (mpfr_cmp_si_2exp (x, 3, i))
578 {
579 printf ("Error in rootn[3] for "
580 "x = 27.0 * 2^(%d), rnd=%s\ngot ",
581 3*i, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
582 mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
583 printf ("\ninstead of 3 * 2^(%d)\n", i);
584 exit (1);
585 }
586 }
587 }
588 }
589 }
590 mpfr_clear (x);
591
592 test_generic_ui (MPFR_PREC_MIN, 200, 30);
593
594 /* The sign of the random value y (used to generate a potential bad case)
595 is negative with a probability 256/512 = 1/2 for odd n, and never
596 negative (probability 0/512) for even n (if y is negative, then
597 (y^(2k))^(1/(2k)) is different from y, so that this would yield
598 an error). */
599 bad_cases (root2, pow2, "rootn[2]", 0, -256, 255, 4, 128, 80, 40);
600 bad_cases (root3, pow3, "rootn[3]", 256, -256, 255, 4, 128, 200, 40);
601 bad_cases (root4, pow4, "rootn[4]", 0, -256, 255, 4, 128, 320, 40);
602 bad_cases (root5, pow5, "rootn[5]", 256, -256, 255, 4, 128, 440, 40);
603 bad_cases (root17, pow17, "rootn[17]", 256, -256, 255, 4, 128, 800, 40);
604 bad_cases (root120, pow120, "rootn[120]", 0, -256, 255, 4, 128, 800, 40);
605
606 tests_end_mpfr ();
607 return 0;
608 }