1/* Test file for mpfr_subnormalize. 2 3Copyright 2005-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 25static const struct { 26 const char *in; 27 int i; 28 mpfr_rnd_t rnd; 29 const char *out; 30 int j; 31} tab[] = { /* 4th field: use the mpfr_dump format, in case of error. */ 32 {"1E1", 0, MPFR_RNDN, "0.100000000E2", 0}, 33 {"1E1", -1, MPFR_RNDZ, "0.100000000E2", -1}, 34 {"1E1", -1, MPFR_RNDD, "0.100000000E2", -1}, 35 {"1E1", 1, MPFR_RNDU, "0.100000000E2", 1}, 36 {"0.10000E-10", 0, MPFR_RNDN, "0.100000000E-10", 0}, 37 {"0.10001E-10", 0, MPFR_RNDN, "0.100000000E-10", -1}, 38 {"0.11001E-10", 0, MPFR_RNDN, "0.100000000E-9", 1}, 39 {"0.11001E-10", 0, MPFR_RNDZ, "0.100000000E-10", -1}, 40 {"0.11001E-10", 0, MPFR_RNDU, "0.100000000E-9", 1}, 41 {"0.11000E-10", 0, MPFR_RNDN, "0.100000000E-9", 1}, 42 {"0.11000E-10", -1, MPFR_RNDN, "0.100000000E-9", 1}, 43 {"0.11000E-10", 1, MPFR_RNDN, "0.100000000E-10", -1}, 44 {"0.11111E-8", 0, MPFR_RNDN, "0.100000000E-7", 1}, 45 {"0.10111E-8", 0, MPFR_RNDN, "0.110000000E-8", 1}, 46 {"0.11110E-8", -1, MPFR_RNDN, "0.100000000E-7", 1}, 47 {"0.10110E-8", 1, MPFR_RNDN, "0.101000000E-8", -1} 48}; 49 50static void 51check1 (void) 52{ 53 mpfr_t x; 54 int i, j, k, s, old_inex, tiny, expj; 55 mpfr_exp_t emin, emax; 56 unsigned int expflags, flags; 57 58 emin = mpfr_get_emin (); 59 emax = mpfr_get_emax (); 60 61 mpfr_set_default_prec (9); 62 set_emin (-10); 63 set_emax (10); 64 65 mpfr_init (x); 66 for (i = 0; i < numberof (tab); i++) 67 for (s = 0; s <= (tab[i].rnd == MPFR_RNDN); s++) 68 for (k = 0; k <= 1; k++) 69 { 70 mpfr_set_str (x, tab[i].in, 2, MPFR_RNDN); 71 old_inex = tab[i].i; 72 expj = tab[i].j; 73 if (s) 74 { 75 mpfr_neg (x, x, MPFR_RNDN); 76 old_inex = - old_inex; 77 expj = - expj; 78 } 79 if (k && old_inex) 80 old_inex = old_inex < 0 ? INT_MIN : INT_MAX; 81 tiny = MPFR_GET_EXP (x) <= -3; 82 mpfr_clear_flags (); 83 j = mpfr_subnormalize (x, old_inex, tab[i].rnd); 84 expflags = 85 (tiny ? MPFR_FLAGS_UNDERFLOW : 0) | 86 (expj ? MPFR_FLAGS_INEXACT : 0); 87 flags = __gmpfr_flags; 88 if (s) 89 mpfr_neg (x, x, MPFR_RNDN); 90 if (mpfr_cmp_str (x, tab[i].out, 2, MPFR_RNDN) != 0 || 91 flags != expflags || ! SAME_SIGN (j, expj)) 92 { 93 const char *sgn = s ? "-" : ""; 94 printf ("Error for i = %d (old_inex = %d), k = %d, x = %s%s\n" 95 "Expected: %s%s\nGot: ", i, old_inex, k, 96 sgn, tab[i].in, sgn, tab[i].out); 97 if (s) 98 mpfr_neg (x, x, MPFR_RNDN); 99 mpfr_dump (x); 100 printf ("Expected flags = %u, got %u\n", expflags, flags); 101 printf ("Expected ternary value = %d, got %d\n", expj, j); 102 exit (1); 103 } 104 } 105 mpfr_clear (x); 106 107 MPFR_ASSERTN (mpfr_get_emin () == -10); 108 MPFR_ASSERTN (mpfr_get_emax () == 10); 109 110 set_emin (emin); 111 set_emax (emax); 112} 113 114/* bug found by Kevin P. Rauch on 22 Oct 2007 */ 115static void 116check2 (void) 117{ 118 mpfr_t x, y, z; 119 int tern; 120 mpfr_exp_t emin; 121 122 emin = mpfr_get_emin (); 123 124 mpfr_init2 (x, 32); 125 mpfr_init2 (y, 32); 126 mpfr_init2 (z, 32); 127 128 mpfr_set_ui (x, 0xC0000000U, MPFR_RNDN); 129 mpfr_neg (x, x, MPFR_RNDN); 130 mpfr_set_ui (y, 0xFFFFFFFEU, MPFR_RNDN); 131 mpfr_set_exp (x, 0); 132 mpfr_set_exp (y, 0); 133 set_emin (-29); 134 135 tern = mpfr_mul (z, x, y, MPFR_RNDN); 136 /* z = -0.BFFFFFFE, tern > 0 */ 137 138 tern = mpfr_subnormalize (z, tern, MPFR_RNDN); 139 /* z should be -0.75 */ 140 MPFR_ASSERTN (tern < 0 && mpfr_cmp_si_2exp (z, -3, -2) == 0); 141 142 mpfr_clear (x); 143 mpfr_clear (y); 144 mpfr_clear (z); 145 146 MPFR_ASSERTN (mpfr_get_emin () == -29); 147 148 set_emin (emin); 149} 150 151/* bug found by Kevin P. Rauch on 22 Oct 2007 */ 152static void 153check3 (void) 154{ 155 mpfr_t x, y, z; 156 int tern; 157 mpfr_exp_t emin; 158 159 emin = mpfr_get_emin (); 160 161 mpfr_init2 (x, 32); 162 mpfr_init2 (y, 32); 163 mpfr_init2 (z, 32); 164 165 mpfr_set_ui (x, 0xBFFFFFFFU, MPFR_RNDN); /* 3221225471/2^32 */ 166 mpfr_set_ui (y, 0x80000001U, MPFR_RNDN); /* 2147483649/2^32 */ 167 mpfr_set_exp (x, 0); 168 mpfr_set_exp (y, 0); 169 set_emin (-1); 170 171 /* the exact product is 6917529028714823679/2^64, which is rounded to 172 3/8 = 0.375, which is smaller, thus tern < 0 */ 173 tern = mpfr_mul (z, x, y, MPFR_RNDN); 174 MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0); 175 176 tern = mpfr_subnormalize (z, tern, MPFR_RNDN); 177 /* since emin = -1, and EXP(z)=-1, z should be rounded to precision 178 EXP(z)-emin+1 = 1, i.e., z should be a multiple of the smallest possible 179 positive representable value with emin=-1, which is 1/4. The two 180 possible values are 1/4 and 2/4, which are at equal distance of z. 181 But since tern < 0, we should choose the largest value, i.e., 2/4. */ 182 MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0); 183 184 /* here is another test for the alternate case, where z was rounded up 185 first, thus we have to round down */ 186 mpfr_set_str_binary (x, "0.11111111111010110101011011011011"); 187 mpfr_set_str_binary (y, "0.01100000000001111100000000001110"); 188 tern = mpfr_mul (z, x, y, MPFR_RNDN); 189 MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 3, -3) == 0); 190 tern = mpfr_subnormalize (z, tern, MPFR_RNDN); 191 MPFR_ASSERTN (tern < 0 && mpfr_cmp_ui_2exp (z, 1, -2) == 0); 192 193 /* finally the case where z was exact, which we simulate here */ 194 mpfr_set_ui_2exp (z, 3, -3, MPFR_RNDN); 195 tern = mpfr_subnormalize (z, 0, MPFR_RNDN); 196 MPFR_ASSERTN (tern > 0 && mpfr_cmp_ui_2exp (z, 1, -1) == 0); 197 198 mpfr_clear (x); 199 mpfr_clear (y); 200 mpfr_clear (z); 201 202 MPFR_ASSERTN (mpfr_get_emin () == -1); 203 204 set_emin (emin); 205} 206 207/* check mpfr_subnormize with RNDNA (experimental) */ 208static void 209check_rndna (void) 210{ 211 mpfr_t x, y; 212 int inex; 213 mpfr_exp_t emin = mpfr_get_emin (); 214 215 mpfr_init2 (x, 11); 216 mpfr_init2 (y, 9); 217 218 mpfr_set_str_binary (x, "0.1111111010E-14"); 219 inex = mpfr_set (y, x, MPFR_RNDN); 220 MPFR_ASSERTN(inex == 0); 221 set_emin (-21); 222 inex = mpfr_subnormalize (y, inex, MPFR_RNDNA); 223 /* mpfr_subnormalize rounds to precision EXP(y)-emin+1 = 8, 224 thus should round to 0.111111110E-14 */ 225 mpfr_set_str_binary (x, "0.111111110E-14"); 226 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 227 MPFR_ASSERTN(inex > 0); 228 229 /* now consider x slightly larger: we should get the same result */ 230 mpfr_set_str_binary (x, "0.1111111011E-14"); 231 inex = mpfr_set (y, x, MPFR_RNDN); 232 MPFR_ASSERTN(inex > 0); 233 inex = mpfr_subnormalize (y, inex, MPFR_RNDNA); 234 mpfr_set_str_binary (x, "0.111111110E-14"); 235 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 236 MPFR_ASSERTN(inex > 0); 237 238 /* now consider x slightly smaller: we should get a different result */ 239 mpfr_set_str_binary (x, "0.11111110001E-14"); 240 inex = mpfr_set (y, x, MPFR_RNDN); 241 MPFR_ASSERTN(inex < 0); 242 inex = mpfr_subnormalize (y, inex, MPFR_RNDNA); 243 mpfr_set_str_binary (x, "0.111111100E-14"); 244 MPFR_ASSERTN(mpfr_cmp (y, x) == 0); 245 MPFR_ASSERTN(inex < 0); 246 247 mpfr_clear (x); 248 mpfr_clear (y); 249 set_emin (emin); 250} 251 252/* exercise a corner case of mpfr_subnormalize: 253 y = 1xxx...xxx0|100000| with old_inexact = -1 */ 254static void 255coverage (void) 256{ 257 mpfr_t y; 258 int inex; 259 260 mpfr_init2 (y, 42); 261 mpfr_set_ui_2exp (y, 131073, mpfr_get_emin () - 2, MPFR_RNDN); 262 MPFR_ASSERTN(mpfr_get_exp (y) == mpfr_get_emin () + 16); 263 /* mpfr_subnormalize rounds y to precision EXP(y) - emin + 1, thus 17 */ 264 inex = mpfr_subnormalize (y, -1, MPFR_RNDN); 265 MPFR_ASSERTN(mpfr_cmp_ui_2exp (y, 65537, mpfr_get_emin () - 1) == 0); 266 MPFR_ASSERTN(inex > 0); 267 mpfr_clear (y); 268} 269 270int 271main (int argc, char *argv[]) 272{ 273 tests_start_mpfr (); 274 275 coverage (); 276 check_rndna (); 277 check1 (); 278 check2 (); 279 check3 (); 280 281 tests_end_mpfr (); 282 return 0; 283} 284