1 /* Test file for mpfr_sqr. 2 3 Copyright 2004-2020 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 #define TEST_FUNCTION mpfr_sqr 26 #include "tgeneric.c" 27 28 static int 29 inexact_sign (int x) 30 { 31 return (x < 0) ? -1 : (x > 0); 32 } 33 34 static void 35 error1 (mpfr_rnd_t rnd, mpfr_prec_t prec, 36 mpfr_t in, mpfr_t outmul, mpfr_t outsqr) 37 { 38 printf("ERROR: for %s and prec=%lu\nINPUT=", mpfr_print_rnd_mode(rnd), 39 (unsigned long) prec); 40 mpfr_dump(in); 41 printf("OutputMul="); mpfr_dump(outmul); 42 printf("OutputSqr="); mpfr_dump(outsqr); 43 exit(1); 44 } 45 46 static void 47 error2 (mpfr_rnd_t rnd, mpfr_prec_t prec, mpfr_t in, mpfr_t out, 48 int inexactmul, int inexactsqr) 49 { 50 printf("ERROR: for %s and prec=%lu\nINPUT=", mpfr_print_rnd_mode(rnd), 51 (unsigned long) prec); 52 mpfr_dump(in); 53 printf("Output="); mpfr_dump(out); 54 printf("InexactMul= %d InexactSqr= %d\n", inexactmul, inexactsqr); 55 exit(1); 56 } 57 58 static void 59 check_random (mpfr_prec_t p) 60 { 61 mpfr_t x,y,z; 62 int r; 63 int i, inexact1, inexact2; 64 65 mpfr_inits2 (p, x, y, z, (mpfr_ptr) 0); 66 for(i = 0 ; i < 500 ; i++) 67 { 68 mpfr_urandomb (x, RANDS); 69 if (MPFR_IS_PURE_FP(x)) 70 RND_LOOP_NO_RNDF (r) 71 { 72 inexact1 = mpfr_mul (y, x, x, (mpfr_rnd_t) r); 73 inexact2 = mpfr_sqr (z, x, (mpfr_rnd_t) r); 74 if (mpfr_cmp (y, z)) 75 error1 ((mpfr_rnd_t) r,p,x,y,z); 76 if (inexact_sign (inexact1) != inexact_sign (inexact2)) 77 error2 ((mpfr_rnd_t) r,p,x,y,inexact1,inexact2); 78 } 79 } 80 mpfr_clears (x, y, z, (mpfr_ptr) 0); 81 } 82 83 static void 84 check_special (void) 85 { 86 mpfr_t x, y; 87 mpfr_exp_t emin; 88 89 mpfr_init (x); 90 mpfr_init (y); 91 92 mpfr_set_nan (x); 93 mpfr_sqr (y, x, MPFR_RNDN); 94 MPFR_ASSERTN (mpfr_nan_p (y)); 95 96 mpfr_set_inf (x, 1); 97 mpfr_sqr (y, x, MPFR_RNDN); 98 MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0); 99 100 mpfr_set_inf (x, -1); 101 mpfr_sqr (y, x, MPFR_RNDN); 102 MPFR_ASSERTN (mpfr_inf_p (y) && mpfr_sgn (y) > 0); 103 104 mpfr_set_ui (x, 0, MPFR_RNDN); 105 mpfr_sqr (y, x, MPFR_RNDN); 106 MPFR_ASSERTN (mpfr_zero_p (y)); 107 108 emin = mpfr_get_emin (); 109 mpfr_set_emin (0); 110 mpfr_set_ui (x, 1, MPFR_RNDN); 111 mpfr_div_2ui (x, x, 1, MPFR_RNDN); 112 MPFR_ASSERTN (!mpfr_zero_p (x)); 113 mpfr_sqr (y, x, MPFR_RNDN); 114 MPFR_ASSERTN (mpfr_zero_p (y)); 115 mpfr_set_emin (emin); 116 117 mpfr_clear (y); 118 mpfr_clear (x); 119 } 120 121 /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */ 122 static void 123 test_underflow (void) 124 { 125 mpfr_t x, y; 126 mpfr_exp_t emin; 127 128 emin = mpfr_get_emin (); 129 mpfr_set_emin (0); 130 131 mpfr_init2 (x, 24); 132 mpfr_init2 (y, 24); 133 134 mpfr_set_ui_2exp (x, 11863283, -24, MPFR_RNDN); 135 /* x^2 = 0.011111111111111111111111101101100111111010101001*2^(-48) 136 thus we have an underflow */ 137 mpfr_clear_underflow (); 138 mpfr_sqr (y, x, MPFR_RNDN); 139 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0); 140 MPFR_ASSERTN(mpfr_underflow_p ()); 141 142 mpfr_set_prec (x, 69); 143 mpfr_set_prec (y, 69); 144 mpfr_set_str_binary (x, "0.101101010000010011110011001100111111100111011110011001001000010001011"); 145 mpfr_clear_underflow (); 146 mpfr_sqr (y, x, MPFR_RNDN); 147 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 1, -1) == 0); 148 MPFR_ASSERTN(mpfr_underflow_p ()); 149 150 mpfr_clear (y); 151 mpfr_clear (x); 152 153 mpfr_set_emin (emin); 154 } 155 156 /* Test of a bug seen with GCC 4.5.2 and GMP 5.0.1 on m68k (m68000 target). 157 https://sympa.inria.fr/sympa/arc/mpfr/2011-05/msg00003.html 158 https://sympa.inria.fr/sympa/arc/mpfr/2011-05/msg00041.html 159 */ 160 static void 161 check_mpn_sqr (void) 162 { 163 #if GMP_NUMB_BITS == 32 && __GNU_MP_VERSION >= 5 164 /* Note: since we test a low-level bug, src is initialized 165 without using a GMP function, just in case. */ 166 mp_limb_t src[5] = 167 { 0x90000000, 0xbaa55f4f, 0x2cbec4d9, 0xfcef3242, 0xda827999 }; 168 mp_limb_t exd[10] = 169 { 0x00000000, 0x31000000, 0xbd4bc59a, 0x41fbe2b5, 0x33471e7e, 170 0x90e826a7, 0xbaa55f4f, 0x2cbec4d9, 0xfcef3242, 0xba827999 }; 171 mp_limb_t dst[10]; 172 int i; 173 174 mpn_sqr (dst, src, 5); /* new in GMP 5 */ 175 for (i = 0; i < 10; i++) 176 { 177 if (dst[i] != exd[i]) 178 { 179 printf ("Error in check_mpn_sqr\n"); 180 printf ("exd[%d] = 0x%08lx\n", i, (unsigned long) exd[i]); 181 printf ("dst[%d] = 0x%08lx\n", i, (unsigned long) dst[i]); 182 printf ("Note: This is not a bug in MPFR, but a bug in" 183 " either GMP or, more\nprobably, in the compiler." 184 " It may cause other tests to fail.\n"); 185 exit (1); 186 } 187 } 188 #endif 189 } 190 191 static void 192 coverage (mpfr_prec_t pmax) 193 { 194 mpfr_prec_t p; 195 196 for (p = MPFR_PREC_MIN; p <= pmax; p++) 197 { 198 mpfr_t a, b, c; 199 int inex; 200 mpfr_exp_t emin; 201 202 mpfr_init2 (a, p); 203 mpfr_init2 (b, p); 204 mpfr_init2 (c, p); 205 206 /* exercise carry in most significant bits of a, with overflow */ 207 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); 208 mpfr_sqrt (b, b, MPFR_RNDU); 209 mpfr_div_2ui (c, b, 1, MPFR_RNDN); 210 mpfr_sqr (c, c, MPFR_RNDN); 211 mpfr_clear_flags (); 212 inex = mpfr_sqr (a, b, MPFR_RNDN); 213 /* if EXP(c) > emax-2, there is overflow */ 214 if (mpfr_get_exp (c) > mpfr_get_emax () - 2) 215 { 216 MPFR_ASSERTN(inex > 0); 217 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); 218 MPFR_ASSERTN(mpfr_overflow_p ()); 219 } 220 else /* no overflow */ 221 { 222 /* 2^p-1 is a square only for p=1 */ 223 MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0)); 224 MPFR_ASSERTN(!mpfr_overflow_p ()); 225 mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ); 226 MPFR_ASSERTN(mpfr_equal_p (a, c)); 227 } 228 229 /* same as above, with RNDU */ 230 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); 231 mpfr_sqrt (b, b, MPFR_RNDU); 232 mpfr_div_2ui (c, b, 1, MPFR_RNDN); 233 mpfr_sqr (c, c, MPFR_RNDU); 234 mpfr_clear_flags (); 235 inex = mpfr_sqr (a, b, MPFR_RNDU); 236 /* if EXP(c) > emax-2, there is overflow */ 237 if (mpfr_get_exp (c) > mpfr_get_emax () - 2) 238 { 239 MPFR_ASSERTN(inex > 0); 240 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); 241 MPFR_ASSERTN(mpfr_overflow_p ()); 242 } 243 else /* no overflow */ 244 { 245 /* 2^p-1 is a square only for p=1 */ 246 MPFR_ASSERTN((p == 1 && inex == 0) || (p > 1 && inex < 0)); 247 MPFR_ASSERTN(!mpfr_overflow_p ()); 248 mpfr_set_ui_2exp (c, 1, mpfr_get_emax (), MPFR_RNDZ); 249 MPFR_ASSERTN(mpfr_equal_p (a, c)); 250 } 251 252 /* exercise trivial overflow */ 253 mpfr_set_ui_2exp (b, 1, mpfr_get_emax (), MPFR_RNDZ); 254 mpfr_sqrt (b, b, MPFR_RNDU); 255 mpfr_mul_2ui (b, b, 1, MPFR_RNDN); 256 mpfr_clear_flags (); 257 inex = mpfr_sqr (a, b, MPFR_RNDN); 258 MPFR_ASSERTN(inex > 0); 259 MPFR_ASSERTN(mpfr_inf_p (a) && mpfr_sgn (a) > 0); 260 MPFR_ASSERTN(mpfr_overflow_p ()); 261 262 /* exercise trivial underflow */ 263 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDZ); 264 mpfr_sqrt (b, b, MPFR_RNDU); 265 mpfr_div_2ui (b, b, 1, MPFR_RNDN); 266 mpfr_clear_flags (); 267 inex = mpfr_sqr (a, b, MPFR_RNDN); 268 MPFR_ASSERTN(inex < 0); 269 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); 270 MPFR_ASSERTN(mpfr_underflow_p ()); 271 272 /* exercise square between 0.5*2^emin and its predecessor (emin even) */ 273 emin = mpfr_get_emin (); 274 mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ 275 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN); 276 inex = mpfr_sqrt (b, b, MPFR_RNDZ); 277 MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */ 278 mpfr_mul_2ui (c, b, 1, MPFR_RNDN); 279 mpfr_sqr (c, c, MPFR_RNDN); 280 mpfr_clear_flags (); 281 inex = mpfr_sqr (a, b, MPFR_RNDN); 282 if (mpfr_get_exp (c) < mpfr_get_emin () + 2) /* underflow */ 283 { 284 /* if c > 0.5*2^(emin+1), we should round to 0.5*2^emin */ 285 if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin ()) > 0) 286 { 287 MPFR_ASSERTN(inex > 0); 288 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); 289 MPFR_ASSERTN(mpfr_underflow_p ()); 290 } 291 else /* we should round to 0 */ 292 { 293 MPFR_ASSERTN(inex < 0); 294 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); 295 MPFR_ASSERTN(mpfr_underflow_p ()); 296 } 297 } 298 else 299 { 300 MPFR_ASSERTN(inex > 0); 301 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); 302 MPFR_ASSERTN(!mpfr_underflow_p ()); 303 } 304 mpfr_set_emin (emin); 305 306 /* exercise exact square root 2^(emin-2) for emin even */ 307 emin = mpfr_get_emin (); 308 mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ 309 mpfr_set_ui_2exp (b, 1, (mpfr_get_emin () - 2) / 2, MPFR_RNDN); 310 inex = mpfr_sqr (a, b, MPFR_RNDN); 311 MPFR_ASSERTN(inex < 0); 312 MPFR_ASSERTN(mpfr_zero_p (a) && mpfr_signbit (a) == 0); 313 MPFR_ASSERTN(mpfr_underflow_p ()); 314 mpfr_set_emin (emin); 315 316 /* same as above, for RNDU */ 317 emin = mpfr_get_emin (); 318 mpfr_set_emin (emin + (emin & 1)); /* now emin is even */ 319 mpfr_set_ui_2exp (b, 1, mpfr_get_emin () - 1, MPFR_RNDN); 320 inex = mpfr_sqrt (b, b, MPFR_RNDZ); 321 MPFR_ASSERTN(inex != 0); /* sqrt(2) is not exact */ 322 mpfr_mul_2ui (c, b, 1, MPFR_RNDN); 323 mpfr_sqr (c, c, MPFR_RNDU); 324 mpfr_clear_flags (); 325 inex = mpfr_sqr (a, b, MPFR_RNDU); 326 MPFR_ASSERTN(inex > 0); 327 MPFR_ASSERTN(mpfr_cmp_ui_2exp (a, 1, mpfr_get_emin () - 1) == 0); 328 /* we have underflow if c < 2^(emin+1) */ 329 if (mpfr_cmp_ui_2exp (c, 1, mpfr_get_emin () + 1) < 0) 330 MPFR_ASSERTN(mpfr_underflow_p ()); 331 else 332 MPFR_ASSERTN(!mpfr_underflow_p ()); 333 mpfr_set_emin (emin); 334 335 mpfr_clear (a); 336 mpfr_clear (b); 337 mpfr_clear (c); 338 } 339 } 340 341 int 342 main (void) 343 { 344 mpfr_prec_t p; 345 346 tests_start_mpfr (); 347 348 coverage (1024); 349 check_mpn_sqr (); 350 check_special (); 351 test_underflow (); 352 353 for (p = MPFR_PREC_MIN; p < 200; p++) 354 check_random (p); 355 356 test_generic (MPFR_PREC_MIN, 200, 15); 357 data_check ("data/sqr", mpfr_sqr, "mpfr_sqr"); 358 bad_cases (mpfr_sqr, mpfr_sqrt, "mpfr_sqr", 8, -256, 255, 4, 128, 800, 50); 359 360 tests_end_mpfr (); 361 return 0; 362 } 363