1/* Test file for mpfr_sqrt. 2 3Copyright 1999, 2001-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 "mpfr-test.h" 24 25#ifdef CHECK_EXTERNAL 26static int 27test_sqrt (mpfr_ptr a, mpfr_srcptr b, mpfr_rnd_t rnd_mode) 28{ 29 int res; 30 int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b); 31 if (ok) 32 { 33 mpfr_print_raw (b); 34 } 35 res = mpfr_sqrt (a, b, rnd_mode); 36 if (ok) 37 { 38 printf (" "); 39 mpfr_print_raw (a); 40 printf ("\n"); 41 } 42 return res; 43} 44#else 45#define test_sqrt mpfr_sqrt 46#endif 47 48static void 49check3 (const char *as, mpfr_rnd_t rnd_mode, const char *qs) 50{ 51 mpfr_t q; 52 53 mpfr_init2 (q, 53); 54 mpfr_set_str1 (q, as); 55 test_sqrt (q, q, rnd_mode); 56 if (mpfr_cmp_str1 (q, qs) ) 57 { 58 printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n", 59 as, mpfr_print_rnd_mode (rnd_mode)); 60 printf ("expected sqrt is %s, got ",qs); 61 mpfr_out_str (stdout, 10, 0, q, MPFR_RNDN); 62 putchar ('\n'); 63 exit (1); 64 } 65 mpfr_clear (q); 66} 67 68static void 69check4 (const char *as, mpfr_rnd_t rnd_mode, const char *qs) 70{ 71 mpfr_t q; 72 73 mpfr_init2 (q, 53); 74 mpfr_set_str1 (q, as); 75 test_sqrt (q, q, rnd_mode); 76 if (mpfr_cmp_str (q, qs, 16, MPFR_RNDN)) 77 { 78 printf ("mpfr_sqrt failed for a=%s, rnd_mode=%s\n", 79 as, mpfr_print_rnd_mode (rnd_mode)); 80 printf ("expected %s\ngot ", qs); 81 mpfr_out_str (stdout, 16, 0, q, MPFR_RNDN); 82 printf ("\n"); 83 mpfr_clear (q); 84 exit (1); 85 } 86 mpfr_clear (q); 87} 88 89static void 90check24 (const char *as, mpfr_rnd_t rnd_mode, const char *qs) 91{ 92 mpfr_t q; 93 94 mpfr_init2 (q, 24); 95 mpfr_set_str1 (q, as); 96 test_sqrt (q, q, rnd_mode); 97 if (mpfr_cmp_str1 (q, qs)) 98 { 99 printf ("mpfr_sqrt failed for a=%s, prec=24, rnd_mode=%s\n", 100 as, mpfr_print_rnd_mode(rnd_mode)); 101 printf ("expected sqrt is %s, got ",qs); 102 mpfr_out_str (stdout, 10, 0, q, MPFR_RNDN); 103 printf ("\n"); 104 exit (1); 105 } 106 mpfr_clear (q); 107} 108 109static void 110check_diverse (const char *as, mpfr_prec_t p, const char *qs) 111{ 112 mpfr_t q; 113 114 mpfr_init2 (q, p); 115 mpfr_set_str1 (q, as); 116 test_sqrt (q, q, MPFR_RNDN); 117 if (mpfr_cmp_str1 (q, qs)) 118 { 119 printf ("mpfr_sqrt failed for a=%s, prec=%lu, rnd_mode=%s\n", 120 as, (unsigned long) p, mpfr_print_rnd_mode (MPFR_RNDN)); 121 printf ("expected sqrt is %s, got ", qs); 122 mpfr_out_str (stdout, 10, 0, q, MPFR_RNDN); 123 printf ("\n"); 124 exit (1); 125 } 126 mpfr_clear (q); 127} 128 129/* the following examples come from the paper "Number-theoretic Test 130 Generation for Directed Rounding" from Michael Parks, Table 3 */ 131static void 132check_float (void) 133{ 134 check24("70368760954880.0", MPFR_RNDN, "8.388609e6"); 135 check24("281474943156224.0", MPFR_RNDN, "1.6777215e7"); 136 check24("70368777732096.0", MPFR_RNDN, "8.388610e6"); 137 check24("281474909601792.0", MPFR_RNDN, "1.6777214e7"); 138 check24("100216216748032.0", MPFR_RNDN, "1.0010805e7"); 139 check24("120137273311232.0", MPFR_RNDN, "1.0960715e7"); 140 check24("229674600890368.0", MPFR_RNDN, "1.5155019e7"); 141 check24("70368794509312.0", MPFR_RNDN, "8.388611e6"); 142 check24("281474876047360.0", MPFR_RNDN, "1.6777213e7"); 143 check24("91214552498176.0", MPFR_RNDN, "9.550631e6"); 144 145 check24("70368760954880.0", MPFR_RNDZ, "8.388608e6"); 146 check24("281474943156224.0", MPFR_RNDZ, "1.6777214e7"); 147 check24("70368777732096.0", MPFR_RNDZ, "8.388609e6"); 148 check24("281474909601792.0", MPFR_RNDZ, "1.6777213e7"); 149 check24("100216216748032.0", MPFR_RNDZ, "1.0010805e7"); 150 check24("120137273311232.0", MPFR_RNDZ, "1.0960715e7"); 151 check24("229674600890368.0", MPFR_RNDZ, "1.5155019e7"); 152 check24("70368794509312.0", MPFR_RNDZ, "8.38861e6"); 153 check24("281474876047360.0", MPFR_RNDZ, "1.6777212e7"); 154 check24("91214552498176.0", MPFR_RNDZ, "9.550631e6"); 155 156 check24("70368760954880.0", MPFR_RNDU, "8.388609e6"); 157 check24("281474943156224.0",MPFR_RNDU, "1.6777215e7"); 158 check24("70368777732096.0", MPFR_RNDU, "8.388610e6"); 159 check24("281474909601792.0", MPFR_RNDU, "1.6777214e7"); 160 check24("100216216748032.0", MPFR_RNDU, "1.0010806e7"); 161 check24("120137273311232.0", MPFR_RNDU, "1.0960716e7"); 162 check24("229674600890368.0", MPFR_RNDU, "1.515502e7"); 163 check24("70368794509312.0", MPFR_RNDU, "8.388611e6"); 164 check24("281474876047360.0", MPFR_RNDU, "1.6777213e7"); 165 check24("91214552498176.0", MPFR_RNDU, "9.550632e6"); 166 167 check24("70368760954880.0", MPFR_RNDD, "8.388608e6"); 168 check24("281474943156224.0", MPFR_RNDD, "1.6777214e7"); 169 check24("70368777732096.0", MPFR_RNDD, "8.388609e6"); 170 check24("281474909601792.0", MPFR_RNDD, "1.6777213e7"); 171 check24("100216216748032.0", MPFR_RNDD, "1.0010805e7"); 172 check24("120137273311232.0", MPFR_RNDD, "1.0960715e7"); 173 check24("229674600890368.0", MPFR_RNDD, "1.5155019e7"); 174 check24("70368794509312.0", MPFR_RNDD, "8.38861e6"); 175 check24("281474876047360.0", MPFR_RNDD, "1.6777212e7"); 176 check24("91214552498176.0", MPFR_RNDD, "9.550631e6"); 177 178 /* check that rounding away is just rounding toward positive infinity */ 179 check24("91214552498176.0", MPFR_RNDA, "9.550632e6"); 180} 181 182static void 183special (void) 184{ 185 mpfr_t x, y, z; 186 int inexact; 187 mpfr_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, "1010000010100011011001010101010010001100001101011101110001011001E-1"); 195 mpfr_set_prec (y, 32); 196 test_sqrt (y, x, MPFR_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, "10100000101000110110010101010100100011000011010111011100010110001E-2"); 205 mpfr_set_prec (y, 32); 206 test_sqrt (y, x, MPFR_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, "101000001010001101100101010101001000110000110101110111000101100011E-3"); 216 mpfr_set_prec (y, 32); 217 test_sqrt (y, x, MPFR_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, "101000001010001101100101010101001000110000110101110111000101100001E-3"); 227 mpfr_set_prec (y, 32); 228 test_sqrt (y, x, MPFR_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, MPFR_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, MPFR_RNDN); 249 mpfr_nexttoinf (z); 250 test_sqrt (x, z, MPFR_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_dump (x); 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, MPFR_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, MPFR_RNDN); 274 mpfr_set_si (z, -1, MPFR_RNDN); 275 test_sqrt (z, x, MPFR_RNDN); 276 if (mpfr_cmp_ui (z, 0) < 0) 277 { 278 printf ("Error: square root of 1 gives "); 279 mpfr_dump (z); 280 exit (1); 281 } 282 283 mpfr_set_prec (x, 192); 284 mpfr_set_prec (z, 160); 285 mpfr_set_str_binary (z, "0.1011010100000100100100100110011001011100100100000011000111011001011101101101110000110100001000100001100001011000E1"); 286 mpfr_set_prec (x, 160); 287 test_sqrt(x, z, MPFR_RNDN); 288 test_sqrt(z, x, MPFR_RNDN); 289 290 mpfr_set_prec (x, 53); 291 mpfr_set_str (x, "8093416094703476.0", 10, MPFR_RNDN); 292 mpfr_div_2ui (x, x, 1075, MPFR_RNDN); 293 test_sqrt (x, x, MPFR_RNDN); 294 mpfr_set_str (z, "1e55596835b5ef@-141", 16, MPFR_RNDN); 295 if (mpfr_cmp (x, z)) 296 { 297 printf ("Error: square root of 8093416094703476*2^(-1075)\n"); 298 printf ("expected "); mpfr_dump (z); 299 printf ("got "); mpfr_dump (x); 300 exit (1); 301 } 302 303 mpfr_set_prec (x, 33); 304 mpfr_set_str_binary (x, "0.111011011011110001100111111001000e-10"); 305 mpfr_set_prec (z, 157); 306 inexact = test_sqrt (z, x, MPFR_RNDN); 307 mpfr_set_prec (x, 157); 308 mpfr_set_str_binary (x, "0.11110110101100101111001011100011100011100001101010111011010000100111011000111110100001001011110011111100101110010110010110011001011011010110010000011001101E-5"); 309 if (mpfr_cmp (x, z)) 310 { 311 printf ("Error: square root (1)\n"); 312 exit (1); 313 } 314 if (inexact <= 0) 315 { 316 printf ("Error: wrong inexact flag (1)\n"); 317 exit (1); 318 } 319 320 /* case prec(result) << prec(input) */ 321 mpfr_set_prec (z, 2); 322 for (p = mpfr_get_prec (z); p < 1000; p++) 323 { 324 mpfr_set_prec (x, p); 325 mpfr_set_ui (x, 1, MPFR_RNDN); 326 mpfr_nextabove (x); 327 /* 1.0 < x <= 1.5 thus 1 < sqrt(x) <= 1.23 */ 328 inexact = test_sqrt (z, x, MPFR_RNDN); 329 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 330 inexact = test_sqrt (z, x, MPFR_RNDZ); 331 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 332 inexact = test_sqrt (z, x, MPFR_RNDU); 333 MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0); 334 inexact = test_sqrt (z, x, MPFR_RNDD); 335 MPFR_ASSERTN(inexact < 0 && mpfr_cmp_ui (z, 1) == 0); 336 inexact = test_sqrt (z, x, MPFR_RNDA); 337 MPFR_ASSERTN(inexact > 0 && mpfr_cmp_ui_2exp (z, 3, -1) == 0); 338 } 339 340 /* corner case rw = 0 in rounding to nearest */ 341 mpfr_set_prec (z, GMP_NUMB_BITS - 1); 342 mpfr_set_prec (y, GMP_NUMB_BITS - 1); 343 mpfr_set_ui (y, 1, MPFR_RNDN); 344 mpfr_mul_2ui (y, y, GMP_NUMB_BITS - 1, MPFR_RNDN); 345 mpfr_nextabove (y); 346 for (p = 2 * GMP_NUMB_BITS - 1; p <= 1000; p++) 347 { 348 mpfr_set_prec (x, p); 349 mpfr_set_ui (x, 1, MPFR_RNDN); 350 mpfr_set_exp (x, GMP_NUMB_BITS); 351 mpfr_add_ui (x, x, 1, MPFR_RNDN); 352 /* now x = 2^(GMP_NUMB_BITS - 1) + 1 (GMP_NUMB_BITS bits) */ 353 inexact = mpfr_sqr (x, x, MPFR_RNDN); 354 MPFR_ASSERTN (inexact == 0); /* exact */ 355 inexact = test_sqrt (z, x, MPFR_RNDN); 356 /* even rule: z should be 2^(GMP_NUMB_BITS - 1) */ 357 MPFR_ASSERTN (inexact < 0); 358 inexact = mpfr_cmp_ui_2exp (z, 1, GMP_NUMB_BITS - 1); 359 MPFR_ASSERTN (inexact == 0); 360 mpfr_nextbelow (x); 361 /* now x is just below [2^(GMP_NUMB_BITS - 1) + 1]^2 */ 362 inexact = test_sqrt (z, x, MPFR_RNDN); 363 MPFR_ASSERTN(inexact < 0 && 364 mpfr_cmp_ui_2exp (z, 1, GMP_NUMB_BITS - 1) == 0); 365 mpfr_nextabove (x); 366 mpfr_nextabove (x); 367 /* now x is just above [2^(GMP_NUMB_BITS - 1) + 1]^2 */ 368 inexact = test_sqrt (z, x, MPFR_RNDN); 369 if (mpfr_cmp (z, y)) 370 { 371 printf ("Error for sqrt(x) in rounding to nearest\n"); 372 printf ("x="); mpfr_dump (x); 373 printf ("Expected "); mpfr_dump (y); 374 printf ("Got "); mpfr_dump (z); 375 exit (1); 376 } 377 if (inexact <= 0) 378 { 379 printf ("Wrong inexact flag in corner case for p = %lu\n", (unsigned long) p); 380 exit (1); 381 } 382 } 383 384 mpfr_set_prec (x, 1000); 385 mpfr_set_ui (x, 9, MPFR_RNDN); 386 mpfr_set_prec (y, 10); 387 inexact = test_sqrt (y, x, MPFR_RNDN); 388 if (mpfr_cmp_ui (y, 3) || inexact != 0) 389 { 390 printf ("Error in sqrt(9:1000) for prec=10\n"); 391 exit (1); 392 } 393 mpfr_set_prec (y, GMP_NUMB_BITS); 394 mpfr_nextabove (x); 395 inexact = test_sqrt (y, x, MPFR_RNDN); 396 if (mpfr_cmp_ui (y, 3) || inexact >= 0) 397 { 398 printf ("Error in sqrt(9:1000) for prec=%d\n", (int) GMP_NUMB_BITS); 399 exit (1); 400 } 401 mpfr_set_prec (x, 2 * GMP_NUMB_BITS); 402 mpfr_set_prec (y, GMP_NUMB_BITS); 403 mpfr_set_ui (y, 1, MPFR_RNDN); 404 mpfr_nextabove (y); 405 mpfr_set (x, y, MPFR_RNDN); 406 inexact = test_sqrt (y, x, MPFR_RNDN); 407 if (mpfr_cmp_ui (y, 1) || inexact >= 0) 408 { 409 printf ("Error in sqrt(1) for prec=%d\n", (int) GMP_NUMB_BITS); 410 mpfr_dump (y); 411 exit (1); 412 } 413 414 mpfr_clear (x); 415 mpfr_clear (y); 416 mpfr_clear (z); 417} 418 419static void 420check_inexact (mpfr_prec_t p) 421{ 422 mpfr_t x, y, z; 423 mpfr_rnd_t rnd; 424 int inexact, sign; 425 426 mpfr_init2 (x, p); 427 mpfr_init2 (y, p); 428 mpfr_init2 (z, 2*p); 429 mpfr_urandomb (x, RANDS); 430 rnd = RND_RAND_NO_RNDF (); 431 inexact = test_sqrt (y, x, rnd); 432 if (mpfr_sqr (z, y, rnd)) /* exact since prec(z) = 2*prec(y) */ 433 { 434 printf ("Error: multiplication should be exact\n"); 435 exit (1); 436 } 437 mpfr_sub (z, z, x, rnd); /* exact also */ 438 sign = mpfr_cmp_ui (z, 0); 439 if (((inexact == 0) && (sign)) || 440 ((inexact > 0) && (sign <= 0)) || 441 ((inexact < 0) && (sign >= 0))) 442 { 443 printf ("Error with rnd=%s: wrong ternary value, expected %d, got %d\n", 444 mpfr_print_rnd_mode (rnd), sign, inexact); 445 printf ("x="); 446 mpfr_dump (x); 447 printf ("y="); 448 mpfr_dump (y); 449 exit (1); 450 } 451 mpfr_clear (x); 452 mpfr_clear (y); 453 mpfr_clear (z); 454} 455 456static void 457check_singular (void) 458{ 459 mpfr_t x, got; 460 461 mpfr_init2 (x, 100L); 462 mpfr_init2 (got, 100L); 463 464 /* sqrt(NaN) == NaN */ 465 MPFR_SET_NAN (x); 466 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 467 MPFR_ASSERTN (mpfr_nan_p (got)); 468 469 /* sqrt(-1) == NaN */ 470 mpfr_set_si (x, -1L, MPFR_RNDZ); 471 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 472 MPFR_ASSERTN (mpfr_nan_p (got)); 473 474 /* sqrt(+inf) == +inf */ 475 MPFR_SET_INF (x); 476 MPFR_SET_POS (x); 477 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 478 MPFR_ASSERTN (mpfr_inf_p (got)); 479 480 /* sqrt(-inf) == NaN */ 481 MPFR_SET_INF (x); 482 MPFR_SET_NEG (x); 483 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 484 MPFR_ASSERTN (mpfr_nan_p (got)); 485 486 /* sqrt(+0) == +0 */ 487 mpfr_set_si (x, 0L, MPFR_RNDZ); 488 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 489 MPFR_ASSERTN (mpfr_number_p (got)); 490 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); 491 MPFR_ASSERTN (MPFR_IS_POS (got)); 492 493 /* sqrt(-0) == -0 */ 494 mpfr_set_si (x, 0L, MPFR_RNDZ); 495 MPFR_SET_NEG (x); 496 MPFR_ASSERTN (test_sqrt (got, x, MPFR_RNDZ) == 0); /* exact */ 497 MPFR_ASSERTN (mpfr_number_p (got)); 498 MPFR_ASSERTN (mpfr_cmp_ui (got, 0L) == 0); 499 MPFR_ASSERTN (MPFR_IS_NEG (got)); 500 501 mpfr_clear (x); 502 mpfr_clear (got); 503} 504 505/* check that -1 <= x/sqrt(x^2+s*y^2) <= 1 for rounding to nearest or up 506 with s = 0 and s = 1 */ 507static void 508test_property1 (mpfr_prec_t p, mpfr_rnd_t r, int s) 509{ 510 mpfr_t x, y, z, t; 511 512 mpfr_init2 (x, p); 513 mpfr_init2 (y, p); 514 mpfr_init2 (z, p); 515 mpfr_init2 (t, p); 516 517 mpfr_urandomb (x, RANDS); 518 mpfr_mul (z, x, x, r); 519 if (s) 520 { 521 mpfr_urandomb (y, RANDS); 522 mpfr_mul (t, y, y, r); 523 mpfr_add (z, z, t, r); 524 } 525 mpfr_sqrt (z, z, r); 526 mpfr_div (z, x, z, r); 527 /* Note: if both x and y are 0, z is NAN, but the test below will 528 be false. So, everything is fine. */ 529 if (mpfr_cmp_si (z, -1) < 0 || mpfr_cmp_ui (z, 1) > 0) 530 { 531 printf ("Error, -1 <= x/sqrt(x^2+y^2) <= 1 does not hold for r=%s\n", 532 mpfr_print_rnd_mode (r)); 533 printf ("x="); mpfr_dump (x); 534 printf ("y="); mpfr_dump (y); 535 printf ("got "); mpfr_dump (z); 536 exit (1); 537 } 538 539 mpfr_clear (x); 540 mpfr_clear (y); 541 mpfr_clear (z); 542 mpfr_clear (t); 543} 544 545/* check sqrt(x^2) = x */ 546static void 547test_property2 (mpfr_prec_t p, mpfr_rnd_t r) 548{ 549 mpfr_t x, y; 550 551 mpfr_init2 (x, p); 552 mpfr_init2 (y, p); 553 554 mpfr_urandomb (x, RANDS); 555 mpfr_mul (y, x, x, r); 556 mpfr_sqrt (y, y, r); 557 if (mpfr_cmp (y, x)) 558 { 559 printf ("Error, sqrt(x^2) = x does not hold for r=%s\n", 560 mpfr_print_rnd_mode (r)); 561 printf ("x="); mpfr_dump (x); 562 printf ("got "); mpfr_dump (y); 563 exit (1); 564 } 565 566 mpfr_clear (x); 567 mpfr_clear (y); 568} 569 570/* Bug reported by Fredrik Johansson, occurring when: 571 - the precision of the result is a multiple of the number of bits 572 per word (GMP_NUMB_BITS), 573 - the rounding mode is to nearest (MPFR_RNDN), 574 - internally, the result has to be rounded up to a power of 2. 575*/ 576static void 577bug20160120 (void) 578{ 579 mpfr_t x, y; 580 581 mpfr_init2 (x, 4 * GMP_NUMB_BITS); 582 mpfr_init2 (y, GMP_NUMB_BITS); 583 584 mpfr_set_ui (x, 1, MPFR_RNDN); 585 mpfr_nextbelow (x); 586 mpfr_sqrt (y, x, MPFR_RNDN); 587 MPFR_ASSERTN(mpfr_check (y)); 588 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); 589 590 mpfr_set_prec (y, 2 * GMP_NUMB_BITS); 591 mpfr_sqrt (y, x, MPFR_RNDN); 592 MPFR_ASSERTN(mpfr_check (y)); 593 MPFR_ASSERTN(mpfr_cmp_ui (y, 1) == 0); 594 595 mpfr_clear(x); 596 mpfr_clear(y); 597} 598 599/* Bug in mpfr_sqrt2 when prec(u) = 2*GMP_NUMB_BITS and the exponent of u is 600 odd: the last bit of u is lost. */ 601static void 602bug20160908 (void) 603{ 604 mpfr_t r, u; 605 int n = GMP_NUMB_BITS, ret; 606 607 mpfr_init2 (r, 2*n - 1); 608 mpfr_init2 (u, 2 * n); 609 mpfr_set_ui_2exp (u, 1, 2*n-2, MPFR_RNDN); /* u=2^(2n-2) with exp(u)=2n-1 */ 610 mpfr_nextabove (u); 611 /* now u = 2^(2n-2) + 1/2 */ 612 ret = mpfr_sqrt (r, u, MPFR_RNDZ); 613 MPFR_ASSERTN(ret == -1 && mpfr_cmp_ui_2exp (r, 1, n-1) == 0); 614 mpfr_clear (r); 615 mpfr_clear (u); 616} 617 618static void 619testall_rndf (mpfr_prec_t pmax) 620{ 621 mpfr_t a, b, d; 622 mpfr_prec_t pa, pb; 623 624 for (pa = MPFR_PREC_MIN; pa <= pmax; pa++) 625 { 626 mpfr_init2 (a, pa); 627 mpfr_init2 (d, pa); 628 for (pb = MPFR_PREC_MIN; pb <= pmax; pb++) 629 { 630 mpfr_init2 (b, pb); 631 mpfr_set_ui (b, 1, MPFR_RNDN); 632 while (mpfr_cmp_ui (b, 4) < 0) 633 { 634 mpfr_sqrt (a, b, MPFR_RNDF); 635 mpfr_sqrt (d, b, MPFR_RNDD); 636 if (!mpfr_equal_p (a, d)) 637 { 638 mpfr_sqrt (d, b, MPFR_RNDU); 639 if (!mpfr_equal_p (a, d)) 640 { 641 printf ("Error: mpfr_sqrt(a,b,RNDF) does not " 642 "match RNDD/RNDU\n"); 643 printf ("b="); mpfr_dump (b); 644 printf ("a="); mpfr_dump (a); 645 exit (1); 646 } 647 } 648 mpfr_nextabove (b); 649 } 650 mpfr_clear (b); 651 } 652 mpfr_clear (a); 653 mpfr_clear (d); 654 } 655} 656 657/* test the case prec = GMP_NUMB_BITS */ 658static void 659test_sqrt1n (void) 660{ 661 mpfr_t r, u; 662 int inex; 663 664 MPFR_ASSERTD(GMP_NUMB_BITS >= 8); /* so that 15^2 is exactly representable */ 665 666 mpfr_init2 (r, GMP_NUMB_BITS); 667 mpfr_init2 (u, GMP_NUMB_BITS); 668 669 inex = mpfr_set_ui_2exp (u, 9 * 9, 2 * GMP_NUMB_BITS - 10, MPFR_RNDN); 670 MPFR_ASSERTN(inex == 0); 671 inex = mpfr_sqrt (r, u, MPFR_RNDN); 672 MPFR_ASSERTN(inex == 0); 673 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 9, GMP_NUMB_BITS - 5) == 0); 674 675 inex = mpfr_set_ui_2exp (u, 15 * 15, 2 * GMP_NUMB_BITS - 10, MPFR_RNDN); 676 MPFR_ASSERTN(inex == 0); 677 inex = mpfr_sqrt (r, u, MPFR_RNDN); 678 MPFR_ASSERTN(inex == 0); 679 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 15, GMP_NUMB_BITS - 5) == 0); 680 681 inex = mpfr_set_ui_2exp (u, 1, GMP_NUMB_BITS - 2, MPFR_RNDN); 682 MPFR_ASSERTN(inex == 0); 683 inex = mpfr_add_ui (u, u, 1, MPFR_RNDN); 684 MPFR_ASSERTN(inex == 0); 685 inex = mpfr_mul_2ui (u, u, GMP_NUMB_BITS, MPFR_RNDN); 686 MPFR_ASSERTN(inex == 0); 687 /* u = 2^(2*GMP_NUMB_BITS-2) + 2^GMP_NUMB_BITS, thus 688 u = r^2 + 2^GMP_NUMB_BITS with r = 2^(GMP_NUMB_BITS-1). 689 Should round to r+1 to nearest. */ 690 inex = mpfr_sqrt (r, u, MPFR_RNDN); 691 MPFR_ASSERTN(inex > 0); 692 inex = mpfr_sub_ui (r, r, 1, MPFR_RNDN); 693 MPFR_ASSERTN(inex == 0); 694 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, GMP_NUMB_BITS - 1) == 0); 695 696 mpfr_clear (r); 697 mpfr_clear (u); 698} 699 700static void 701check_overflow (void) 702{ 703 mpfr_t r, u; 704 mpfr_prec_t p; 705 mpfr_exp_t emax; 706 int inex; 707 708 emax = mpfr_get_emax (); 709 for (p = MPFR_PREC_MIN; p <= 1024; p++) 710 { 711 mpfr_init2 (r, p); 712 mpfr_init2 (u, p); 713 714 set_emax (-1); 715 mpfr_set_ui_2exp (u, 1, mpfr_get_emax () - 1, MPFR_RNDN); 716 mpfr_nextbelow (u); 717 mpfr_mul_2ui (u, u, 1, MPFR_RNDN); 718 /* now u = (1 - 2^(-p))*2^emax is the largest number < +Inf, 719 it square root is near 0.707 and has exponent 0 > emax */ 720 /* for RNDN, the result should be +Inf */ 721 inex = mpfr_sqrt (r, u, MPFR_RNDN); 722 MPFR_ASSERTN(inex > 0); 723 MPFR_ASSERTN(mpfr_inf_p (r) && mpfr_sgn (r) > 0); 724 /* for RNDA, the result should be +Inf */ 725 inex = mpfr_sqrt (r, u, MPFR_RNDA); 726 MPFR_ASSERTN(inex > 0); 727 MPFR_ASSERTN(mpfr_inf_p (r) && mpfr_sgn (r) > 0); 728 /* for RNDZ, the result should be u */ 729 inex = mpfr_sqrt (r, u, MPFR_RNDZ); 730 MPFR_ASSERTN(inex < 0); 731 MPFR_ASSERTN(mpfr_equal_p (r, u)); 732 733 set_emax (0); 734 mpfr_set_ui_2exp (u, 1, mpfr_get_emax () - 1, MPFR_RNDN); 735 mpfr_nextbelow (u); 736 mpfr_mul_2ui (u, u, 1, MPFR_RNDN); 737 /* u = 1-2^(-p), its square root is > u, and should thus give +Inf when 738 rounding away */ 739 inex = mpfr_sqrt (r, u, MPFR_RNDA); 740 MPFR_ASSERTN(inex > 0); 741 MPFR_ASSERTN(mpfr_inf_p (r) && mpfr_sgn (r) > 0); 742 743 mpfr_clear (r); 744 mpfr_clear (u); 745 } 746 set_emax (emax); 747} 748 749static void 750check_underflow (void) 751{ 752 mpfr_t r, u; 753 mpfr_prec_t p; 754 mpfr_exp_t emin; 755 int inex; 756 757 emin = mpfr_get_emin (); 758 for (p = MPFR_PREC_MIN; p <= 1024; p++) 759 { 760 mpfr_init2 (r, p); 761 mpfr_init2 (u, p); 762 763 set_emin (2); 764 mpfr_set_ui_2exp (u, 1, mpfr_get_emin () - 1, MPFR_RNDN); /* u = 2 */ 765 /* for RNDN, since sqrt(2) is closer from 2 than 0, the result is 2 */ 766 mpfr_clear_flags (); 767 inex = mpfr_sqrt (r, u, MPFR_RNDN); 768 MPFR_ASSERTN(inex > 0); 769 MPFR_ASSERTN(mpfr_equal_p (r, u)); 770 MPFR_ASSERTN(mpfr_underflow_p ()); 771 /* for RNDA, the result should be u, and there is underflow for p > 1, 772 since for p=1 we have 1 < sqrt(2) < 2, but for p >= 2, sqrt(2) should 773 be rounded to a number <= 1.5, which is representable */ 774 mpfr_clear_flags (); 775 inex = mpfr_sqrt (r, u, MPFR_RNDA); 776 MPFR_ASSERTN(inex > 0); 777 MPFR_ASSERTN(mpfr_equal_p (r, u)); 778 MPFR_ASSERTN((p == 1 && !mpfr_underflow_p ()) || 779 (p != 1 && mpfr_underflow_p ())); 780 /* for RNDZ, the result should be +0 */ 781 mpfr_clear_flags (); 782 inex = mpfr_sqrt (r, u, MPFR_RNDZ); 783 MPFR_ASSERTN(inex < 0); 784 MPFR_ASSERTN(mpfr_zero_p (r) && mpfr_signbit (r) == 0); 785 MPFR_ASSERTN(mpfr_underflow_p ()); 786 787 /* generate an input u such that sqrt(u) < 0.5*2^emin but there is no 788 underflow since sqrt(u) >= pred(0.5*2^emin), thus u >= 2^(2emin-2) */ 789 mpfr_set_ui_2exp (u, 1, 2 * mpfr_get_emin () - 2, MPFR_RNDN); 790 mpfr_clear_flags (); 791 inex = mpfr_sqrt (r, u, MPFR_RNDN); 792 MPFR_ASSERTN(inex == 0); 793 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 794 MPFR_ASSERTN(!mpfr_underflow_p ()); 795 mpfr_clear_flags (); 796 inex = mpfr_sqrt (r, u, MPFR_RNDA); 797 MPFR_ASSERTN(inex == 0); 798 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 799 MPFR_ASSERTN(!mpfr_underflow_p ()); 800 mpfr_clear_flags (); 801 inex = mpfr_sqrt (r, u, MPFR_RNDZ); 802 MPFR_ASSERTN(inex == 0); 803 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 804 MPFR_ASSERTN(!mpfr_underflow_p ()); 805 806 /* next number */ 807 mpfr_set_ui_2exp (u, 1, 2 * mpfr_get_emin () - 2, MPFR_RNDN); 808 mpfr_nextabove (u); 809 mpfr_clear_flags (); 810 inex = mpfr_sqrt (r, u, MPFR_RNDN); 811 MPFR_ASSERTN(inex < 0); 812 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 813 MPFR_ASSERTN(!mpfr_underflow_p ()); 814 mpfr_clear_flags (); 815 inex = mpfr_sqrt (r, u, MPFR_RNDA); 816 MPFR_ASSERTN(inex > 0); 817 mpfr_nextbelow (r); 818 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 819 MPFR_ASSERTN(!mpfr_underflow_p ()); 820 mpfr_clear_flags (); 821 inex = mpfr_sqrt (r, u, MPFR_RNDZ); 822 MPFR_ASSERTN(inex < 0); 823 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 824 MPFR_ASSERTN(!mpfr_underflow_p ()); 825 826 /* previous number */ 827 mpfr_set_ui_2exp (u, 1, 2 * mpfr_get_emin () - 2, MPFR_RNDN); 828 mpfr_nextbelow (u); 829 mpfr_clear_flags (); 830 inex = mpfr_sqrt (r, u, MPFR_RNDN); 831 MPFR_ASSERTN(inex > 0); 832 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 833 /* since sqrt(u) is just below the middle between 0.5*2^emin and 834 the previous number (with unbounded exponent range), there is 835 underflow */ 836 MPFR_ASSERTN(mpfr_underflow_p ()); 837 mpfr_clear_flags (); 838 inex = mpfr_sqrt (r, u, MPFR_RNDA); 839 MPFR_ASSERTN(inex > 0); 840 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 841 MPFR_ASSERTN(!mpfr_underflow_p ()); 842 mpfr_clear_flags (); 843 inex = mpfr_sqrt (r, u, MPFR_RNDZ); 844 MPFR_ASSERTN(inex < 0); 845 mpfr_nextabove (r); 846 MPFR_ASSERTN(mpfr_cmp_ui_2exp (r, 1, mpfr_get_emin () - 1) == 0); 847 MPFR_ASSERTN(mpfr_underflow_p ()); 848 849 set_emin (3); 850 mpfr_set_ui_2exp (u, 1, mpfr_get_emin () - 1, MPFR_RNDN); /* u = 4 */ 851 /* sqrt(u) = 2 = 0.5^2^(emin-1) should be rounded to +0 */ 852 mpfr_clear_flags (); 853 inex = mpfr_sqrt (r, u, MPFR_RNDN); 854 MPFR_ASSERTN(inex < 0); 855 MPFR_ASSERTN(mpfr_zero_p (r) && mpfr_signbit (r) == 0); 856 MPFR_ASSERTN(mpfr_underflow_p ()); 857 858 /* next number */ 859 mpfr_set_ui_2exp (u, 1, mpfr_get_emin () - 1, MPFR_RNDN); /* u = 4 */ 860 mpfr_nextabove (u); 861 /* sqrt(u) should be rounded to 4 */ 862 mpfr_clear_flags (); 863 inex = mpfr_sqrt (r, u, MPFR_RNDN); 864 MPFR_ASSERTN(inex > 0); 865 MPFR_ASSERTN(mpfr_cmp_ui (r, 4) == 0); 866 MPFR_ASSERTN(mpfr_underflow_p ()); 867 868 set_emin (4); 869 mpfr_set_ui_2exp (u, 1, mpfr_get_emin () - 1, MPFR_RNDN); /* u = 8 */ 870 /* sqrt(u) should be rounded to +0 */ 871 mpfr_clear_flags (); 872 inex = mpfr_sqrt (r, u, MPFR_RNDN); 873 MPFR_ASSERTN(inex < 0); 874 MPFR_ASSERTN(mpfr_zero_p (r) && mpfr_signbit (r) == 0); 875 MPFR_ASSERTN(mpfr_underflow_p ()); 876 877 mpfr_clear (r); 878 mpfr_clear (u); 879 } 880 set_emin (emin); 881} 882 883static void 884coverage (void) 885{ 886 mpfr_t r, t, u, v, w; 887 mpfr_prec_t p; 888 int inex; 889 890 /* exercise even rule */ 891 for (p = MPFR_PREC_MIN; p <= 1024; p++) 892 { 893 mpfr_init2 (r, p); 894 mpfr_init2 (t, p + 1); 895 mpfr_init2 (u, 2 * p + 2); 896 mpfr_init2 (v, p); 897 mpfr_init2 (w, p); 898 do 899 mpfr_urandomb (v, RANDS); 900 while (mpfr_zero_p (v)); 901 mpfr_set (w, v, MPFR_RNDN); 902 mpfr_nextabove (w); /* w = nextabove(v) */ 903 mpfr_set (t, v, MPFR_RNDN); 904 mpfr_nextabove (t); 905 mpfr_sqr (u, t, MPFR_RNDN); 906 inex = mpfr_sqrt (r, u, MPFR_RNDN); 907 if (mpfr_min_prec (v) < p) /* v is even */ 908 { 909 MPFR_ASSERTN(inex < 0); 910 MPFR_ASSERTN(mpfr_equal_p (r, v)); 911 } 912 else /* v is odd */ 913 { 914 MPFR_ASSERTN(inex > 0); 915 MPFR_ASSERTN(mpfr_equal_p (r, w)); 916 } 917 mpfr_clear (r); 918 mpfr_clear (t); 919 mpfr_clear (u); 920 mpfr_clear (v); 921 mpfr_clear (w); 922 } 923} 924 925#define TEST_FUNCTION test_sqrt 926#define TEST_RANDOM_POS 8 927#include "tgeneric.c" 928 929int 930main (void) 931{ 932 mpfr_prec_t p; 933 int k; 934 935 tests_start_mpfr (); 936 937 coverage (); 938 check_underflow (); 939 check_overflow (); 940 testall_rndf (16); 941 for (p = MPFR_PREC_MIN; p <= 128; p++) 942 { 943 test_property1 (p, MPFR_RNDN, 0); 944 test_property1 (p, MPFR_RNDU, 0); 945 test_property1 (p, MPFR_RNDA, 0); 946 test_property1 (p, MPFR_RNDN, 1); 947 test_property1 (p, MPFR_RNDU, 1); 948 test_property1 (p, MPFR_RNDA, 1); 949 test_property2 (p, MPFR_RNDN); 950 } 951 952 check_diverse ("635030154261163106768013773815762607450069292760790610550915652722277604820131530404842415587328", 160, "796887792767063979679855997149887366668464780637"); 953 special (); 954 check_singular (); 955 956 for (p=2; p<200; p++) 957 for (k=0; k<200; k++) 958 check_inexact (p); 959 check_float(); 960 961 check3 ("-0.0", MPFR_RNDN, "0.0"); 962 check4 ("6.37983013646045901440e+32", MPFR_RNDN, "5.9bc5036d09e0c@13"); 963 check4 ("1.0", MPFR_RNDN, "1"); 964 check4 ("1.0", MPFR_RNDZ, "1"); 965 check4 ("3.725290298461914062500000e-9", MPFR_RNDN, "4@-4"); 966 check4 ("3.725290298461914062500000e-9", MPFR_RNDZ, "4@-4"); 967 check4 ("1190456976439861.0", MPFR_RNDZ, "2.0e7957873529a@6"); 968 check4 ("1219027943874417664.0", MPFR_RNDZ, "4.1cf2af0e6a534@7"); 969 /* the following examples are bugs in Cygnus compiler/system, found by 970 Fabrice Rouillier while porting mpfr to Windows */ 971 check4 ("9.89438396044940256501e-134", MPFR_RNDU, "8.7af7bf0ebbee@-56"); 972 check4 ("7.86528588050363751914e+31", MPFR_RNDZ, "1.f81fc40f32062@13"); 973 check4 ("0.99999999999999988897", MPFR_RNDN, "f.ffffffffffff8@-1"); 974 check4 ("1.00000000000000022204", MPFR_RNDN, "1"); 975 /* the following examples come from the paper "Number-theoretic Test 976 Generation for Directed Rounding" from Michael Parks, Table 4 */ 977 978 check4 ("78652858805036375191418371571712.0", MPFR_RNDN, 979 "1.f81fc40f32063@13"); 980 check4 ("38510074998589467860312736661504.0", MPFR_RNDN, 981 "1.60c012a92fc65@13"); 982 check4 ("35318779685413012908190921129984.0", MPFR_RNDN, 983 "1.51d17526c7161@13"); 984 check4 ("26729022595358440976973142425600.0", MPFR_RNDN, 985 "1.25e19302f7e51@13"); 986 check4 ("22696567866564242819241453027328.0", MPFR_RNDN, 987 "1.0ecea7dd2ec3d@13"); 988 check4 ("22696888073761729132924856434688.0", MPFR_RNDN, 989 "1.0ecf250e8e921@13"); 990 check4 ("36055652513981905145251657416704.0", MPFR_RNDN, 991 "1.5552f3eedcf33@13"); 992 check4 ("30189856268896404997497182748672.0", MPFR_RNDN, 993 "1.3853ee10c9c99@13"); 994 check4 ("36075288240584711210898775080960.0", MPFR_RNDN, 995 "1.556abe212b56f@13"); 996 check4 ("72154663483843080704304789585920.0", MPFR_RNDN, 997 "1.e2d9a51977e6e@13"); 998 999 check4 ("78652858805036375191418371571712.0", MPFR_RNDZ, 1000 "1.f81fc40f32062@13"); 1001 check4 ("38510074998589467860312736661504.0", MPFR_RNDZ, 1002 "1.60c012a92fc64@13"); 1003 check4 ("35318779685413012908190921129984.0", MPFR_RNDZ, "1.51d17526c716@13"); 1004 check4 ("26729022595358440976973142425600.0", MPFR_RNDZ, "1.25e19302f7e5@13"); 1005 check4 ("22696567866564242819241453027328.0", MPFR_RNDZ, 1006 "1.0ecea7dd2ec3c@13"); 1007 check4 ("22696888073761729132924856434688.0", MPFR_RNDZ, "1.0ecf250e8e92@13"); 1008 check4 ("36055652513981905145251657416704.0", MPFR_RNDZ, 1009 "1.5552f3eedcf32@13"); 1010 check4 ("30189856268896404997497182748672.0", MPFR_RNDZ, 1011 "1.3853ee10c9c98@13"); 1012 check4 ("36075288240584711210898775080960.0", MPFR_RNDZ, 1013 "1.556abe212b56e@13"); 1014 check4 ("72154663483843080704304789585920.0", MPFR_RNDZ, 1015 "1.e2d9a51977e6d@13"); 1016 1017 check4 ("78652858805036375191418371571712.0", MPFR_RNDU, 1018 "1.f81fc40f32063@13"); 1019 check4 ("38510074998589467860312736661504.0", MPFR_RNDU, 1020 "1.60c012a92fc65@13"); 1021 check4 ("35318779685413012908190921129984.0", MPFR_RNDU, 1022 "1.51d17526c7161@13"); 1023 check4 ("26729022595358440976973142425600.0", MPFR_RNDU, 1024 "1.25e19302f7e51@13"); 1025 check4 ("22696567866564242819241453027328.0", MPFR_RNDU, 1026 "1.0ecea7dd2ec3d@13"); 1027 check4 ("22696888073761729132924856434688.0", MPFR_RNDU, 1028 "1.0ecf250e8e921@13"); 1029 check4 ("36055652513981905145251657416704.0", MPFR_RNDU, 1030 "1.5552f3eedcf33@13"); 1031 check4 ("30189856268896404997497182748672.0", MPFR_RNDU, 1032 "1.3853ee10c9c99@13"); 1033 check4 ("36075288240584711210898775080960.0", MPFR_RNDU, 1034 "1.556abe212b56f@13"); 1035 check4 ("72154663483843080704304789585920.0", MPFR_RNDU, 1036 "1.e2d9a51977e6e@13"); 1037 1038 check4 ("78652858805036375191418371571712.0", MPFR_RNDD, 1039 "1.f81fc40f32062@13"); 1040 check4 ("38510074998589467860312736661504.0", MPFR_RNDD, 1041 "1.60c012a92fc64@13"); 1042 check4 ("35318779685413012908190921129984.0", MPFR_RNDD, "1.51d17526c716@13"); 1043 check4 ("26729022595358440976973142425600.0", MPFR_RNDD, "1.25e19302f7e5@13"); 1044 check4 ("22696567866564242819241453027328.0", MPFR_RNDD, 1045 "1.0ecea7dd2ec3c@13"); 1046 check4 ("22696888073761729132924856434688.0", MPFR_RNDD, "1.0ecf250e8e92@13"); 1047 check4 ("36055652513981905145251657416704.0", MPFR_RNDD, 1048 "1.5552f3eedcf32@13"); 1049 check4 ("30189856268896404997497182748672.0", MPFR_RNDD, 1050 "1.3853ee10c9c98@13"); 1051 check4 ("36075288240584711210898775080960.0", MPFR_RNDD, 1052 "1.556abe212b56e@13"); 1053 check4 ("72154663483843080704304789585920.0", MPFR_RNDD, 1054 "1.e2d9a51977e6d@13"); 1055 1056 /* check that rounding away is just rounding toward positive infinity */ 1057 check4 ("72154663483843080704304789585920.0", MPFR_RNDA, 1058 "1.e2d9a51977e6e@13"); 1059 1060 test_generic (MPFR_PREC_MIN, 300, 15); 1061 data_check ("data/sqrt", mpfr_sqrt, "mpfr_sqrt"); 1062 bad_cases (mpfr_sqrt, mpfr_sqr, "mpfr_sqrt", 0, -256, 255, 4, 128, 80, 50); 1063 1064 bug20160120 (); 1065 bug20160908 (); 1066 test_sqrt1n (); 1067 1068 tests_end_mpfr (); 1069 return 0; 1070} 1071