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