1 /* Test file for mpfr_ui_pow and mpfr_ui_pow_ui. 2 3 Copyright 2001-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 void 26 test1 (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 82 static void 83 check1 (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 145 int 146 main (int argc, char *argv[]) 147 { 148 mpfr_t x, y; 149 unsigned long int n; 150 151 tests_start_mpfr (); 152 153 mpfr_init (x); 154 mpfr_init (y); 155 156 do n = randlimb (); while (n <= 1); 157 158 MPFR_SET_INF(x); 159 mpfr_ui_pow (y, n, x, MPFR_RNDN); 160 if(!MPFR_IS_INF(y)) 161 { 162 printf ("evaluation of function at INF does not return INF\n"); 163 exit (1); 164 } 165 166 MPFR_CHANGE_SIGN(x); 167 mpfr_ui_pow (y, n, x, MPFR_RNDN); 168 if(!MPFR_IS_ZERO(y)) 169 { 170 printf ("evaluation of function in -INF does not return 0"); 171 exit (1); 172 } 173 174 MPFR_SET_NAN(x); 175 mpfr_ui_pow (y, n, x, MPFR_RNDN); 176 if(!MPFR_IS_NAN(y)) 177 { 178 printf ("evaluation of function in NAN does not return NAN"); 179 exit (1); 180 } 181 182 test1 (); 183 184 { 185 mpfr_t z, t; 186 mpfr_prec_t prec; 187 mpfr_rnd_t rnd; 188 unsigned int n; 189 190 mpfr_prec_t p0=2, p1=100; 191 unsigned int N=20; 192 193 mpfr_init2 (z, 38); 194 mpfr_init2 (t, 6); 195 196 /* check exact power */ 197 mpfr_set_str_binary (t, "0.110000E5"); 198 mpfr_ui_pow (z, 3, t, MPFR_RNDN); 199 200 mpfr_set_prec (x, 2); 201 mpfr_set_prec (y, 2); 202 mpfr_set_str (x, "-0.5", 10, MPFR_RNDZ); 203 mpfr_ui_pow (y, 4, x, MPFR_RNDD); 204 if (mpfr_cmp_ui_2exp(y, 1, -1)) 205 { 206 fprintf (stderr, "Error for 4^(-0.5), prec=2, MPFR_RNDD\n"); 207 fprintf (stderr, "expected 0.5, got "); 208 mpfr_out_str (stderr, 2, 0, y, MPFR_RNDN); 209 fprintf (stderr, "\n"); 210 exit (1); 211 } 212 213 /* problem found by Kevin on spe175.testdrive.compaq.com 214 (03 Sep 2003), ia64 under HP-UX */ 215 mpfr_set_prec (x, 2); 216 mpfr_set_prec (y, 2); 217 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 218 mpfr_ui_pow (y, 398441521, x, MPFR_RNDN); 219 if (mpfr_cmp_ui_2exp(y, 1, 14)) 220 { 221 fprintf (stderr, "Error for 398441521^(0.5), prec=2, MPFR_RNDN\n"); 222 fprintf (stderr, "expected 1.0e14, got "); 223 mpfr_out_str (stderr, 2, 0, y, MPFR_RNDN); 224 fprintf (stderr, "\n"); 225 exit (1); 226 } 227 228 mpfr_clear (z); 229 mpfr_clear (t); 230 231 mpfr_set_prec (x, 2); 232 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 233 check1 (x, 2, 398441521, MPFR_RNDN); /* 398441521 = 19961^2 */ 234 235 /* generic test */ 236 for (prec = p0; prec <= p1; prec++) 237 { 238 mpfr_set_prec (x, prec); 239 for (n = 0; n < N; n++) 240 { 241 int nt; 242 nt = randlimb () & INT_MAX; 243 mpfr_urandomb (x, RANDS); 244 rnd = RND_RAND_NO_RNDF (); 245 check1 (x, prec, nt, rnd); 246 } 247 } 248 } 249 250 mpfr_clear (x); 251 mpfr_clear (y); 252 253 tests_end_mpfr (); 254 return 0; 255 } 256