1 /* Test file for mpfr_sqrt.
2
3 Copyright 1999, 2001-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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 "mpfr-test.h"
24
25 #ifdef CHECK_EXTERNAL
26 static int
test_sqrt(mpfr_ptr a,mpfr_srcptr b,mpfr_rnd_t rnd_mode)27 test_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
48 static void
check3(const char * as,mpfr_rnd_t rnd_mode,const char * qs)49 check3 (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
68 static void
check4(const char * as,mpfr_rnd_t rnd_mode,const char * qs)69 check4 (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
89 static void
check24(const char * as,mpfr_rnd_t rnd_mode,const char * qs)90 check24 (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
109 static void
check_diverse(const char * as,mpfr_prec_t p,const char * qs)110 check_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 */
131 static void
check_float(void)132 check_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
182 static void
special(void)183 special (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
419 static void
check_inexact(mpfr_prec_t p)420 check_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
456 static void
check_singular(void)457 check_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 */
507 static void
test_property1(mpfr_prec_t p,mpfr_rnd_t r,int s)508 test_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 */
546 static void
test_property2(mpfr_prec_t p,mpfr_rnd_t r)547 test_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 */
576 static void
bug20160120(void)577 bug20160120 (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. */
601 static void
bug20160908(void)602 bug20160908 (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
618 static void
testall_rndf(mpfr_prec_t pmax)619 testall_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 */
658 static void
test_sqrt1n(void)659 test_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
700 static void
check_overflow(void)701 check_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
749 static void
check_underflow(void)750 check_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
883 static void
coverage(void)884 coverage (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
929 int
main(void)930 main (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