1 /* Test file for mpfr_sqrt. 2 3 Copyright 1999, 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 26 #include "mpfr-test.h" 27 28 #ifdef CHECK_EXTERNAL 29 static int 30 test_sqrt (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode) 31 { 32 int res; 33 int ok = rnd_mode == MPFR_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, mpfr_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, MPFR_RNDN); 65 putchar ('\n'); 66 exit (1); 67 } 68 mpfr_clear (q); 69 } 70 71 static void 72 check4 (const char *as, mpfr_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, MPFR_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, MPFR_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, mpfr_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, MPFR_RNDN); 106 printf ("\n"); 107 exit (1); 108 } 109 mpfr_clear (q); 110 } 111 112 static void 113 check_diverse (const char *as, mpfr_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, MPFR_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, (unsigned long) p, mpfr_print_rnd_mode (MPFR_RNDN)); 124 printf ("expected sqrt is %s, got ", qs); 125 mpfr_out_str (stdout, 10, 0, q, MPFR_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", MPFR_RNDN, "8.388609e6"); 138 check24("281474943156224.0", MPFR_RNDN, "1.6777215e7"); 139 check24("70368777732096.0", MPFR_RNDN, "8.388610e6"); 140 check24("281474909601792.0", MPFR_RNDN, "1.6777214e7"); 141 check24("100216216748032.0", MPFR_RNDN, "1.0010805e7"); 142 check24("120137273311232.0", MPFR_RNDN, "1.0960715e7"); 143 check24("229674600890368.0", MPFR_RNDN, "1.5155019e7"); 144 check24("70368794509312.0", MPFR_RNDN, "8.388611e6"); 145 check24("281474876047360.0", MPFR_RNDN, "1.6777213e7"); 146 check24("91214552498176.0", MPFR_RNDN, "9.550631e6"); 147 148 check24("70368760954880.0", MPFR_RNDZ, "8.388608e6"); 149 check24("281474943156224.0", MPFR_RNDZ, "1.6777214e7"); 150 check24("70368777732096.0", MPFR_RNDZ, "8.388609e6"); 151 check24("281474909601792.0", MPFR_RNDZ, "1.6777213e7"); 152 check24("100216216748032.0", MPFR_RNDZ, "1.0010805e7"); 153 check24("120137273311232.0", MPFR_RNDZ, "1.0960715e7"); 154 check24("229674600890368.0", MPFR_RNDZ, "1.5155019e7"); 155 check24("70368794509312.0", MPFR_RNDZ, "8.38861e6"); 156 check24("281474876047360.0", MPFR_RNDZ, "1.6777212e7"); 157 check24("91214552498176.0", MPFR_RNDZ, "9.550631e6"); 158 159 check24("70368760954880.0", MPFR_RNDU, "8.388609e6"); 160 check24("281474943156224.0",MPFR_RNDU, "1.6777215e7"); 161 check24("70368777732096.0", MPFR_RNDU, "8.388610e6"); 162 check24("281474909601792.0", MPFR_RNDU, "1.6777214e7"); 163 check24("100216216748032.0", MPFR_RNDU, "1.0010806e7"); 164 check24("120137273311232.0", MPFR_RNDU, "1.0960716e7"); 165 check24("229674600890368.0", MPFR_RNDU, "1.515502e7"); 166 check24("70368794509312.0", MPFR_RNDU, "8.388611e6"); 167 check24("281474876047360.0", MPFR_RNDU, "1.6777213e7"); 168 check24("91214552498176.0", MPFR_RNDU, "9.550632e6"); 169 170 check24("70368760954880.0", MPFR_RNDD, "8.388608e6"); 171 check24("281474943156224.0", MPFR_RNDD, "1.6777214e7"); 172 check24("70368777732096.0", MPFR_RNDD, "8.388609e6"); 173 check24("281474909601792.0", MPFR_RNDD, "1.6777213e7"); 174 check24("100216216748032.0", MPFR_RNDD, "1.0010805e7"); 175 check24("120137273311232.0", MPFR_RNDD, "1.0960715e7"); 176 check24("229674600890368.0", MPFR_RNDD, "1.5155019e7"); 177 check24("70368794509312.0", MPFR_RNDD, "8.38861e6"); 178 check24("281474876047360.0", MPFR_RNDD, "1.6777212e7"); 179 check24("91214552498176.0", MPFR_RNDD, "9.550631e6"); 180 181 /* check that rounding away is just rounding toward plus infinity */ 182 check24("91214552498176.0", MPFR_RNDA, "9.550632e6"); 183 } 184 185 static void 186 special (void) 187 { 188 mpfr_t x, y, z; 189 int inexact; 190 mpfr_prec_t p; 191 192 mpfr_init (x); 193 mpfr_init (y); 194 mpfr_init (z); 195 196 mpfr_set_prec (x, 64); 197 mpfr_set_str_binary (x, "1010000010100011011001010101010010001100001101011101110001011001E-1"); 198 mpfr_set_prec (y, 32); 199 test_sqrt (y, x, MPFR_RNDN); 200 if (mpfr_cmp_ui (y, 2405743844UL)) 201 { 202 printf ("Error for n^2+n+1/2 with n=2405743843\n"); 203 exit (1); 204 } 205 206 mpfr_set_prec (x, 65); 207 mpfr_set_str_binary (x, "10100000101000110110010101010100100011000011010111011100010110001E-2"); 208 mpfr_set_prec (y, 32); 209 test_sqrt (y, x, MPFR_RNDN); 210 if (mpfr_cmp_ui (y, 2405743844UL)) 211 { 212 printf ("Error for n^2+n+1/4 with n=2405743843\n"); 213 mpfr_dump (y); 214 exit (1); 215 } 216 217 mpfr_set_prec (x, 66); 218 mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100011E-3"); 219 mpfr_set_prec (y, 32); 220 test_sqrt (y, x, MPFR_RNDN); 221 if (mpfr_cmp_ui (y, 2405743844UL)) 222 { 223 printf ("Error for n^2+n+1/4+1/8 with n=2405743843\n"); 224 mpfr_dump (y); 225 exit (1); 226 } 227 228 mpfr_set_prec (x, 66); 229 mpfr_set_str_binary (x, "101000001010001101100101010101001000110000110101110111000101100001E-3"); 230 mpfr_set_prec (y, 32); 231 test_sqrt (y, x, MPFR_RNDN); 232 if (mpfr_cmp_ui (y, 2405743843UL)) 233 { 234 printf ("Error for n^2+n+1/8 with n=2405743843\n"); 235 mpfr_dump (y); 236 exit (1); 237 } 238 239 mpfr_set_prec (x, 27); 240 mpfr_set_str_binary (x, "0.110100111010101000010001011"); 241 if ((inexact = test_sqrt (x, x, MPFR_RNDZ)) >= 0) 242 { 243 printf ("Wrong inexact flag: expected -1, got %d\n", inexact); 244 exit (1); 245 } 246 247 mpfr_set_prec (x, 2); 248 for (p=2; p<1000; p++) 249 { 250 mpfr_set_prec (z, p); 251 mpfr_set_ui (z, 1, MPFR_RNDN); 252 mpfr_nexttoinf (z); 253 test_sqrt (x, z, MPFR_RNDU); 254 if (mpfr_cmp_ui_2exp(x, 3, -1)) 255 { 256 printf ("Error: sqrt(1+ulp(1), up) should give 1.5 (prec=%u)\n", 257 (unsigned int) p); 258 printf ("got "); mpfr_print_binary (x); puts (""); 259 exit (1); 260 } 261 } 262 263 /* check inexact flag */ 264 mpfr_set_prec (x, 5); 265 mpfr_set_str_binary (x, "1.1001E-2"); 266 if ((inexact = test_sqrt (x, x, MPFR_RNDN))) 267 { 268 printf ("Wrong inexact flag: expected 0, got %d\n", inexact); 269 exit (1); 270 } 271 272 mpfr_set_prec (x, 2); 273 mpfr_set_prec (z, 2); 274 275 /* checks the sign is correctly set */ 276 mpfr_set_si (x, 1, MPFR_RNDN); 277 mpfr_set_si (z, -1, MPFR_RNDN); 278 test_sqrt (z, x, MPFR_RNDN); 279 if (mpfr_cmp_ui (z, 0) < 0) 280 { 281 printf ("Error: square root of 1 gives "); 282 mpfr_print_binary(z); 283 putchar('\n'); 284 exit (1); 285 } 286 287 mpfr_set_prec (x, 192); 288 mpfr_set_prec (z, 160); 289 mpfr_set_str_binary (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1"); 290 mpfr_set_prec (x, 160); 291 test_sqrt(x, z, MPFR_RNDN); 292 test_sqrt(z, x, MPFR_RNDN); 293 294 mpfr_set_prec (x, 53); 295 mpfr_set_str (x, "8093416094703476.0", 10, MPFR_RNDN); 296 mpfr_div_2exp (x, x, 1075, MPFR_RNDN); 297 test_sqrt (x, x, MPFR_RNDN); 298 mpfr_set_str (z, "1e55596835b5ef@-141", 16, MPFR_RNDN); 299 if (mpfr_cmp (x, z)) 300 { 301 printf ("Error: square root of 8093416094703476*2^(-1075)\n"); 302 printf ("expected "); mpfr_dump (z); 303 printf ("got "); mpfr_dump (x); 304 exit (1); 305 } 306 307 mpfr_set_prec (x, 33); 308 mpfr_set_str_binary (x, "0.111011011011110001100111111001000e-10"); 309 mpfr_set_prec (z, 157); 310 inexact = test_sqrt (z, x, MPFR_RNDN); 311 mpfr_set_prec (x, 157); 312 mpfr_set_str_binary (x, "0.11110110101100101111001011100011100011100001101010111011010000100111011000111110100001001011110011111100101110010110010110011001011011010110010000011001101E-5"); 313 if (mpfr_cmp (x, z)) 314 { 315 printf ("Error: square root (1)\n"); 316 exit (1); 317 } 318 if (inexact <= 0) 319 { 320 printf ("Error: wrong inexact flag (1)\n"); 321 exit (1); 322 } 323 324 /* case prec(result) << prec(input) */ 325 mpfr_set_prec (z, 2); 326 for (p = 2; p < 1000; p++) 327 { 328 mpfr_set_prec (x, p); 329 mpfr_set_ui (x, 1, MPFR_RNDN); 330 mpfr_nextabove (x); 331 /* 1.0 < x <= 1.5 thus 1 < sqrt(x) <= 1.23 */ 332 inexact = test_sqrt (z, x, MPFR_RNDN); 333 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 334 inexact = test_sqrt (z, x, MPFR_RNDZ); 335 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 336 inexact = test_sqrt (z, x, MPFR_RNDU); 337 MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0); 338 inexact = test_sqrt (z, x, MPFR_RNDD); 339 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 340 inexact = test_sqrt (z, x, MPFR_RNDA); 341 MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0); 342 } 343 344 /* corner case rw = 0 in rounding to nearest */ 345 mpfr_set_prec (z, GMP_NUMB_BITS - 1); 346 mpfr_set_prec (y, GMP_NUMB_BITS - 1); 347 mpfr_set_ui (y, 1, MPFR_RNDN); 348 mpfr_mul_2exp (y, y, GMP_NUMB_BITS - 1, MPFR_RNDN); 349 mpfr_nextabove (y); 350 for (p = 2 * GMP_NUMB_BITS - 1; p <= 1000; p++) 351 { 352 mpfr_set_prec (x, p); 353 mpfr_set_ui (x, 1, MPFR_RNDN); 354 mpfr_set_exp (x, GMP_NUMB_BITS); 355 mpfr_add_ui (x, x, 1, MPFR_RNDN); 356 /* now x = 2^(GMP_NUMB_BITS - 1) + 1 (GMP_NUMB_BITS bits) */ 357 MPFR_ASSERTN (mpfr_mul (x, x, x, MPFR_RNDN) == 0); /* exact */ 358 inexact = test_sqrt (z, x, MPFR_RNDN); 359 /* even rule: z should be 2^(GMP_NUMB_BITS - 1) */ 360 MPFR_ASSERTN (inexact < 0); 361 MPFR_ASSERTN (mpfr_cmp_ui_2exp (z, 1, GMP_NUMB_BITS - 1) == 0); 362 mpfr_nextbelow (x); 363 /* now x is just below [2^(GMP_NUMB_BITS - 1) + 1]^2 */ 364 inexact = test_sqrt (z, x, MPFR_RNDN); 365 MPFR_ASSERTN(inexact < 0 && 366 mpfr_cmp_ui_2exp (z, 1, GMP_NUMB_BITS - 1) == 0); 367 mpfr_nextabove (x); 368 mpfr_nextabove (x); 369 /* now x is just above [2^(GMP_NUMB_BITS - 1) + 1]^2 */ 370 inexact = test_sqrt (z, x, MPFR_RNDN); 371 if (mpfr_cmp (z, y)) 372 { 373 printf ("Error for sqrt(x) in rounding to nearest\n"); 374 printf ("x="); mpfr_dump (x); 375 printf ("Expected "); mpfr_dump (y); 376 printf ("Got "); mpfr_dump (z); 377 exit (1); 378 } 379 if (inexact <= 0) 380 { 381 printf ("Wrong inexact flag in corner case for p = %lu\n", (unsigned long) p); 382 exit (1); 383 } 384 } 385 386 mpfr_set_prec (x, 1000); 387 mpfr_set_ui (x, 9, MPFR_RNDN); 388 mpfr_set_prec (y, 10); 389 inexact = test_sqrt (y, x, MPFR_RNDN); 390 if (mpfr_cmp_ui (y, 3) || inexact != 0) 391 { 392 printf ("Error in sqrt(9:1000) for prec=10\n"); 393 exit (1); 394 } 395 mpfr_set_prec (y, GMP_NUMB_BITS); 396 mpfr_nextabove (x); 397 inexact = test_sqrt (y, x, MPFR_RNDN); 398 if (mpfr_cmp_ui (y, 3) || inexact >= 0) 399 { 400 printf ("Error in sqrt(9:1000) for prec=%d\n", (int) GMP_NUMB_BITS); 401 exit (1); 402 } 403 mpfr_set_prec (x, 2 * GMP_NUMB_BITS); 404 mpfr_set_prec (y, GMP_NUMB_BITS); 405 mpfr_set_ui (y, 1, MPFR_RNDN); 406 mpfr_nextabove (y); 407 mpfr_set (x, y, MPFR_RNDN); 408 inexact = test_sqrt (y, x, MPFR_RNDN); 409 if (mpfr_cmp_ui (y, 1) || inexact >= 0) 410 { 411 printf ("Error in sqrt(1) for prec=%d\n", (int) GMP_NUMB_BITS); 412 mpfr_dump (y); 413 exit (1); 414 } 415 416 mpfr_clear (x); 417 mpfr_clear (y); 418 mpfr_clear (z); 419 } 420 421 static void 422 check_inexact (mpfr_prec_t p) 423 { 424 mpfr_t x, y, z; 425 mpfr_rnd_t rnd; 426 int inexact, sign; 427 428 mpfr_init2 (x, p); 429 mpfr_init2 (y, p); 430 mpfr_init2 (z, 2*p); 431 mpfr_urandomb (x, RANDS); 432 rnd = RND_RAND (); 433 inexact = test_sqrt (y, x, rnd); 434 if (mpfr_mul (z, y, y, rnd)) /* exact since prec(z) = 2*prec(y) */ 435 { 436 printf ("Error: multiplication should be exact\n"); 437 exit (1); 438 } 439 mpfr_sub (z, z, x, rnd); /* exact also */ 440 sign = mpfr_cmp_ui (z, 0); 441 if (((inexact == 0) && (sign)) || 442 ((inexact > 0) && (sign <= 0)) || 443 ((inexact < 0) && (sign >= 0))) 444 { 445 printf ("Error: wrong inexact flag, expected %d, got %d\n", 446 sign, inexact); 447 printf ("x="); 448 mpfr_print_binary (x); 449 printf (" rnd=%s\n", mpfr_print_rnd_mode (rnd)); 450 printf ("y="); mpfr_print_binary (y); puts (""); 451 exit (1); 452 } 453 mpfr_clear (x); 454 mpfr_clear (y); 455 mpfr_clear (z); 456 } 457 458 static void 459 check_singular (void) 460 { 461 mpfr_t x, got; 462 463 mpfr_init2 (x, 100L); 464 mpfr_init2 (got, 100L); 465 466 /* sqrt(NaN) == NaN */ 467 MPFR_SET_NAN (x); 468 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 469 MPFR_ASSERTN (mpfr_nan_p (got)); 470 471 /* sqrt(-1) == NaN */ 472 mpfr_set_si (x, -1L, MPFR_RNDZ); 473 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 474 MPFR_ASSERTN (mpfr_nan_p (got)); 475 476 /* sqrt(+inf) == +inf */ 477 MPFR_SET_INF (x); 478 MPFR_SET_POS (x); 479 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 480 MPFR_ASSERTN (mpfr_inf_p (got)); 481 482 /* sqrt(-inf) == NaN */ 483 MPFR_SET_INF (x); 484 MPFR_SET_NEG (x); 485 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 486 MPFR_ASSERTN (mpfr_nan_p (got)); 487 488 /* sqrt(+0) == +0 */ 489 mpfr_set_si (x, 0L, MPFR_RNDZ); 490 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 491 MPFR_ASSERTN (mpfr_number_p (got)); 492 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); 493 MPFR_ASSERTN (MPFR_IS_POS (got)); 494 495 /* sqrt(-0) == -0 */ 496 mpfr_set_si (x, 0L, MPFR_RNDZ); 497 MPFR_SET_NEG (x); 498 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 499 MPFR_ASSERTN (mpfr_number_p (got)); 500 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); 501 MPFR_ASSERTN (MPFR_IS_NEG (got)); 502 503 mpfr_clear (x); 504 mpfr_clear (got); 505 } 506 507 /* check that -1 <= x/sqrt(x^2+s*y^2) <= 1 for rounding to nearest or up 508 with s = 0 and s = 1 */ 509 static void 510 test_property1 (mpfr_prec_t p, mpfr_rnd_t r, int s) 511 { 512 mpfr_t x, y, z, t; 513 514 mpfr_init2 (x, p); 515 mpfr_init2 (y, p); 516 mpfr_init2 (z, p); 517 mpfr_init2 (t, p); 518 519 mpfr_urandomb (x, RANDS); 520 mpfr_mul (z, x, x, r); 521 if (s) 522 { 523 mpfr_urandomb (y, RANDS); 524 mpfr_mul (t, y, y, r); 525 mpfr_add (z, z, t, r); 526 } 527 mpfr_sqrt (z, z, r); 528 mpfr_div (z, x, z, r); 529 /* Note: if both x and y are 0, z is NAN, but the test below will 530 be false. So, everything is fine. */ 531 if (mpfr_cmp_si (z, -1) < 0 || mpfr_cmp_ui (z, 1) > 0) 532 { 533 printf ("Error, -1 <= x/sqrt(x^2+y^2) <= 1 does not hold for r=%s\n", 534 mpfr_print_rnd_mode (r)); 535 printf ("x="); mpfr_dump (x); 536 printf ("y="); mpfr_dump (y); 537 printf ("got "); mpfr_dump (z); 538 exit (1); 539 } 540 541 mpfr_clear (x); 542 mpfr_clear (y); 543 mpfr_clear (z); 544 mpfr_clear (t); 545 } 546 547 /* check sqrt(x^2) = x */ 548 static void 549 test_property2 (mpfr_prec_t p, mpfr_rnd_t r) 550 { 551 mpfr_t x, y; 552 553 mpfr_init2 (x, p); 554 mpfr_init2 (y, p); 555 556 mpfr_urandomb (x, RANDS); 557 mpfr_mul (y, x, x, r); 558 mpfr_sqrt (y, y, r); 559 if (mpfr_cmp (y, x)) 560 { 561 printf ("Error, sqrt(x^2) = x does not hold for r=%s\n", 562 mpfr_print_rnd_mode (r)); 563 printf ("x="); mpfr_dump (x); 564 printf ("got "); mpfr_dump (y); 565 exit (1); 566 } 567 568 mpfr_clear (x); 569 mpfr_clear (y); 570 } 571 572 #define TEST_FUNCTION test_sqrt 573 #define TEST_RANDOM_POS 8 574 #include "tgeneric.c" 575 576 int 577 main (void) 578 { 579 mpfr_prec_t p; 580 int k; 581 582 tests_start_mpfr (); 583 584 for (p = MPFR_PREC_MIN; p <= 128; p++) 585 { 586 test_property1 (p, MPFR_RNDN, 0); 587 test_property1 (p, MPFR_RNDU, 0); 588 test_property1 (p, MPFR_RNDA, 0); 589 test_property1 (p, MPFR_RNDN, 1); 590 test_property1 (p, MPFR_RNDU, 1); 591 test_property1 (p, MPFR_RNDA, 1); 592 test_property2 (p, MPFR_RNDN); 593 } 594 595 check_diverse ("635030154261163106768013773815762607450069292760790610550915652722277604820131530404842415587328", 160, "796887792767063979679855997149887366668464780637"); 596 special (); 597 check_singular (); 598 599 for (p=2; p<200; p++) 600 for (k=0; k<200; k++) 601 check_inexact (p); 602 check_float(); 603 604 check3 ("-0.0", MPFR_RNDN, "0.0"); 605 check4 ("6.37983013646045901440e+32", MPFR_RNDN, "5.9bc5036d09e0c@13"); 606 check4 ("1.0", MPFR_RNDN, "1"); 607 check4 ("1.0", MPFR_RNDZ, "1"); 608 check4 ("3.725290298461914062500000e-9", MPFR_RNDN, "4@-4"); 609 check4 ("3.725290298461914062500000e-9", MPFR_RNDZ, "4@-4"); 610 check4 ("1190456976439861.0", MPFR_RNDZ, "2.0e7957873529a@6"); 611 check4 ("1219027943874417664.0", MPFR_RNDZ, "4.1cf2af0e6a534@7"); 612 /* the following examples are bugs in Cygnus compiler/system, found by 613 Fabrice Rouillier while porting mpfr to Windows */ 614 check4 ("9.89438396044940256501e-134", MPFR_RNDU, "8.7af7bf0ebbee@-56"); 615 check4 ("7.86528588050363751914e+31", MPFR_RNDZ, "1.f81fc40f32062@13"); 616 check4 ("0.99999999999999988897", MPFR_RNDN, "f.ffffffffffff8@-1"); 617 check4 ("1.00000000000000022204", MPFR_RNDN, "1"); 618 /* the following examples come from the paper "Number-theoretic Test 619 Generation for Directed Rounding" from Michael Parks, Table 4 */ 620 621 check4 ("78652858805036375191418371571712.0", MPFR_RNDN, 622 "1.f81fc40f32063@13"); 623 check4 ("38510074998589467860312736661504.0", MPFR_RNDN, 624 "1.60c012a92fc65@13"); 625 check4 ("35318779685413012908190921129984.0", MPFR_RNDN, 626 "1.51d17526c7161@13"); 627 check4 ("26729022595358440976973142425600.0", MPFR_RNDN, 628 "1.25e19302f7e51@13"); 629 check4 ("22696567866564242819241453027328.0", MPFR_RNDN, 630 "1.0ecea7dd2ec3d@13"); 631 check4 ("22696888073761729132924856434688.0", MPFR_RNDN, 632 "1.0ecf250e8e921@13"); 633 check4 ("36055652513981905145251657416704.0", MPFR_RNDN, 634 "1.5552f3eedcf33@13"); 635 check4 ("30189856268896404997497182748672.0", MPFR_RNDN, 636 "1.3853ee10c9c99@13"); 637 check4 ("36075288240584711210898775080960.0", MPFR_RNDN, 638 "1.556abe212b56f@13"); 639 check4 ("72154663483843080704304789585920.0", MPFR_RNDN, 640 "1.e2d9a51977e6e@13"); 641 642 check4 ("78652858805036375191418371571712.0", MPFR_RNDZ, 643 "1.f81fc40f32062@13"); 644 check4 ("38510074998589467860312736661504.0", MPFR_RNDZ, 645 "1.60c012a92fc64@13"); 646 check4 ("35318779685413012908190921129984.0", MPFR_RNDZ, "1.51d17526c716@13"); 647 check4 ("26729022595358440976973142425600.0", MPFR_RNDZ, "1.25e19302f7e5@13"); 648 check4 ("22696567866564242819241453027328.0", MPFR_RNDZ, 649 "1.0ecea7dd2ec3c@13"); 650 check4 ("22696888073761729132924856434688.0", MPFR_RNDZ, "1.0ecf250e8e92@13"); 651 check4 ("36055652513981905145251657416704.0", MPFR_RNDZ, 652 "1.5552f3eedcf32@13"); 653 check4 ("30189856268896404997497182748672.0", MPFR_RNDZ, 654 "1.3853ee10c9c98@13"); 655 check4 ("36075288240584711210898775080960.0", MPFR_RNDZ, 656 "1.556abe212b56e@13"); 657 check4 ("72154663483843080704304789585920.0", MPFR_RNDZ, 658 "1.e2d9a51977e6d@13"); 659 660 check4 ("78652858805036375191418371571712.0", MPFR_RNDU, 661 "1.f81fc40f32063@13"); 662 check4 ("38510074998589467860312736661504.0", MPFR_RNDU, 663 "1.60c012a92fc65@13"); 664 check4 ("35318779685413012908190921129984.0", MPFR_RNDU, 665 "1.51d17526c7161@13"); 666 check4 ("26729022595358440976973142425600.0", MPFR_RNDU, 667 "1.25e19302f7e51@13"); 668 check4 ("22696567866564242819241453027328.0", MPFR_RNDU, 669 "1.0ecea7dd2ec3d@13"); 670 check4 ("22696888073761729132924856434688.0", MPFR_RNDU, 671 "1.0ecf250e8e921@13"); 672 check4 ("36055652513981905145251657416704.0", MPFR_RNDU, 673 "1.5552f3eedcf33@13"); 674 check4 ("30189856268896404997497182748672.0", MPFR_RNDU, 675 "1.3853ee10c9c99@13"); 676 check4 ("36075288240584711210898775080960.0", MPFR_RNDU, 677 "1.556abe212b56f@13"); 678 check4 ("72154663483843080704304789585920.0", MPFR_RNDU, 679 "1.e2d9a51977e6e@13"); 680 681 check4 ("78652858805036375191418371571712.0", MPFR_RNDD, 682 "1.f81fc40f32062@13"); 683 check4 ("38510074998589467860312736661504.0", MPFR_RNDD, 684 "1.60c012a92fc64@13"); 685 check4 ("35318779685413012908190921129984.0", MPFR_RNDD, "1.51d17526c716@13"); 686 check4 ("26729022595358440976973142425600.0", MPFR_RNDD, "1.25e19302f7e5@13"); 687 check4 ("22696567866564242819241453027328.0", MPFR_RNDD, 688 "1.0ecea7dd2ec3c@13"); 689 check4 ("22696888073761729132924856434688.0", MPFR_RNDD, "1.0ecf250e8e92@13"); 690 check4 ("36055652513981905145251657416704.0", MPFR_RNDD, 691 "1.5552f3eedcf32@13"); 692 check4 ("30189856268896404997497182748672.0", MPFR_RNDD, 693 "1.3853ee10c9c98@13"); 694 check4 ("36075288240584711210898775080960.0", MPFR_RNDD, 695 "1.556abe212b56e@13"); 696 check4 ("72154663483843080704304789585920.0", MPFR_RNDD, 697 "1.e2d9a51977e6d@13"); 698 699 /* check that rounding away is just rounding toward plus infinity */ 700 check4 ("72154663483843080704304789585920.0", MPFR_RNDA, 701 "1.e2d9a51977e6e@13"); 702 703 test_generic (2, 300, 15); 704 data_check ("data/sqrt", mpfr_sqrt, "mpfr_sqrt"); 705 bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 8, -256, 255, 4, 128, 800, 50); 706 707 tests_end_mpfr (); 708 return 0; 709 } 710