1/* Test file for mpfr_ui_pow and mpfr_ui_pow_ui. 2 3Copyright 2001-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-test.h" 24 25static void 26test1 (void) 27{ 28 mpfr_t x, y, z, a; 29 int res1, res2; 30 31 mpfr_init2 (x, 32); 32 mpfr_init2 (y, 65); 33 mpfr_init2 (z, 17); 34 mpfr_init2 (a, 17); 35 36 mpfr_set_str_binary (x, "-0.101110001001011011011e-9"); 37 mpfr_ui_pow (y, 7, x, MPFR_RNDN); 38 39 mpfr_set_prec (x, 40); 40 mpfr_set_str_binary (x, "-0.1100101100101111011001010010110011110110E-1"); 41 mpfr_set_prec (y, 74); 42 mpfr_ui_pow (y, 8, x, MPFR_RNDN); 43 mpfr_set_prec (x, 74); 44 mpfr_set_str_binary (x, "0.11100000010100111101000011111011011010011000011000101011010011010101000011E-1"); 45 if (mpfr_cmp (x, y)) 46 { 47 printf ("Error for input of 40 bits, output of 74 bits\n"); 48 exit (1); 49 } 50 51 /* Check for ui_pow_ui */ 52 mpfr_ui_pow_ui (x, 0, 1, MPFR_RNDN); 53 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x)); 54 mpfr_ui_pow_ui (x, 0, 4, MPFR_RNDN); 55 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x)); 56 res1 = mpfr_ui_pow_ui (z, 17, 42, MPFR_RNDD); 57 mpfr_set_ui (x, 17, MPFR_RNDN); 58 mpfr_set_ui (y, 42, MPFR_RNDN); 59 res2 = mpfr_pow (a, x, y, MPFR_RNDD); 60 if (mpfr_cmp (z, a) || res1 != res2) 61 { 62 printf ("Error for ui_pow_ui for 17^42\n" 63 "Inexact1 = %d Inexact2 = %d\n", res1, res2); 64 mpfr_dump (z); 65 mpfr_dump (a); 66 exit (1); 67 } 68 mpfr_set_prec (x, 2); 69 mpfr_ui_pow_ui (x, 65537, 65535, MPFR_RNDN); 70 if (mpfr_cmp_str (x, "0.11E1048562", 2, MPFR_RNDN) != 0) 71 { 72 printf ("Error for ui_pow_ui for 65537 ^65535 with 2 bits of precision\n"); 73 mpfr_dump (x); 74 exit (1); 75 } 76 mpfr_clear (x); 77 mpfr_clear (y); 78 mpfr_clear (z); 79 mpfr_clear (a); 80} 81 82static void 83check1 (mpfr_ptr x, mpfr_prec_t prec, unsigned long nt, mpfr_rnd_t rnd) 84{ 85 mpfr_t y, z, t; 86 int inexact, compare, compare2; 87 mpfr_prec_t yprec; 88 mpfr_exp_t err; 89 90 yprec = prec + 10; 91 92 mpfr_init (y); 93 mpfr_init (z); 94 mpfr_init (t); 95 mpfr_set_prec (y, yprec); 96 mpfr_set_prec (z, prec); 97 mpfr_set_prec (t, prec); 98 99 compare = mpfr_ui_pow (y, nt, x, rnd); 100 err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec; 101 if (mpfr_can_round (y, err, rnd, rnd, prec)) 102 { 103 mpfr_set (t, y, rnd); 104 inexact = mpfr_ui_pow (z, nt, x, rnd); 105 if (mpfr_cmp (t, z)) 106 { 107 printf ("results differ for x="); 108 mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN); 109 printf (" n=%lu", nt); 110 printf (" prec=%u rnd_mode=%s\n", (unsigned) prec, 111 mpfr_print_rnd_mode (rnd)); 112 printf ("got "); 113 mpfr_dump (z); 114 printf ("expected "); 115 mpfr_dump (t); 116 printf ("approx "); 117 mpfr_dump (y); 118 exit (1); 119 } 120 compare2 = mpfr_cmp (t, y); 121 /* if rounding to nearest, cannot know the sign of t - f(x) 122 because of composed rounding: y = o(f(x)) and t = o(y) */ 123 if ((rnd != MPFR_RNDN) && (compare * compare2 >= 0)) 124 compare = compare + compare2; 125 else 126 compare = inexact; /* cannot determine sign(t-f(x)) */ 127 if (((inexact == 0) && (compare != 0)) || 128 ((inexact > 0) && (compare <= 0)) || 129 ((inexact < 0) && (compare >= 0))) 130 { 131 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n", 132 mpfr_print_rnd_mode (rnd), compare, inexact); 133 printf ("x="); mpfr_dump (x); 134 printf ("y="); mpfr_dump (y); 135 printf ("t="); mpfr_dump (t); 136 exit (1); 137 } 138 } 139 140 mpfr_clear (y); 141 mpfr_clear (z); 142 mpfr_clear (t); 143} 144 145static void 146huge (void) 147{ 148 mpfr_exp_t old_emin, old_emax; 149 mpfr_t x; 150 151 old_emin = mpfr_get_emin (); 152 old_emax = mpfr_get_emax (); 153 154 set_emin (MPFR_EMIN_MIN); 155 set_emax (MPFR_EMAX_MAX); 156 157 mpfr_init2 (x, 8); 158 159 /* The purpose of this test is more to check that mpfr_ui_pow_ui 160 terminates (without taking much memory) rather than checking 161 the value of x. On 2023-02-13, the +Inf case was not handled 162 in the Ziv iteration, yielding an infinite loop, affecting 163 mpfr_log10 in particular. See 164 commit 90de094f0d9c309daca707aa227470d810866616 165 */ 166 mpfr_ui_pow_ui (x, 5, ULONG_MAX, MPFR_RNDN); 167 if (MPFR_EMAX_MAX <= ULONG_MAX) /* true with default _MPFR_EXP_FORMAT */ 168 MPFR_ASSERTN (MPFR_IS_INF (x)); 169 170 mpfr_clear (x); 171 172 set_emin (old_emin); 173 set_emax (old_emax); 174} 175 176static void 177test2 (void) 178{ 179 mpfr_t x, y, z, t; 180 mpfr_prec_t prec; 181 mpfr_rnd_t rnd; 182 unsigned int n; 183 mpfr_prec_t p0 = 2, p1 = 100; 184 unsigned int N = 20; 185 186 mpfr_init2 (x, 2); 187 mpfr_init2 (y, 2); 188 mpfr_init2 (z, 38); 189 mpfr_init2 (t, 6); 190 191 /* check exact power */ 192 mpfr_set_str_binary (t, "0.110000E5"); 193 mpfr_ui_pow (z, 3, t, MPFR_RNDN); 194 195 mpfr_set_str (x, "-0.5", 10, MPFR_RNDZ); 196 mpfr_ui_pow (y, 4, x, MPFR_RNDD); 197 if (mpfr_cmp_ui_2exp (y, 1, -1)) 198 { 199 printf ("Error for 4^(-0.5), prec=2, MPFR_RNDD\n"); 200 printf ("expected 0.5, got "); 201 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 202 printf ("\n"); 203 exit (1); 204 } 205 206 /* problem found by Kevin on spe175.testdrive.compaq.com 207 (03 Sep 2003), ia64 under HP-UX */ 208 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 209 mpfr_ui_pow (y, 398441521, x, MPFR_RNDN); 210 if (mpfr_cmp_ui_2exp (y, 1, 14)) 211 { 212 printf ("Error for 398441521^(0.5), prec=2, MPFR_RNDN\n"); 213 printf ("expected 1.0e14, got "); 214 mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN); 215 printf ("\n"); 216 exit (1); 217 } 218 219 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 220 check1 (x, 2, 398441521, MPFR_RNDN); /* 398441521 = 19961^2 */ 221 222 /* generic test */ 223 for (prec = p0; prec <= p1; prec++) 224 { 225 mpfr_set_prec (x, prec); 226 for (n = 0; n < N; n++) 227 { 228 int nt; 229 nt = randlimb () & INT_MAX; 230 mpfr_urandomb (x, RANDS); 231 rnd = RND_RAND_NO_RNDF (); 232 check1 (x, prec, nt, rnd); 233 } 234 } 235 236 mpfr_clears (x, y, z, t, (mpfr_ptr) 0); 237} 238 239/* reverse the arguments n and x for tgeneric_ui.c */ 240static int 241ui_pow_rev (mpfr_ptr y, mpfr_srcptr x, unsigned long n, mpfr_rnd_t rnd) 242{ 243 return mpfr_ui_pow (y, n, x, rnd); 244} 245 246#define TEST_FUNCTION ui_pow_rev 247#define INTEGER_TYPE unsigned long 248#define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), 1, RANDS) 249#define INT_RAND_FUNCTION() \ 250 (randlimb () % 16 == 0 ? randulong () : (unsigned long) (randlimb () % 32)) 251#include "tgeneric_ui.c" 252 253int 254main (int argc, char *argv[]) 255{ 256 mpfr_t x, y; 257 unsigned long int n; 258 259 tests_start_mpfr (); 260 261 mpfr_init (x); 262 mpfr_init (y); 263 264 do n = randlimb (); while (n <= 1); 265 266 MPFR_SET_INF (x); 267 mpfr_ui_pow (y, n, x, MPFR_RNDN); 268 if (! MPFR_IS_INF (y)) 269 { 270 printf ("evaluation of function at INF does not return INF\n"); 271 exit (1); 272 } 273 274 MPFR_CHANGE_SIGN (x); 275 mpfr_ui_pow (y, n, x, MPFR_RNDN); 276 if (! MPFR_IS_ZERO (y)) 277 { 278 printf ("evaluation of function in -INF does not return 0\n"); 279 exit (1); 280 } 281 282 MPFR_SET_NAN (x); 283 mpfr_ui_pow (y, n, x, MPFR_RNDN); 284 if (! MPFR_IS_NAN (y)) 285 { 286 printf ("evaluation of function in NAN does not return NAN\n"); 287 exit (1); 288 } 289 290 mpfr_clear (x); 291 mpfr_clear (y); 292 293 test1 (); 294 test2 (); 295 huge (); 296 297 test_generic_ui (MPFR_PREC_MIN, 100, 100); 298 299 tests_end_mpfr (); 300 return 0; 301} 302