| OLD | NEW |
|---|
| (Empty) |
| 1 /* Test file for mpfr_sqrt. | |
| 2 | |
| 3 Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Softwa
re Foundation, Inc. | |
| 4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to | |
| 20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | |
| 21 MA 02110-1301, USA. */ | |
| 22 | |
| 23 #include <stdio.h> | |
| 24 #include <stdlib.h> | |
| 25 | |
| 26 #include "mpfr-test.h" | |
| 27 | |
| 28 #ifdef CHECK_EXTERNAL | |
| 29 static int | |
| 30 test_sqrt (mpfr_ptr a, mpfr_srcptr b, mp_rnd_t rnd_mode) | |
| 31 { | |
| 32 int res; | |
| 33 int ok = rnd_mode == GMP_RNDN && mpfr_number_p (b); | |
| 34 if (ok) | |
| 35 { | |
| 36 mpfr_print_raw (b); | |
| 37 } | |
| 38 res = mpfr_sqrt (a, b, rnd_mode); | |
| 39 if (ok) | |
| 40 { | |
| 41 printf (" "); | |
| 42 mpfr_print_raw (a); | |
| 43 printf ("\n"); | |
| 44 } | |
| 45 return res; | |
| 46 } | |
| 47 #else | |
| 48 #define test_sqrt mpfr_sqrt | |
| 49 #endif | |
| 50 | |
| 51 static void | |
| 52 check3 (const char *as, mp_rnd_t rnd_mode, const char *qs) | |
| 53 { | |
| 54 mpfr_t q; | |
| 55 | |
| 56 mpfr_init2 (q, 53); | |
| 57 mpfr_set_str1 (q, as); | |
| 58 test_sqrt (q, q, rnd_mode); | |
| 59 if (mpfr_cmp_str1 (q, qs) ) | |
| 60 { | |
| 61 printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n", | |
| 62 as, mpfr_print_rnd_mode (rnd_mode)); | |
| 63 printf ("expected sqrt is %s, got ",qs); | |
| 64 mpfr_out_str (stdout, 10, 0, q, GMP_RNDN); | |
| 65 putchar ('\n'); | |
| 66 exit (1); | |
| 67 } | |
| 68 mpfr_clear (q); | |
| 69 } | |
| 70 | |
| 71 static void | |
| 72 check4 (const char *as, mp_rnd_t rnd_mode, const char *Qs) | |
| 73 { | |
| 74 mpfr_t q; | |
| 75 | |
| 76 mpfr_init2 (q, 53); | |
| 77 mpfr_set_str1 (q, as); | |
| 78 test_sqrt (q, q, rnd_mode); | |
| 79 if (mpfr_cmp_str (q, Qs, 16, GMP_RNDN)) | |
| 80 { | |
| 81 printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n", | |
| 82 as, mpfr_print_rnd_mode(rnd_mode)); | |
| 83 printf ("expected "); | |
| 84 mpfr_out_str (stdout, 16, 0, q, GMP_RNDN); | |
| 85 printf ("\ngot %s\n", Qs); | |
| 86 mpfr_clear (q); | |
| 87 exit (1); | |
| 88 } | |
| 89 mpfr_clear (q); | |
| 90 } | |
| 91 | |
| 92 static void | |
| 93 check24 (const char *as, mp_rnd_t rnd_mode, const char *qs) | |
| 94 { | |
| 95 mpfr_t q; | |
| 96 | |
| 97 mpfr_init2 (q, 24); | |
| 98 mpfr_set_str1 (q, as); | |
| 99 test_sqrt (q, q, rnd_mode); | |
| 100 if (mpfr_cmp_str1 (q, qs)) | |
| 101 { | |
| 102 printf ("mpfr_sqrt failed for a=%s, prec=24, rnd_mode=%s\n", | |
| 103 as, mpfr_print_rnd_mode(rnd_mode)); | |
| 104 printf ("expected sqrt is %s, got ",qs); | |
| 105 mpfr_out_str (stdout, 10, 0, q, GMP_RNDN); | |
| 106 printf ("\n"); | |
| 107 exit (1); | |
| 108 } | |
| 109 mpfr_clear (q); | |
| 110 } | |
| 111 | |
| 112 static void | |
| 113 check_diverse (const char *as, mp_prec_t p, const char *qs) | |
| 114 { | |
| 115 mpfr_t q; | |
| 116 | |
| 117 mpfr_init2 (q, p); | |
| 118 mpfr_set_str1 (q, as); | |
| 119 test_sqrt (q, q, GMP_RNDN); | |
| 120 if (mpfr_cmp_str1 (q, qs)) | |
| 121 { | |
| 122 printf ("mpfr_sqrt failed for a=%s, prec=%lu, rnd_mode=%s\n", | |
| 123 as, p, mpfr_print_rnd_mode (GMP_RNDN)); | |
| 124 printf ("expected sqrt is %s, got ", qs); | |
| 125 mpfr_out_str (stdout, 10, 0, q, GMP_RNDN); | |
| 126 printf ("\n"); | |
| 127 exit (1); | |
| 128 } | |
| 129 mpfr_clear (q); | |
| 130 } | |
| 131 | |
| 132 /* the following examples come from the paper "Number-theoretic Test | |
| 133 Generation for Directed Rounding" from Michael Parks, Table 3 */ | |
| 134 static void | |
| 135 check_float (void) | |
| 136 { | |
| 137 check24("70368760954880.0", GMP_RNDN, "8.388609e6"); | |
| 138 check24("281474943156224.0", GMP_RNDN, "1.6777215e7"); | |
| 139 check24("70368777732096.0", GMP_RNDN, "8.388610e6"); | |
| 140 check24("281474909601792.0", GMP_RNDN, "1.6777214e7"); | |
| 141 check24("100216216748032.0", GMP_RNDN, "1.0010805e7"); | |
| 142 check24("120137273311232.0", GMP_RNDN, "1.0960715e7"); | |
| 143 check24("229674600890368.0", GMP_RNDN, "1.5155019e7"); | |
| 144 check24("70368794509312.0", GMP_RNDN, "8.388611e6"); | |
| 145 check24("281474876047360.0", GMP_RNDN, "1.6777213e7"); | |
| 146 check24("91214552498176.0", GMP_RNDN, "9.550631e6"); | |
| 147 | |
| 148 check24("70368760954880.0", GMP_RNDZ, "8.388608e6"); | |
| 149 check24("281474943156224.0", GMP_RNDZ, "1.6777214e7"); | |
| 150 check24("70368777732096.0", GMP_RNDZ, "8.388609e6"); | |
| 151 check24("281474909601792.0", GMP_RNDZ, "1.6777213e7"); | |
| 152 check24("100216216748032.0", GMP_RNDZ, "1.0010805e7"); | |
| 153 check24("120137273311232.0", GMP_RNDZ, "1.0960715e7"); | |
| 154 check24("229674600890368.0", GMP_RNDZ, "1.5155019e7"); | |
| 155 check24("70368794509312.0", GMP_RNDZ, "8.38861e6"); | |
| 156 check24("281474876047360.0", GMP_RNDZ, "1.6777212e7"); | |
| 157 check24("91214552498176.0", GMP_RNDZ, "9.550631e6"); | |
| 158 | |
| 159 check24("70368760954880.0", GMP_RNDU, "8.388609e6"); | |
| 160 check24("281474943156224.0",GMP_RNDU, "1.6777215e7"); | |
| 161 check24("70368777732096.0", GMP_RNDU, "8.388610e6"); | |
| 162 check24("281474909601792.0", GMP_RNDU, "1.6777214e7"); | |
| 163 check24("100216216748032.0", GMP_RNDU, "1.0010806e7"); | |
| 164 check24("120137273311232.0", GMP_RNDU, "1.0960716e7"); | |
| 165 check24("229674600890368.0", GMP_RNDU, "1.515502e7"); | |
| 166 check24("70368794509312.0", GMP_RNDU, "8.388611e6"); | |
| 167 check24("281474876047360.0", GMP_RNDU, "1.6777213e7"); | |
| 168 check24("91214552498176.0", GMP_RNDU, "9.550632e6"); | |
| 169 | |
| 170 check24("70368760954880.0", GMP_RNDD, "8.388608e6"); | |
| 171 check24("281474943156224.0", GMP_RNDD, "1.6777214e7"); | |
| 172 check24("70368777732096.0", GMP_RNDD, "8.388609e6"); | |
| 173 check24("281474909601792.0", GMP_RNDD, "1.6777213e7"); | |
| 174 check24("100216216748032.0", GMP_RNDD, "1.0010805e7"); | |
| 175 check24("120137273311232.0", GMP_RNDD, "1.0960715e7"); | |
| 176 check24("229674600890368.0", GMP_RNDD, "1.5155019e7"); | |
| 177 check24("70368794509312.0", GMP_RNDD, "8.38861e6"); | |
| 178 check24("281474876047360.0", GMP_RNDD, "1.6777212e7"); | |
| 179 check24("91214552498176.0", GMP_RNDD, "9.550631e6"); | |
| 180 } | |
| 181 | |
| 182 static void | |
| 183 special (void) | |
| 184 { | |
| 185 mpfr_t x, y, z; | |
| 186 int inexact; | |
| 187 mp_prec_t p; | |
| 188 | |
| 189 mpfr_init (x); | |
| 190 mpfr_init (y); | |
| 191 mpfr_init (z); | |
| 192 | |
| 193 mpfr_set_prec (x, 64); | |
| 194 mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011
10001011001E-1"); | |
| 195 mpfr_set_prec (y, 32); | |
| 196 test_sqrt (y, x, GMP_RNDN); | |
| 197 if (mpfr_cmp_ui (y, 2405743844UL)) | |
| 198 { | |
| 199 printf ("Error for n^2+n+1/2 with n=2405743843\n"); | |
| 200 exit (1); | |
| 201 } | |
| 202 | |
| 203 mpfr_set_prec (x, 65); | |
| 204 mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011
100010110001E-2"); | |
| 205 mpfr_set_prec (y, 32); | |
| 206 test_sqrt (y, x, GMP_RNDN); | |
| 207 if (mpfr_cmp_ui (y, 2405743844UL)) | |
| 208 { | |
| 209 printf ("Error for n^2+n+1/4 with n=2405743843\n"); | |
| 210 mpfr_dump (y); | |
| 211 exit (1); | |
| 212 } | |
| 213 | |
| 214 mpfr_set_prec (x, 66); | |
| 215 mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011
1000101100011E-3"); | |
| 216 mpfr_set_prec (y, 32); | |
| 217 test_sqrt (y, x, GMP_RNDN); | |
| 218 if (mpfr_cmp_ui (y, 2405743844UL)) | |
| 219 { | |
| 220 printf ("Error for n^2+n+1/4+1/8 with n=2405743843\n"); | |
| 221 mpfr_dump (y); | |
| 222 exit (1); | |
| 223 } | |
| 224 | |
| 225 mpfr_set_prec (x, 66); | |
| 226 mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011
1000101100001E-3"); | |
| 227 mpfr_set_prec (y, 32); | |
| 228 test_sqrt (y, x, GMP_RNDN); | |
| 229 if (mpfr_cmp_ui (y, 2405743843UL)) | |
| 230 { | |
| 231 printf ("Error for n^2+n+1/8 with n=2405743843\n"); | |
| 232 mpfr_dump (y); | |
| 233 exit (1); | |
| 234 } | |
| 235 | |
| 236 mpfr_set_prec (x, 27); | |
| 237 mpfr_set_str_binary (x, "0.110100111010101000010001011"); | |
| 238 if ((inexact = test_sqrt (x, x, GMP_RNDZ)) >= 0) | |
| 239 { | |
| 240 printf ("Wrong inexact flag: expected -1, got %d\n", inexact); | |
| 241 exit (1); | |
| 242 } | |
| 243 | |
| 244 mpfr_set_prec (x, 2); | |
| 245 for (p=2; p<1000; p++) | |
| 246 { | |
| 247 mpfr_set_prec (z, p); | |
| 248 mpfr_set_ui (z, 1, GMP_RNDN); | |
| 249 mpfr_nexttoinf (z); | |
| 250 test_sqrt (x, z, GMP_RNDU); | |
| 251 if (mpfr_cmp_ui_2exp(x, 3, -1)) | |
| 252 { | |
| 253 printf ("Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n", | |
| 254 (unsigned int) p); | |
| 255 printf ("got "); mpfr_print_binary (x); puts (""); | |
| 256 exit (1); | |
| 257 } | |
| 258 } | |
| 259 | |
| 260 /* check inexact flag */ | |
| 261 mpfr_set_prec (x, 5); | |
| 262 mpfr_set_str_binary (x, "1.1001E-2"); | |
| 263 if ((inexact = test_sqrt (x, x, GMP_RNDN))) | |
| 264 { | |
| 265 printf ("Wrong inexact flag: expected 0, got %d\n", inexact); | |
| 266 exit (1); | |
| 267 } | |
| 268 | |
| 269 mpfr_set_prec (x, 2); | |
| 270 mpfr_set_prec (z, 2); | |
| 271 | |
| 272 /* checks the sign is correctly set */ | |
| 273 mpfr_set_si (x, 1, GMP_RNDN); | |
| 274 mpfr_set_si (z, -1, GMP_RNDN); | |
| 275 test_sqrt (z, x, GMP_RNDN); | |
| 276 if (mpfr_cmp_ui (z, 0) < 0) | |
| 277 { | |
| 278 printf ("Error: square root of 1 gives "); | |
| 279 mpfr_print_binary(z); | |
| 280 putchar('\n'); | |
| 281 exit (1); | |
| 282 } | |
| 283 | |
| 284 mpfr_set_prec (x, 192); | |
| 285 mpfr_set_prec (z, 160); | |
| 286 mpfr_set_str_binary (z, "0.101101010000010010010010011001100101110010010000001
1000111011001011101101101110000110100001000100001100001011000E1"); | |
| 287 mpfr_set_prec (x, 160); | |
| 288 test_sqrt(x, z, GMP_RNDN); | |
| 289 test_sqrt(z, x, GMP_RNDN); | |
| 290 | |
| 291 mpfr_set_prec (x, 53); | |
| 292 mpfr_set_str (x, "8093416094703476.0", 10, GMP_RNDN); | |
| 293 mpfr_div_2exp (x, x, 1075, GMP_RNDN); | |
| 294 test_sqrt (x, x, GMP_RNDN); | |
| 295 mpfr_set_str (z, "1e55596835b5ef@-141", 16, GMP_RNDN); | |
| 296 if (mpfr_cmp (x, z)) | |
| 297 { | |
| 298 printf ("Error: square root of 8093416094703476*2^(-1075)\n"); | |
| 299 printf ("expected "); mpfr_dump (z); | |
| 300 printf ("got "); mpfr_dump (x); | |
| 301 exit (1); | |
| 302 } | |
| 303 | |
| 304 mpfr_set_prec (x, 33); | |
| 305 mpfr_set_str_binary (x, "0.111011011011110001100111111001000e-10"); | |
| 306 mpfr_set_prec (z, 157); | |
| 307 inexact = test_sqrt (z, x, GMP_RNDN); | |
| 308 mpfr_set_prec (x, 157); | |
| 309 mpfr_set_str_binary (x, "0.111101101011001011110010111000111000111000011010101
11011010000100111011000111110100001001011110011111100101110010110010110011001011
011010110010000011001101E-5"); | |
| 310 if (mpfr_cmp (x, z)) | |
| 311 { | |
| 312 printf ("Error: square root (1)\n"); | |
| 313 exit (1); | |
| 314 } | |
| 315 if (inexact <= 0) | |
| 316 { | |
| 317 printf ("Error: wrong inexact flag (1)\n"); | |
| 318 exit (1); | |
| 319 } | |
| 320 | |
| 321 /* case prec(result) << prec(input) */ | |
| 322 mpfr_set_prec (z, 2); | |
| 323 for (p = 2; p < 1000; p++) | |
| 324 { | |
| 325 mpfr_set_prec (x, p); | |
| 326 mpfr_set_ui (x, 1, GMP_RNDN); | |
| 327 mpfr_nextabove (x); | |
| 328 /* 1.0 < x <= 1.5 thus 1 < sqrt(x) <= 1.23 */ | |
| 329 inexact = test_sqrt (z, x, GMP_RNDN); | |
| 330 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); | |
| 331 inexact = test_sqrt (z, x, GMP_RNDZ); | |
| 332 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); | |
| 333 inexact = test_sqrt (z, x, GMP_RNDU); | |
| 334 MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0); | |
| 335 inexact = test_sqrt (z, x, GMP_RNDD); | |
| 336 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); | |
| 337 } | |
| 338 | |
| 339 /* corner case rw = 0 in rounding to nearest */ | |
| 340 mpfr_set_prec (z, BITS_PER_MP_LIMB - 1); | |
| 341 mpfr_set_prec (y, BITS_PER_MP_LIMB - 1); | |
| 342 mpfr_set_ui (y, 1, GMP_RNDN); | |
| 343 mpfr_mul_2exp (y, y, BITS_PER_MP_LIMB - 1, GMP_RNDN); | |
| 344 mpfr_nextabove (y); | |
| 345 for (p = 2 * BITS_PER_MP_LIMB - 1; p <= 1000; p++) | |
| 346 { | |
| 347 mpfr_set_prec (x, p); | |
| 348 mpfr_set_ui (x, 1, GMP_RNDN); | |
| 349 mpfr_set_exp (x, BITS_PER_MP_LIMB); | |
| 350 mpfr_add_ui (x, x, 1, GMP_RNDN); | |
| 351 /* now x = 2^(BITS_PER_MP_LIMB - 1) + 1 (BITS_PER_MP_LIMB bits) */ | |
| 352 MPFR_ASSERTN (mpfr_mul (x, x, x, GMP_RNDN) == 0); /* exact */ | |
| 353 inexact = test_sqrt (z, x, GMP_RNDN); | |
| 354 /* even rule: z should be 2^(BITS_PER_MP_LIMB - 1) */ | |
| 355 MPFR_ASSERTN (inexact < 0); | |
| 356 MPFR_ASSERTN (mpfr_cmp_ui_2exp (z, 1, BITS_PER_MP_LIMB - 1) == 0); | |
| 357 mpfr_nextbelow (x); | |
| 358 /* now x is just below [2^(BITS_PER_MP_LIMB - 1) + 1]^2 */ | |
| 359 inexact = test_sqrt (z, x, GMP_RNDN); | |
| 360 MPFR_ASSERTN(inexact < 0 && | |
| 361 mpfr_cmp_ui_2exp (z, 1, BITS_PER_MP_LIMB - 1) == 0); | |
| 362 mpfr_nextabove (x); | |
| 363 mpfr_nextabove (x); | |
| 364 /* now x is just above [2^(BITS_PER_MP_LIMB - 1) + 1]^2 */ | |
| 365 inexact = test_sqrt (z, x, GMP_RNDN); | |
| 366 if (mpfr_cmp (z, y)) | |
| 367 { | |
| 368 printf ("Error for sqrt(x) in rounding to nearest\n"); | |
| 369 printf ("x="); mpfr_dump (x); | |
| 370 printf ("Expected "); mpfr_dump (y); | |
| 371 printf ("Got "); mpfr_dump (z); | |
| 372 exit (1); | |
| 373 } | |
| 374 if (inexact <= 0) | |
| 375 { | |
| 376 printf ("Wrong inexact flag in corner case for p = %lu\n", (unsigned l
ong) p); | |
| 377 exit (1); | |
| 378 } | |
| 379 } | |
| 380 | |
| 381 mpfr_set_prec (x, 1000); | |
| 382 mpfr_set_ui (x, 9, GMP_RNDN); | |
| 383 mpfr_set_prec (y, 10); | |
| 384 inexact = test_sqrt (y, x, GMP_RNDN); | |
| 385 if (mpfr_cmp_ui (y, 3) || inexact != 0) | |
| 386 { | |
| 387 printf ("Error in sqrt(9:1000) for prec=10\n"); | |
| 388 exit (1); | |
| 389 } | |
| 390 mpfr_set_prec (y, BITS_PER_MP_LIMB); | |
| 391 mpfr_nextabove (x); | |
| 392 inexact = test_sqrt (y, x, GMP_RNDN); | |
| 393 if (mpfr_cmp_ui (y, 3) || inexact >= 0) | |
| 394 { | |
| 395 printf ("Error in sqrt(9:1000) for prec=%d\n", (int) BITS_PER_MP_LIMB); | |
| 396 exit (1); | |
| 397 } | |
| 398 mpfr_set_prec (x, 2 * BITS_PER_MP_LIMB); | |
| 399 mpfr_set_prec (y, BITS_PER_MP_LIMB); | |
| 400 mpfr_set_ui (y, 1, GMP_RNDN); | |
| 401 mpfr_nextabove (y); | |
| 402 mpfr_set (x, y, GMP_RNDN); | |
| 403 inexact = test_sqrt (y, x, GMP_RNDN); | |
| 404 if (mpfr_cmp_ui (y, 1) || inexact >= 0) | |
| 405 { | |
| 406 printf ("Error in sqrt(1) for prec=%d\n", (int) BITS_PER_MP_LIMB); | |
| 407 mpfr_dump (y); | |
| 408 exit (1); | |
| 409 } | |
| 410 | |
| 411 mpfr_clear (x); | |
| 412 mpfr_clear (y); | |
| 413 mpfr_clear (z); | |
| 414 } | |
| 415 | |
| 416 static void | |
| 417 check_inexact (mp_prec_t p) | |
| 418 { | |
| 419 mpfr_t x, y, z; | |
| 420 mp_rnd_t rnd; | |
| 421 int inexact, sign; | |
| 422 | |
| 423 mpfr_init2 (x, p); | |
| 424 mpfr_init2 (y, p); | |
| 425 mpfr_init2 (z, 2*p); | |
| 426 mpfr_urandomb (x, RANDS); | |
| 427 rnd = RND_RAND (); | |
| 428 inexact = test_sqrt (y, x, rnd); | |
| 429 if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */ | |
| 430 { | |
| 431 printf ("Error: multiplication should be exact\n"); | |
| 432 exit (1); | |
| 433 } | |
| 434 mpfr_sub (z, z, x, rnd); /* exact also */ | |
| 435 sign = mpfr_cmp_ui (z, 0); | |
| 436 if (((inexact == 0) && (sign)) || | |
| 437 ((inexact > 0) && (sign <= 0)) || | |
| 438 ((inexact < 0) && (sign >= 0))) | |
| 439 { | |
| 440 printf ("Error: wrong inexact flag, expected %d, got %d\n", | |
| 441 sign, inexact); | |
| 442 printf ("x="); | |
| 443 mpfr_print_binary (x); | |
| 444 printf (" rnd=%s\n", mpfr_print_rnd_mode (rnd)); | |
| 445 printf ("y="); mpfr_print_binary (y); puts (""); | |
| 446 exit (1); | |
| 447 } | |
| 448 mpfr_clear (x); | |
| 449 mpfr_clear (y); | |
| 450 mpfr_clear (z); | |
| 451 } | |
| 452 | |
| 453 static void | |
| 454 check_nan (void) | |
| 455 { | |
| 456 mpfr_t x, got; | |
| 457 | |
| 458 mpfr_init2 (x, 100L); | |
| 459 mpfr_init2 (got, 100L); | |
| 460 | |
| 461 /* sqrt(NaN) == NaN */ | |
| 462 MPFR_CLEAR_FLAGS (x); | |
| 463 MPFR_SET_NAN (x); | |
| 464 MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */ | |
| 465 MPFR_ASSERTN (mpfr_nan_p (got)); | |
| 466 | |
| 467 /* sqrt(-1) == NaN */ | |
| 468 mpfr_set_si (x, -1L, GMP_RNDZ); | |
| 469 MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */ | |
| 470 MPFR_ASSERTN (mpfr_nan_p (got)); | |
| 471 | |
| 472 /* sqrt(+inf) == +inf */ | |
| 473 MPFR_CLEAR_FLAGS (x); | |
| 474 MPFR_SET_INF (x); | |
| 475 MPFR_SET_POS (x); | |
| 476 MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */ | |
| 477 MPFR_ASSERTN (mpfr_inf_p (got)); | |
| 478 | |
| 479 /* sqrt(-inf) == NaN */ | |
| 480 MPFR_CLEAR_FLAGS (x); | |
| 481 MPFR_SET_INF (x); | |
| 482 MPFR_SET_NEG (x); | |
| 483 MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */ | |
| 484 MPFR_ASSERTN (mpfr_nan_p (got)); | |
| 485 | |
| 486 /* sqrt(-0) == 0 */ | |
| 487 mpfr_set_si (x, 0L, GMP_RNDZ); | |
| 488 MPFR_SET_NEG (x); | |
| 489 MPFR_ASSERTN (test_sqrt (got, x, GMP_RNDZ) == 0); /* exact */ | |
| 490 MPFR_ASSERTN (mpfr_number_p (got)); | |
| 491 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); | |
| 492 | |
| 493 mpfr_clear (x); | |
| 494 mpfr_clear (got); | |
| 495 } | |
| 496 | |
| 497 /* check that -1 <= x/sqrt(x^2+y^2) <= 1 for rounding to nearest or up */ | |
| 498 static void | |
| 499 test_property1 (mp_prec_t p, mp_rnd_t r) | |
| 500 { | |
| 501 mpfr_t x, y, z, t; | |
| 502 | |
| 503 mpfr_init2 (x, p); | |
| 504 mpfr_init2 (y, p); | |
| 505 mpfr_init2 (z, p); | |
| 506 mpfr_init2 (t, p); | |
| 507 | |
| 508 mpfr_urandomb (x, RANDS); | |
| 509 mpfr_urandomb (y, RANDS); | |
| 510 mpfr_mul (z, x, x, r); | |
| 511 mpfr_mul (t, x, x, r); | |
| 512 mpfr_add (z, z, t, r); | |
| 513 mpfr_sqrt (z, z, r); | |
| 514 mpfr_div (z, x, z, r); | |
| 515 if (mpfr_cmp_si (z, -1) < 0 || mpfr_cmp_ui (z, 1) > 0) | |
| 516 { | |
| 517 printf ("Error, -1 <= x/sqrt(x^2+y^2) <= 1 does not hold for r=%s\n", | |
| 518 mpfr_print_rnd_mode (r)); | |
| 519 printf ("x="); mpfr_dump (x); | |
| 520 printf ("y="); mpfr_dump (y); | |
| 521 printf ("got "); mpfr_dump (z); | |
| 522 exit (1); | |
| 523 } | |
| 524 | |
| 525 mpfr_clear (x); | |
| 526 mpfr_clear (y); | |
| 527 mpfr_clear (z); | |
| 528 mpfr_clear (t); | |
| 529 } | |
| 530 | |
| 531 /* check sqrt(x^2) = x */ | |
| 532 static void | |
| 533 test_property2 (mp_prec_t p, mp_rnd_t r) | |
| 534 { | |
| 535 mpfr_t x, y; | |
| 536 | |
| 537 mpfr_init2 (x, p); | |
| 538 mpfr_init2 (y, p); | |
| 539 | |
| 540 mpfr_urandomb (x, RANDS); | |
| 541 mpfr_mul (y, x, x, r); | |
| 542 mpfr_sqrt (y, y, r); | |
| 543 if (mpfr_cmp (y, x)) | |
| 544 { | |
| 545 printf ("Error, sqrt(x^2) = x does not hold for r=%s\n", | |
| 546 mpfr_print_rnd_mode (r)); | |
| 547 printf ("x="); mpfr_dump (x); | |
| 548 printf ("got "); mpfr_dump (y); | |
| 549 exit (1); | |
| 550 } | |
| 551 | |
| 552 mpfr_clear (x); | |
| 553 mpfr_clear (y); | |
| 554 } | |
| 555 | |
| 556 #define TEST_FUNCTION test_sqrt | |
| 557 #define TEST_RANDOM_POS 8 | |
| 558 #include "tgeneric.c" | |
| 559 | |
| 560 int | |
| 561 main (void) | |
| 562 { | |
| 563 mp_prec_t p; | |
| 564 int k; | |
| 565 | |
| 566 tests_start_mpfr (); | |
| 567 | |
| 568 for (p = MPFR_PREC_MIN; p <= 128; p++) | |
| 569 { | |
| 570 test_property1 (p, GMP_RNDN); | |
| 571 test_property1 (p, GMP_RNDU); | |
| 572 test_property2 (p, GMP_RNDN); | |
| 573 } | |
| 574 | |
| 575 check_diverse ("63503015426116310676801377381576260745006929276079061055091565
2722277604820131530404842415587328", 160, "7968877927670639796798559971498873666
68464780637"); | |
| 576 special (); | |
| 577 check_nan (); | |
| 578 | |
| 579 for (p=2; p<200; p++) | |
| 580 for (k=0; k<200; k++) | |
| 581 check_inexact (p); | |
| 582 check_float(); | |
| 583 | |
| 584 check3 ("-0.0", GMP_RNDN, "0.0"); | |
| 585 check4 ("6.37983013646045901440e+32", GMP_RNDN, "5.9bc5036d09e0c@13"); | |
| 586 check4 ("1.0", GMP_RNDN, "1"); | |
| 587 check4 ("1.0", GMP_RNDZ, "1"); | |
| 588 check4 ("3.725290298461914062500000e-9", GMP_RNDN, "4@-4"); | |
| 589 check4 ("3.725290298461914062500000e-9", GMP_RNDZ, "4@-4"); | |
| 590 check4 ("1190456976439861.0", GMP_RNDZ, "2.0e7957873529a@6"); | |
| 591 check4 ("1219027943874417664.0", GMP_RNDZ, "4.1cf2af0e6a534@7"); | |
| 592 /* the following examples are bugs in Cygnus compiler/system, found by | |
| 593 Fabrice Rouillier while porting mpfr to Windows */ | |
| 594 check4 ("9.89438396044940256501e-134", GMP_RNDU, "8.7af7bf0ebbee@-56"); | |
| 595 check4 ("7.86528588050363751914e+31", GMP_RNDZ, "1.f81fc40f32062@13"); | |
| 596 check4 ("0.99999999999999988897", GMP_RNDN, "f.ffffffffffff8@-1"); | |
| 597 check4 ("1.00000000000000022204", GMP_RNDN, "1"); | |
| 598 /* the following examples come from the paper "Number-theoretic Test | |
| 599 Generation for Directed Rounding" from Michael Parks, Table 4 */ | |
| 600 | |
| 601 check4 ("78652858805036375191418371571712.0", GMP_RNDN, | |
| 602 "1.f81fc40f32063@13"); | |
| 603 check4 ("38510074998589467860312736661504.0", GMP_RNDN, | |
| 604 "1.60c012a92fc65@13"); | |
| 605 check4 ("35318779685413012908190921129984.0", GMP_RNDN, | |
| 606 "1.51d17526c7161@13"); | |
| 607 check4 ("26729022595358440976973142425600.0", GMP_RNDN, | |
| 608 "1.25e19302f7e51@13"); | |
| 609 check4 ("22696567866564242819241453027328.0", GMP_RNDN, | |
| 610 "1.0ecea7dd2ec3d@13"); | |
| 611 check4 ("22696888073761729132924856434688.0", GMP_RNDN, | |
| 612 "1.0ecf250e8e921@13"); | |
| 613 check4 ("36055652513981905145251657416704.0", GMP_RNDN, | |
| 614 "1.5552f3eedcf33@13"); | |
| 615 check4 ("30189856268896404997497182748672.0", GMP_RNDN, | |
| 616 "1.3853ee10c9c99@13"); | |
| 617 check4 ("36075288240584711210898775080960.0", GMP_RNDN, | |
| 618 "1.556abe212b56f@13"); | |
| 619 check4 ("72154663483843080704304789585920.0", GMP_RNDN, | |
| 620 "1.e2d9a51977e6e@13"); | |
| 621 | |
| 622 check4 ("78652858805036375191418371571712.0", GMP_RNDZ, | |
| 623 "1.f81fc40f32062@13"); | |
| 624 check4 ("38510074998589467860312736661504.0", GMP_RNDZ, | |
| 625 "1.60c012a92fc64@13"); | |
| 626 check4 ("35318779685413012908190921129984.0", GMP_RNDZ, "1.51d17526c716@13"); | |
| 627 check4 ("26729022595358440976973142425600.0", GMP_RNDZ, "1.25e19302f7e5@13"); | |
| 628 check4 ("22696567866564242819241453027328.0", GMP_RNDZ, | |
| 629 "1.0ecea7dd2ec3c@13"); | |
| 630 check4 ("22696888073761729132924856434688.0", GMP_RNDZ, "1.0ecf250e8e92@13"); | |
| 631 check4 ("36055652513981905145251657416704.0", GMP_RNDZ, | |
| 632 "1.5552f3eedcf32@13"); | |
| 633 check4 ("30189856268896404997497182748672.0", GMP_RNDZ, | |
| 634 "1.3853ee10c9c98@13"); | |
| 635 check4 ("36075288240584711210898775080960.0", GMP_RNDZ, | |
| 636 "1.556abe212b56e@13"); | |
| 637 check4 ("72154663483843080704304789585920.0", GMP_RNDZ, | |
| 638 "1.e2d9a51977e6d@13"); | |
| 639 | |
| 640 check4 ("78652858805036375191418371571712.0", GMP_RNDU, | |
| 641 "1.f81fc40f32063@13"); | |
| 642 check4 ("38510074998589467860312736661504.0", GMP_RNDU, | |
| 643 "1.60c012a92fc65@13"); | |
| 644 check4 ("35318779685413012908190921129984.0", GMP_RNDU, | |
| 645 "1.51d17526c7161@13"); | |
| 646 check4 ("26729022595358440976973142425600.0", GMP_RNDU, | |
| 647 "1.25e19302f7e51@13"); | |
| 648 check4 ("22696567866564242819241453027328.0", GMP_RNDU, | |
| 649 "1.0ecea7dd2ec3d@13"); | |
| 650 check4 ("22696888073761729132924856434688.0", GMP_RNDU, | |
| 651 "1.0ecf250e8e921@13"); | |
| 652 check4 ("36055652513981905145251657416704.0", GMP_RNDU, | |
| 653 "1.5552f3eedcf33@13"); | |
| 654 check4 ("30189856268896404997497182748672.0", GMP_RNDU, | |
| 655 "1.3853ee10c9c99@13"); | |
| 656 check4 ("36075288240584711210898775080960.0", GMP_RNDU, | |
| 657 "1.556abe212b56f@13"); | |
| 658 check4 ("72154663483843080704304789585920.0", GMP_RNDU, | |
| 659 "1.e2d9a51977e6e@13"); | |
| 660 | |
| 661 check4 ("78652858805036375191418371571712.0", GMP_RNDD, | |
| 662 "1.f81fc40f32062@13"); | |
| 663 check4 ("38510074998589467860312736661504.0", GMP_RNDD, | |
| 664 "1.60c012a92fc64@13"); | |
| 665 check4 ("35318779685413012908190921129984.0", GMP_RNDD, "1.51d17526c716@13"); | |
| 666 check4 ("26729022595358440976973142425600.0", GMP_RNDD, "1.25e19302f7e5@13"); | |
| 667 check4 ("22696567866564242819241453027328.0", GMP_RNDD, | |
| 668 "1.0ecea7dd2ec3c@13"); | |
| 669 check4 ("22696888073761729132924856434688.0", GMP_RNDD, "1.0ecf250e8e92@13"); | |
| 670 check4 ("36055652513981905145251657416704.0", GMP_RNDD, | |
| 671 "1.5552f3eedcf32@13"); | |
| 672 check4 ("30189856268896404997497182748672.0", GMP_RNDD, | |
| 673 "1.3853ee10c9c98@13"); | |
| 674 check4 ("36075288240584711210898775080960.0", GMP_RNDD, | |
| 675 "1.556abe212b56e@13"); | |
| 676 check4 ("72154663483843080704304789585920.0", GMP_RNDD, | |
| 677 "1.e2d9a51977e6d@13"); | |
| 678 | |
| 679 test_generic (2, 300, 15); | |
| 680 data_check ("data/sqrt", mpfr_sqrt, "mpfr_sqrt"); | |
| 681 bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 8, -256, 255, 4, 128, 800, 50); | |
| 682 | |
| 683 tests_end_mpfr (); | |
| 684 return 0; | |
| 685 } | |
| OLD | NEW |
|---|
Chromium Code Reviews
Issue 3050029:
[gcc] GCC 4.5.0=>4.5.1 (Closed)
Base URL: ssh://git@gitrw.chromium.org:9222/nacl-toolchain.git