1 /* Test file for mpfr_set_ld and mpfr_get_ld. 2 3 Copyright 2002-2016 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 http://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 <stdio.h> 24 #include <stdlib.h> 25 #include <float.h> 26 #include <limits.h> 27 #ifdef WITH_FPU_CONTROL 28 #include <fpu_control.h> 29 #endif 30 31 #include "mpfr-test.h" 32 33 static void 34 check_gcc33_bug (void) 35 { 36 volatile long double x; 37 x = (long double) 9007199254740992.0 + 1.0; 38 if (x != 0.0) 39 return; /* OK */ 40 printf 41 ("Detected optimization bug of gcc 3.3 on Alpha concerning long double\n" 42 "comparisons; set_ld tests might fail (set_ld won't work correctly).\n" 43 "See http://gcc.gnu.org/ml/gcc-bugs/2003-10/msg00853.html for more\n" 44 "information.\n"); 45 } 46 47 static int 48 Isnan_ld (long double d) 49 { 50 /* Do not convert d to double as this can give an overflow, which 51 may confuse compilers without IEEE 754 support (such as clang 52 -fsanitize=undefined), or trigger a trap if enabled. 53 The DOUBLE_ISNAN macro should work fine on long double. */ 54 if (DOUBLE_ISNAN (d)) 55 return 1; 56 LONGDOUBLE_NAN_ACTION (d, goto yes); 57 return 0; 58 yes: 59 return 1; 60 } 61 62 /* checks that a long double converted to a mpfr (with precision >=113), 63 then converted back to a long double gives the initial value, 64 or in other words mpfr_get_ld(mpfr_set_ld(d)) = d. 65 */ 66 static void 67 check_set_get (long double d, mpfr_t x) 68 { 69 int r; 70 long double e; 71 int inex; 72 73 for (r = 0; r < MPFR_RND_MAX; r++) 74 { 75 inex = mpfr_set_ld (x, d, (mpfr_rnd_t) r); 76 if (inex != 0) 77 { 78 mpfr_exp_t emin, emax; 79 emin = mpfr_get_emin (); 80 emax = mpfr_get_emax (); 81 printf ("Error: mpfr_set_ld should be exact\n"); 82 printf ("d=%1.30Le inex=%d\n", d, inex); 83 if (emin >= LONG_MIN) 84 printf ("emin=%ld\n", (long) emin); 85 if (emax <= LONG_MAX) 86 printf ("emax=%ld\n", (long) emax); 87 mpfr_dump (x); 88 exit (1); 89 } 90 e = mpfr_get_ld (x, (mpfr_rnd_t) r); 91 if ((Isnan_ld(d) && ! Isnan_ld(e)) || 92 (Isnan_ld(e) && ! Isnan_ld(d)) || 93 (e != d && !(Isnan_ld(e) && Isnan_ld(d)))) 94 { 95 printf ("Error: mpfr_get_ld o mpfr_set_ld <> Id\n"); 96 printf (" r=%d\n", r); 97 printf (" d=%1.30Le get_ld(set_ld(d))=%1.30Le\n", d, e); 98 ld_trace (" d", d); 99 printf (" x="); mpfr_out_str (NULL, 16, 0, x, MPFR_RNDN); 100 printf ("\n"); 101 ld_trace (" e", e); 102 #ifdef MPFR_NANISNAN 103 if (Isnan_ld(d) || Isnan_ld(e)) 104 printf ("The reason is that NAN == NAN. Please look at the " 105 "configure output\nand Section \"In case of problem\" " 106 "of the INSTALL file.\n"); 107 #endif 108 exit (1); 109 } 110 } 111 } 112 113 static void 114 test_small (void) 115 { 116 mpfr_t x, y, z; 117 long double d; 118 119 mpfr_init2 (x, 64); 120 mpfr_init2 (y, 64); 121 mpfr_init2 (z, 64); 122 123 /* x = 11906603631607553907/2^(16381+64) */ 124 mpfr_set_str (x, "0.1010010100111100110000001110101101000111010110000001111101110011E-16381", 2, MPFR_RNDN); 125 d = mpfr_get_ld (x, MPFR_RNDN); /* infinite loop? */ 126 mpfr_set_ld (y, d, MPFR_RNDN); 127 mpfr_sub (z, x, y, MPFR_RNDN); 128 mpfr_abs (z, z, MPFR_RNDN); 129 mpfr_clear_erangeflag (); 130 /* If long double = double, d should be equal to 0; 131 in this case, everything is OK. */ 132 if (d != 0 && (mpfr_cmp_str (z, "1E-16434", 2, MPFR_RNDN) > 0 || 133 mpfr_erangeflag_p ())) 134 { 135 printf ("Error with x = "); 136 mpfr_out_str (NULL, 10, 21, x, MPFR_RNDN); 137 printf (" = "); 138 mpfr_out_str (NULL, 16, 0, x, MPFR_RNDN); 139 printf ("\n -> d = %.21Lg", d); 140 printf ("\n -> y = "); 141 mpfr_out_str (NULL, 10, 21, y, MPFR_RNDN); 142 printf (" = "); 143 mpfr_out_str (NULL, 16, 0, y, MPFR_RNDN); 144 printf ("\n -> |x-y| = "); 145 mpfr_out_str (NULL, 16, 0, z, MPFR_RNDN); 146 printf ("\n"); 147 exit (1); 148 } 149 150 mpfr_clear (x); 151 mpfr_clear (y); 152 mpfr_clear (z); 153 } 154 155 static void 156 test_fixed_bugs (void) 157 { 158 mpfr_t x; 159 long double l, m; 160 161 /* bug found by Steve Kargl (2009-03-14) */ 162 mpfr_init2 (x, 64); 163 mpfr_set_ui_2exp (x, 1, -16447, MPFR_RNDN); 164 mpfr_get_ld (x, MPFR_RNDN); /* an assertion failed in init2.c:50 */ 165 166 /* bug reported by Jakub Jelinek (2010-10-17) 167 https://gforge.inria.fr/tracker/?func=detail&aid=11300 */ 168 mpfr_set_prec (x, MPFR_LDBL_MANT_DIG); 169 /* l = 0x1.23456789abcdef0123456789abcdp-914L; */ 170 l = 8.215640181713713164092636634579e-276; 171 mpfr_set_ld (x, l, MPFR_RNDN); 172 m = mpfr_get_ld (x, MPFR_RNDN); 173 if (m != l) 174 { 175 printf ("Error in get_ld o set_ld for l=%Le\n", l); 176 printf ("Got m=%Le instead of l\n", m); 177 exit (1); 178 } 179 180 /* another similar test which failed with extended double precision and the 181 generic code for mpfr_set_ld */ 182 /* l = 0x1.23456789abcdef0123456789abcdp-968L; */ 183 l = 4.560596445887084662336528403703e-292; 184 mpfr_set_ld (x, l, MPFR_RNDN); 185 m = mpfr_get_ld (x, MPFR_RNDN); 186 if (m != l) 187 { 188 printf ("Error in get_ld o set_ld for l=%Le\n", l); 189 printf ("Got m=%Le instead of l\n", m); 190 exit (1); 191 } 192 193 mpfr_clear (x); 194 } 195 196 /* bug reported by Walter Mascarenhas 197 https://sympa.inria.fr/sympa/arc/mpfr/2016-09/msg00005.html */ 198 static void 199 bug_20160907 (void) 200 { 201 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE 202 long double dn, ld; 203 mpfr_t mp; 204 long e; 205 mpfr_long_double_t x; 206 207 /* the following is the encoding of the smallest subnormal number 208 for HAVE_LDOUBLE_IEEE_EXT_LITTLE */ 209 x.s.manl = 1; 210 x.s.manh = 0; 211 x.s.expl = 0; 212 x.s.exph = 0; 213 x.s.sign= 0; 214 dn = x.ld; 215 e = -16445; 216 /* dn=2^e is now the smallest subnormal. */ 217 218 mpfr_init2 (mp, 64); 219 mpfr_set_ui_2exp (mp, 1, e - 1, MPFR_RNDN); 220 ld = mpfr_get_ld (mp, MPFR_RNDU); 221 /* since mp = 2^(e-1) and ld is rounded upwards, we should have 222 ld = 2^e */ 223 if (ld != dn) 224 { 225 printf ("Error, ld = %Le <> dn = %Le\n", ld, dn); 226 printf ("mp="); 227 mpfr_out_str (stdout, 10, 0, mp, MPFR_RNDN); 228 printf ("\n"); 229 exit (1); 230 } 231 232 /* check a few more numbers */ 233 for (e = -16446; e <= -16381; e++) 234 { 235 mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); 236 ld = mpfr_get_ld (mp, MPFR_RNDU); 237 mpfr_set_ld (mp, ld, MPFR_RNDU); 238 /* mp is 2^e rounded up, thus should be >= 2^e */ 239 MPFR_ASSERTN(mpfr_cmp_ui_2exp (mp, 1, e) >= 0); 240 241 mpfr_set_ui_2exp (mp, 1, e, MPFR_RNDN); 242 ld = mpfr_get_ld (mp, MPFR_RNDD); 243 mpfr_set_ld (mp, ld, MPFR_RNDD); 244 /* mp is 2^e rounded down, thus should be <= 2^e */ 245 if (mpfr_cmp_ui_2exp (mp, 3, e) > 0) 246 { 247 printf ("Error, expected value <= 2^%ld\n", e); 248 printf ("got "); mpfr_dump (mp); 249 exit (1); 250 } 251 } 252 253 mpfr_clear (mp); 254 #endif 255 } 256 257 int 258 main (int argc, char *argv[]) 259 { 260 long double d, e; 261 mpfr_t x; 262 int i; 263 mpfr_exp_t emax; 264 #ifdef WITH_FPU_CONTROL 265 fpu_control_t cw; 266 267 if (argc > 1) 268 { 269 cw = strtol(argv[1], NULL, 0); 270 printf ("FPU control word: 0x%x\n", (unsigned int) cw); 271 _FPU_SETCW (cw); 272 } 273 #endif 274 275 tests_start_mpfr (); 276 mpfr_test_init (); 277 278 check_gcc33_bug (); 279 test_fixed_bugs (); 280 281 mpfr_init2 (x, MPFR_LDBL_MANT_DIG); 282 283 #if !defined(MPFR_ERRDIVZERO) 284 /* check NaN */ 285 mpfr_set_nan (x); 286 d = mpfr_get_ld (x, MPFR_RNDN); 287 check_set_get (d, x); 288 #endif 289 290 /* check +0.0 and -0.0 */ 291 d = 0.0; 292 check_set_get (d, x); 293 d = DBL_NEG_ZERO; 294 check_set_get (d, x); 295 296 /* check that the sign of -0.0 is set */ 297 mpfr_set_ld (x, DBL_NEG_ZERO, MPFR_RNDN); 298 if (MPFR_SIGN(x) > 0) 299 { 300 printf ("Error: sign of -0.0 is not set correctly\n"); 301 #if _GMP_IEEE_FLOATS 302 exit (1); 303 /* Non IEEE doesn't support negative zero yet */ 304 #endif 305 } 306 307 #if !defined(MPFR_ERRDIVZERO) 308 /* check +Inf */ 309 mpfr_set_inf (x, 1); 310 d = mpfr_get_ld (x, MPFR_RNDN); 311 check_set_get (d, x); 312 313 /* check -Inf */ 314 mpfr_set_inf (x, -1); 315 d = mpfr_get_ld (x, MPFR_RNDN); 316 check_set_get (d, x); 317 #endif 318 319 /* check the largest power of two */ 320 d = 1.0; while (d < LDBL_MAX / 2.0) d += d; 321 check_set_get (d, x); 322 check_set_get (-d, x); 323 324 /* check largest long double */ 325 d = LDBL_MAX; 326 check_set_get (d, x); 327 check_set_get (-d, x); 328 329 /* check the smallest power of two */ 330 d = 1.0; 331 while ((e = d / 2.0) != (long double) 0.0 && e != d) 332 d = e; 333 check_set_get (d, x); 334 check_set_get (-d, x); 335 336 /* check largest 2^(2^k) that is representable as a long double */ 337 d = (LDBL_MAX / 2) + (LDBL_MAX / 4 * LDBL_EPSILON); 338 check_set_get (d, x); 339 340 /* check that 2^i, 2^i+1 and 2^i-1 are correctly converted */ 341 d = 1.0; 342 for (i = 1; i < MPFR_LDBL_MANT_DIG; i++) 343 { 344 d = 2.0 * d; /* d = 2^i */ 345 check_set_get (d, x); 346 check_set_get (d + 1.0, x); 347 check_set_get (d - 1.0, x); 348 } 349 350 for (i = 0; i < 10000; i++) 351 { 352 mpfr_urandomb (x, RANDS); 353 d = mpfr_get_ld (x, MPFR_RNDN); 354 check_set_get (d, x); 355 } 356 357 /* check with reduced emax to exercise overflow */ 358 emax = mpfr_get_emax (); 359 mpfr_set_prec (x, 2); 360 set_emax (1); 361 mpfr_set_ld (x, (long double) 2.0, MPFR_RNDN); 362 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 363 for (d = (long double) 2.0, i = 0; i < 13; i++, d *= d); 364 /* now d = 2^8192 */ 365 mpfr_set_ld (x, d, MPFR_RNDN); 366 MPFR_ASSERTN(mpfr_inf_p (x) && mpfr_sgn (x) > 0); 367 set_emax (emax); 368 369 mpfr_clear (x); 370 371 test_small (); 372 373 bug_20160907 (); 374 375 tests_end_mpfr (); 376 377 return 0; 378 } 379