1/* Test file for mpfr_sqrt. 2 3Copyright 1999, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc. 4Contributed by the Arenaire and Cacao 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 20http://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 <stdio.h> 24#include <stdlib.h> 25 26#include "mpfr-test.h" 27 28#ifdef CHECK_EXTERNAL 29static int 30test_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 51static void 52check3 (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 71static void 72check4 (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 92static void 93check24 (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 112static void 113check_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, 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 */ 134static void 135check_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 185static void 186special (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 421static void 422check_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 458static void 459check_nan (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_SET_NEG (x); 491 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 492 MPFR_ASSERTN (mpfr_number_p (got)); 493 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); 494 495 mpfr_clear (x); 496 mpfr_clear (got); 497} 498 499/* check that -1 <= x/sqrt(x^2+s*y^2) <= 1 for rounding to nearest or up 500 with s = 0 and s = 1 */ 501static void 502test_property1 (mpfr_prec_t p, mpfr_rnd_t r, int s) 503{ 504 mpfr_t x, y, z, t; 505 506 mpfr_init2 (x, p); 507 mpfr_init2 (y, p); 508 mpfr_init2 (z, p); 509 mpfr_init2 (t, p); 510 511 mpfr_urandomb (x, RANDS); 512 mpfr_mul (z, x, x, r); 513 if (s) 514 { 515 mpfr_urandomb (y, RANDS); 516 mpfr_mul (t, y, y, r); 517 mpfr_add (z, z, t, r); 518 } 519 mpfr_sqrt (z, z, r); 520 mpfr_div (z, x, z, r); 521 /* Note: if both x and y are 0, z is NAN, but the test below will 522 be false. So, everything is fine. */ 523 if (mpfr_cmp_si (z, -1) < 0 || mpfr_cmp_ui (z, 1) > 0) 524 { 525 printf ("Error, -1 <= x/sqrt(x^2+y^2) <= 1 does not hold for r=%s\n", 526 mpfr_print_rnd_mode (r)); 527 printf ("x="); mpfr_dump (x); 528 printf ("y="); mpfr_dump (y); 529 printf ("got "); mpfr_dump (z); 530 exit (1); 531 } 532 533 mpfr_clear (x); 534 mpfr_clear (y); 535 mpfr_clear (z); 536 mpfr_clear (t); 537} 538 539/* check sqrt(x^2) = x */ 540static void 541test_property2 (mpfr_prec_t p, mpfr_rnd_t r) 542{ 543 mpfr_t x, y; 544 545 mpfr_init2 (x, p); 546 mpfr_init2 (y, p); 547 548 mpfr_urandomb (x, RANDS); 549 mpfr_mul (y, x, x, r); 550 mpfr_sqrt (y, y, r); 551 if (mpfr_cmp (y, x)) 552 { 553 printf ("Error, sqrt(x^2) = x does not hold for r=%s\n", 554 mpfr_print_rnd_mode (r)); 555 printf ("x="); mpfr_dump (x); 556 printf ("got "); mpfr_dump (y); 557 exit (1); 558 } 559 560 mpfr_clear (x); 561 mpfr_clear (y); 562} 563 564#define TEST_FUNCTION test_sqrt 565#define TEST_RANDOM_POS 8 566#include "tgeneric.c" 567 568int 569main (void) 570{ 571 mpfr_prec_t p; 572 int k; 573 574 tests_start_mpfr (); 575 576 for (p = MPFR_PREC_MIN; p <= 128; p++) 577 { 578 test_property1 (p, MPFR_RNDN, 0); 579 test_property1 (p, MPFR_RNDU, 0); 580 test_property1 (p, MPFR_RNDA, 0); 581 test_property1 (p, MPFR_RNDN, 1); 582 test_property1 (p, MPFR_RNDU, 1); 583 test_property1 (p, MPFR_RNDA, 1); 584 test_property2 (p, MPFR_RNDN); 585 } 586 587 check_diverse ("635030154261163106768013773815762607450069292760790610550915652722277604820131530404842415587328", 160, "796887792767063979679855997149887366668464780637"); 588 special (); 589 check_nan (); 590 591 for (p=2; p<200; p++) 592 for (k=0; k<200; k++) 593 check_inexact (p); 594 check_float(); 595 596 check3 ("-0.0", MPFR_RNDN, "0.0"); 597 check4 ("6.37983013646045901440e+32", MPFR_RNDN, "5.9bc5036d09e0c@13"); 598 check4 ("1.0", MPFR_RNDN, "1"); 599 check4 ("1.0", MPFR_RNDZ, "1"); 600 check4 ("3.725290298461914062500000e-9", MPFR_RNDN, "4@-4"); 601 check4 ("3.725290298461914062500000e-9", MPFR_RNDZ, "4@-4"); 602 check4 ("1190456976439861.0", MPFR_RNDZ, "2.0e7957873529a@6"); 603 check4 ("1219027943874417664.0", MPFR_RNDZ, "4.1cf2af0e6a534@7"); 604 /* the following examples are bugs in Cygnus compiler/system, found by 605 Fabrice Rouillier while porting mpfr to Windows */ 606 check4 ("9.89438396044940256501e-134", MPFR_RNDU, "8.7af7bf0ebbee@-56"); 607 check4 ("7.86528588050363751914e+31", MPFR_RNDZ, "1.f81fc40f32062@13"); 608 check4 ("0.99999999999999988897", MPFR_RNDN, "f.ffffffffffff8@-1"); 609 check4 ("1.00000000000000022204", MPFR_RNDN, "1"); 610 /* the following examples come from the paper "Number-theoretic Test 611 Generation for Directed Rounding" from Michael Parks, Table 4 */ 612 613 check4 ("78652858805036375191418371571712.0", MPFR_RNDN, 614 "1.f81fc40f32063@13"); 615 check4 ("38510074998589467860312736661504.0", MPFR_RNDN, 616 "1.60c012a92fc65@13"); 617 check4 ("35318779685413012908190921129984.0", MPFR_RNDN, 618 "1.51d17526c7161@13"); 619 check4 ("26729022595358440976973142425600.0", MPFR_RNDN, 620 "1.25e19302f7e51@13"); 621 check4 ("22696567866564242819241453027328.0", MPFR_RNDN, 622 "1.0ecea7dd2ec3d@13"); 623 check4 ("22696888073761729132924856434688.0", MPFR_RNDN, 624 "1.0ecf250e8e921@13"); 625 check4 ("36055652513981905145251657416704.0", MPFR_RNDN, 626 "1.5552f3eedcf33@13"); 627 check4 ("30189856268896404997497182748672.0", MPFR_RNDN, 628 "1.3853ee10c9c99@13"); 629 check4 ("36075288240584711210898775080960.0", MPFR_RNDN, 630 "1.556abe212b56f@13"); 631 check4 ("72154663483843080704304789585920.0", MPFR_RNDN, 632 "1.e2d9a51977e6e@13"); 633 634 check4 ("78652858805036375191418371571712.0", MPFR_RNDZ, 635 "1.f81fc40f32062@13"); 636 check4 ("38510074998589467860312736661504.0", MPFR_RNDZ, 637 "1.60c012a92fc64@13"); 638 check4 ("35318779685413012908190921129984.0", MPFR_RNDZ, "1.51d17526c716@13"); 639 check4 ("26729022595358440976973142425600.0", MPFR_RNDZ, "1.25e19302f7e5@13"); 640 check4 ("22696567866564242819241453027328.0", MPFR_RNDZ, 641 "1.0ecea7dd2ec3c@13"); 642 check4 ("22696888073761729132924856434688.0", MPFR_RNDZ, "1.0ecf250e8e92@13"); 643 check4 ("36055652513981905145251657416704.0", MPFR_RNDZ, 644 "1.5552f3eedcf32@13"); 645 check4 ("30189856268896404997497182748672.0", MPFR_RNDZ, 646 "1.3853ee10c9c98@13"); 647 check4 ("36075288240584711210898775080960.0", MPFR_RNDZ, 648 "1.556abe212b56e@13"); 649 check4 ("72154663483843080704304789585920.0", MPFR_RNDZ, 650 "1.e2d9a51977e6d@13"); 651 652 check4 ("78652858805036375191418371571712.0", MPFR_RNDU, 653 "1.f81fc40f32063@13"); 654 check4 ("38510074998589467860312736661504.0", MPFR_RNDU, 655 "1.60c012a92fc65@13"); 656 check4 ("35318779685413012908190921129984.0", MPFR_RNDU, 657 "1.51d17526c7161@13"); 658 check4 ("26729022595358440976973142425600.0", MPFR_RNDU, 659 "1.25e19302f7e51@13"); 660 check4 ("22696567866564242819241453027328.0", MPFR_RNDU, 661 "1.0ecea7dd2ec3d@13"); 662 check4 ("22696888073761729132924856434688.0", MPFR_RNDU, 663 "1.0ecf250e8e921@13"); 664 check4 ("36055652513981905145251657416704.0", MPFR_RNDU, 665 "1.5552f3eedcf33@13"); 666 check4 ("30189856268896404997497182748672.0", MPFR_RNDU, 667 "1.3853ee10c9c99@13"); 668 check4 ("36075288240584711210898775080960.0", MPFR_RNDU, 669 "1.556abe212b56f@13"); 670 check4 ("72154663483843080704304789585920.0", MPFR_RNDU, 671 "1.e2d9a51977e6e@13"); 672 673 check4 ("78652858805036375191418371571712.0", MPFR_RNDD, 674 "1.f81fc40f32062@13"); 675 check4 ("38510074998589467860312736661504.0", MPFR_RNDD, 676 "1.60c012a92fc64@13"); 677 check4 ("35318779685413012908190921129984.0", MPFR_RNDD, "1.51d17526c716@13"); 678 check4 ("26729022595358440976973142425600.0", MPFR_RNDD, "1.25e19302f7e5@13"); 679 check4 ("22696567866564242819241453027328.0", MPFR_RNDD, 680 "1.0ecea7dd2ec3c@13"); 681 check4 ("22696888073761729132924856434688.0", MPFR_RNDD, "1.0ecf250e8e92@13"); 682 check4 ("36055652513981905145251657416704.0", MPFR_RNDD, 683 "1.5552f3eedcf32@13"); 684 check4 ("30189856268896404997497182748672.0", MPFR_RNDD, 685 "1.3853ee10c9c98@13"); 686 check4 ("36075288240584711210898775080960.0", MPFR_RNDD, 687 "1.556abe212b56e@13"); 688 check4 ("72154663483843080704304789585920.0", MPFR_RNDD, 689 "1.e2d9a51977e6d@13"); 690 691 /* check that rounding away is just rounding toward plus infinity */ 692 check4 ("72154663483843080704304789585920.0", MPFR_RNDA, 693 "1.e2d9a51977e6e@13"); 694 695 test_generic (2, 300, 15); 696 data_check ("data/sqrt", mpfr_sqrt, "mpfr_sqrt"); 697 bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 8, -256, 255, 4, 128, 800, 50); 698 699 tests_end_mpfr (); 700 return 0; 701} 702