1 /* Test file for mpfr_ui_pow and mpfr_ui_pow_ui. 2 3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 http://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 <stdio.h> 24 #include <stdlib.h> 25 #include <limits.h> 26 27 #include "mpfr-test.h" 28 29 static void 30 test1 (void) 31 { 32 mpfr_t x, y, z, a; 33 int res1, res2; 34 35 mpfr_init2 (x, 32); 36 mpfr_init2 (y, 65); 37 mpfr_init2 (z, 17); 38 mpfr_init2 (a, 17); 39 40 mpfr_set_str_binary (x, "-0.101110001001011011011e-9"); 41 mpfr_ui_pow (y, 7, x, MPFR_RNDN); 42 43 mpfr_set_prec (x, 40); 44 mpfr_set_str_binary (x, "-0.1100101100101111011001010010110011110110E-1"); 45 mpfr_set_prec (y, 74); 46 mpfr_ui_pow (y, 8, x, MPFR_RNDN); 47 mpfr_set_prec (x, 74); 48 mpfr_set_str_binary (x, "0.11100000010100111101000011111011011010011000011000101011010011010101000011E-1"); 49 if (mpfr_cmp (x, y)) 50 { 51 printf ("Error for input of 40 bits, output of 74 bits\n"); 52 exit (1); 53 } 54 55 /* Check for ui_pow_ui */ 56 mpfr_ui_pow_ui (x, 0, 1, MPFR_RNDN); 57 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x)); 58 mpfr_ui_pow_ui (x, 0, 4, MPFR_RNDN); 59 MPFR_ASSERTN(mpfr_cmp_ui (x, 0) == 0 && MPFR_IS_POS (x)); 60 res1 = mpfr_ui_pow_ui (z, 17, 42, MPFR_RNDD); 61 mpfr_set_ui (x, 17, MPFR_RNDN); 62 mpfr_set_ui (y, 42, MPFR_RNDN); 63 res2 = mpfr_pow (a, x, y, MPFR_RNDD); 64 if (mpfr_cmp (z, a) || res1 != res2) 65 { 66 printf ("Error for ui_pow_ui for 17^42\n" 67 "Inexact1 = %d Inexact2 = %d\n", res1, res2); 68 mpfr_dump (z); 69 mpfr_dump (a); 70 exit (1); 71 } 72 mpfr_set_prec (x, 2); 73 mpfr_ui_pow_ui (x, 65537, 65535, MPFR_RNDN); 74 if (mpfr_cmp_str (x, "0.11E1048562", 2, MPFR_RNDN) != 0) 75 { 76 printf ("Error for ui_pow_ui for 65537 ^65535 with 2 bits of precision\n"); 77 mpfr_dump (x); 78 exit (1); 79 } 80 mpfr_clear (x); 81 mpfr_clear (y); 82 mpfr_clear (z); 83 mpfr_clear (a); 84 } 85 86 static void 87 check1 (mpfr_ptr x, mpfr_prec_t prec, unsigned long nt, mpfr_rnd_t rnd) 88 { 89 mpfr_t y, z, t; 90 int inexact, compare, compare2; 91 mpfr_prec_t yprec; 92 mpfr_exp_t err; 93 94 yprec = prec + 10; 95 96 mpfr_init (y); 97 mpfr_init (z); 98 mpfr_init (t); 99 mpfr_set_prec (y, yprec); 100 mpfr_set_prec (z, prec); 101 mpfr_set_prec (t, prec); 102 103 compare = mpfr_ui_pow (y, nt, x, rnd); 104 err = (rnd == MPFR_RNDN) ? yprec + 1 : yprec; 105 if (mpfr_can_round (y, err, rnd, rnd, prec)) 106 { 107 mpfr_set (t, y, rnd); 108 inexact = mpfr_ui_pow (z, nt, x, rnd); 109 if (mpfr_cmp (t, z)) 110 { 111 printf ("results differ for x="); 112 mpfr_out_str (stdout, 2, prec, x, MPFR_RNDN); 113 printf (" n=%lu", nt); 114 printf (" prec=%u rnd_mode=%s\n", (unsigned) prec, 115 mpfr_print_rnd_mode (rnd)); 116 printf ("got "); 117 mpfr_out_str (stdout, 2, prec, z, MPFR_RNDN); 118 puts (""); 119 printf ("expected "); 120 mpfr_out_str (stdout, 2, prec, t, MPFR_RNDN); 121 puts (""); 122 printf ("approx "); 123 mpfr_print_binary (y); 124 puts (""); 125 exit (1); 126 } 127 compare2 = mpfr_cmp (t, y); 128 /* if rounding to nearest, cannot know the sign of t - f(x) 129 because of composed rounding: y = o(f(x)) and t = o(y) */ 130 if ((rnd != MPFR_RNDN) && (compare * compare2 >= 0)) 131 compare = compare + compare2; 132 else 133 compare = inexact; /* cannot determine sign(t-f(x)) */ 134 if (((inexact == 0) && (compare != 0)) || 135 ((inexact > 0) && (compare <= 0)) || 136 ((inexact < 0) && (compare >= 0))) 137 { 138 printf ("Wrong inexact flag for rnd=%s: expected %d, got %d\n", 139 mpfr_print_rnd_mode (rnd), compare, inexact); 140 printf ("x="); mpfr_print_binary (x); puts (""); 141 printf ("y="); mpfr_print_binary (y); puts (""); 142 printf ("t="); mpfr_print_binary (t); puts (""); 143 exit (1); 144 } 145 } 146 147 mpfr_clear (y); 148 mpfr_clear (z); 149 mpfr_clear (t); 150 } 151 152 int 153 main (int argc, char *argv[]) 154 { 155 mpfr_t x, y; 156 unsigned long int n; 157 158 tests_start_mpfr (); 159 160 mpfr_init (x); 161 mpfr_init (y); 162 163 n = randlimb (); 164 165 MPFR_SET_INF(x); 166 mpfr_ui_pow (y, n, x, MPFR_RNDN); 167 if(!MPFR_IS_INF(y)) 168 { 169 printf ("evaluation of function in INF does not return INF\n"); 170 exit (1); 171 } 172 173 MPFR_CHANGE_SIGN(x); 174 mpfr_ui_pow (y, n, x, MPFR_RNDN); 175 if(!MPFR_IS_ZERO(y)) 176 { 177 printf ("evaluation of function in -INF does not return 0"); 178 exit (1); 179 } 180 181 MPFR_SET_NAN(x); 182 mpfr_ui_pow (y, n, x, MPFR_RNDN); 183 if(!MPFR_IS_NAN(y)) 184 { 185 printf ("evaluation of function in NAN does not return NAN"); 186 exit (1); 187 } 188 189 test1 (); 190 191 { 192 mpfr_t z, t; 193 mpfr_prec_t prec; 194 mpfr_rnd_t rnd; 195 unsigned int n; 196 197 mpfr_prec_t p0=2, p1=100; 198 unsigned int N=20; 199 200 mpfr_init2 (z, 38); 201 mpfr_init2 (t, 6); 202 203 /* check exact power */ 204 mpfr_set_str_binary (t, "0.110000E5"); 205 mpfr_ui_pow (z, 3, t, MPFR_RNDN); 206 207 mpfr_set_prec (x, 2); 208 mpfr_set_prec (y, 2); 209 mpfr_set_str (x, "-0.5", 10, MPFR_RNDZ); 210 mpfr_ui_pow (y, 4, x, MPFR_RNDD); 211 if (mpfr_cmp_ui_2exp(y, 1, -1)) 212 { 213 fprintf (stderr, "Error for 4^(-0.5), prec=2, MPFR_RNDD\n"); 214 fprintf (stderr, "expected 0.5, got "); 215 mpfr_out_str (stderr, 2, 0, y, MPFR_RNDN); 216 fprintf (stderr, "\n"); 217 exit (1); 218 } 219 220 /* problem found by Kevin on spe175.testdrive.compaq.com 221 (03 Sep 2003), ia64 under HP-UX */ 222 mpfr_set_prec (x, 2); 223 mpfr_set_prec (y, 2); 224 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 225 mpfr_ui_pow (y, 398441521, x, MPFR_RNDN); 226 if (mpfr_cmp_ui_2exp(y, 1, 14)) 227 { 228 fprintf (stderr, "Error for 398441521^(0.5), prec=2, MPFR_RNDN\n"); 229 fprintf (stderr, "expected 1.0e14, got "); 230 mpfr_out_str (stderr, 2, 0, y, MPFR_RNDN); 231 fprintf (stderr, "\n"); 232 exit (1); 233 } 234 235 mpfr_clear (z); 236 mpfr_clear (t); 237 238 mpfr_set_prec (x, 2); 239 mpfr_set_str (x, "0.5", 10, MPFR_RNDN); 240 check1 (x, 2, 398441521, MPFR_RNDN); /* 398441521 = 19961^2 */ 241 242 /* generic test */ 243 for (prec = p0; prec <= p1; prec++) 244 { 245 mpfr_set_prec (x, prec); 246 for (n=0; n<N; n++) 247 { 248 int nt; 249 nt = randlimb () & INT_MAX; 250 mpfr_urandomb (x, RANDS); 251 rnd = RND_RAND (); 252 check1 (x, prec, nt, rnd); 253 } 254 } 255 } 256 257 mpfr_clear (x); 258 mpfr_clear (y); 259 260 tests_end_mpfr (); 261 return 0; 262 } 263