1 /* tfmod -- test file for mpfr_fmod 2 3 Copyright 2007-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_fmod 26 #define TWO_ARGS 27 #include "tgeneric.c" 28 29 /* compute remainder as in definition: 30 r = x - n * y, where n = trunc(x/y). 31 warning: may change flags. */ 32 static int 33 slow_fmod (mpfr_ptr r, mpfr_srcptr x, mpfr_srcptr y, mpfr_rnd_t rnd) 34 { 35 mpfr_t q; 36 int inexact; 37 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x) || MPFR_IS_SINGULAR (y))) 38 { 39 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y) || MPFR_IS_INF (x) 40 || MPFR_IS_ZERO (y)) 41 { 42 MPFR_SET_NAN (r); 43 MPFR_RET_NAN; 44 } 45 else /* either y is Inf and x is 0 or non-special, 46 or x is 0 and y is non-special, 47 in both cases the quotient is zero. */ 48 return mpfr_set (r, x, rnd); 49 } 50 /* regular cases */ 51 /* if 2^(ex-1) <= |x| < 2^ex, and 2^(ey-1) <= |y| < 2^ey, 52 then |x/y| < 2^(ex-ey+1) */ 53 mpfr_init2 (q, 54 MAX (MPFR_PREC_MIN, mpfr_get_exp (x) - mpfr_get_exp (y) + 1)); 55 mpfr_div (q, x, y, MPFR_RNDZ); 56 mpfr_trunc (q, q); /* may change inexact flag */ 57 mpfr_prec_round (q, mpfr_get_prec (q) + mpfr_get_prec (y), MPFR_RNDZ); 58 inexact = mpfr_mul (q, q, y, MPFR_RNDZ); /* exact */ 59 inexact = mpfr_sub (r, x, q, rnd); 60 mpfr_clear (q); 61 return inexact; 62 } 63 64 static void 65 test_failed (mpfr_ptr erem, mpfr_ptr grem, int eret, int gret, 66 mpfr_ptr x, mpfr_ptr y, mpfr_rnd_t rnd) 67 { 68 printf ("error: mpfr_fmod (r, x, y, rnd)\n x = "); 69 mpfr_out_str (stdout, 10, 0, x, MPFR_RNDD); 70 printf ("\n y = "); 71 mpfr_out_str (stdout, 10, 0, y, MPFR_RNDD); 72 printf ("\nrnd = %s", mpfr_print_rnd_mode (rnd)); 73 if (eret != gret) 74 printf ("\nexpected %s return value, got %d", 75 (eret < 0 ? "negative" : eret > 0 ? "positive" : "zero"), gret); 76 printf ("\n expected r = "); 77 mpfr_out_str (stdout, 10, 0, erem, MPFR_RNDD); 78 printf ("\n got r = "); 79 mpfr_out_str (stdout, 10, 0, grem, MPFR_RNDD); 80 putchar ('\n'); 81 82 exit (1); 83 } 84 85 static void 86 check (mpfr_ptr r0, mpfr_ptr x, mpfr_ptr y, mpfr_rnd_t rnd) 87 { 88 int inex0, inex1; 89 mpfr_t r1; 90 mpfr_init2 (r1, mpfr_get_prec (r0)); 91 92 inex0 = mpfr_fmod (r0, x, y, rnd); 93 inex1 = slow_fmod (r1, x, y, rnd); 94 if (!mpfr_equal_p (r0, r1) || inex0 != inex1) 95 test_failed (r1, r0, inex1, inex0, x, y, rnd); 96 mpfr_clear (r1); 97 } 98 99 static void 100 regular (void) 101 { 102 mpfr_t x, y, r; 103 mpfr_inits (x, y, r, (mpfr_ptr) 0); 104 105 /* remainder = 0 */ 106 mpfr_set_str (y, "FEDCBA987654321p-64", 16, MPFR_RNDN); 107 mpfr_pow_ui (x, y, 42, MPFR_RNDN); 108 check (r, x, y, MPFR_RNDN); 109 110 /* x < y */ 111 mpfr_set_ui_2exp (x, 64723, -19, MPFR_RNDN); 112 mpfr_mul (x, x, y, MPFR_RNDN); 113 check (r, x, y, MPFR_RNDN); 114 115 /* sign(x) = sign (r) */ 116 mpfr_set_ui (x, 123798, MPFR_RNDN); 117 mpfr_set_ui (y, 10, MPFR_RNDN); 118 check (r, x, y, MPFR_RNDN); 119 120 /* huge difference between precisions */ 121 mpfr_set_prec (x, 314); 122 mpfr_set_prec (y, 8); 123 mpfr_set_prec (r, 123); 124 mpfr_const_pi (x, MPFR_RNDD); /* x = pi */ 125 mpfr_set_ui_2exp (y, 1, 3, MPFR_RNDD); /* y = 1/8 */ 126 check (r, x, y, MPFR_RNDD); 127 128 mpfr_clears (x, y, r, (mpfr_ptr) 0); 129 } 130 131 static void 132 special (void) 133 { 134 int inexact; 135 mpfr_t x, y, r, t; 136 137 mpfr_inits (x, y, r, t, (mpfr_ptr) 0); 138 139 mpfr_set_nan (t); 140 141 /* fmod (NaN, NaN) is NaN */ 142 mpfr_set_nan (x); 143 mpfr_set_nan (y); 144 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 145 if (!mpfr_nan_p (r) || inexact != 0) 146 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 147 148 /* fmod (NaN, +0) is NaN */ 149 mpfr_set_ui (y, 0, MPFR_RNDN); 150 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 151 if (!mpfr_nan_p (r) || inexact != 0) 152 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 153 154 /* fmod (+1, 0) is NaN */ 155 mpfr_set_ui (x, 1, MPFR_RNDN); 156 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 157 if (!mpfr_nan_p (r) || inexact != 0) 158 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 159 160 /* fmod (0, 0) is NaN */ 161 mpfr_set_ui (x, 0, MPFR_RNDN); 162 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 163 if (!mpfr_nan_p (r) || inexact != 0) 164 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 165 166 /* fmod (+inf, +0) is NaN */ 167 mpfr_set_inf (x, +1); 168 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 169 if (!mpfr_nan_p (r) || inexact != 0) 170 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 171 172 /* fmod (-inf, +0) is NaN */ 173 mpfr_set_inf (x, -1); 174 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 175 if (!mpfr_nan_p (r) || inexact != 0) 176 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 177 178 /* fmod (-inf, -0) is NaN */ 179 mpfr_neg (x, x, MPFR_RNDN); 180 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 181 if (!mpfr_nan_p (r) || inexact != 0) 182 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 183 184 /* fmod (-inf, +1) is NaN */ 185 mpfr_set_ui (y, +1, MPFR_RNDN); 186 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 187 if (!mpfr_nan_p (r) || inexact != 0) 188 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 189 190 /* fmod (+inf, +1) is NaN */ 191 mpfr_neg (x, x, MPFR_RNDN); 192 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 193 if (!mpfr_nan_p (r) || inexact != 0) 194 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 195 196 /* fmod (+inf, -inf) is NaN */ 197 mpfr_set_inf (y, -1); 198 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 199 if (!mpfr_nan_p (r) || inexact != 0) 200 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 201 202 /* fmod (-inf, -inf) is NaN */ 203 mpfr_neg (x, x, MPFR_RNDN); 204 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 205 if (!mpfr_nan_p (r) || inexact != 0) 206 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 207 208 /* fmod (-inf, +inf) is NaN */ 209 mpfr_neg (y, y, MPFR_RNDN); 210 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 211 if (!mpfr_nan_p (r) || inexact != 0) 212 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 213 214 /* fmod (+inf, +inf) is NaN */ 215 mpfr_neg (x, x, MPFR_RNDN); 216 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 217 if (!mpfr_nan_p (r) || inexact != 0) 218 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 219 220 /* fmod (x, +inf) = x, if x is finite */ 221 mpfr_set_ui (x, 1, MPFR_RNDN); 222 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 223 if (!mpfr_equal_p (r, x) || inexact != 0) 224 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 225 mpfr_neg (x, x, MPFR_RNDN); 226 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 227 if (!mpfr_equal_p (r, x) || inexact != 0) 228 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 229 230 /* fmod (+0, +inf) = +0 */ 231 mpfr_set_ui (x, 0, MPFR_RNDN); 232 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 233 if (!mpfr_equal_p (r, x) || inexact != 0) 234 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 235 236 /* fmod (-0, +inf) = -0 */ 237 mpfr_neg (x, x, MPFR_RNDN); 238 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 239 if (!mpfr_equal_p (r, x) || inexact != 0) 240 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 241 242 /* fmod (x, -inf) = x, if x is finite */ 243 mpfr_set_inf (y, -1); 244 mpfr_set_ui (x, 1, MPFR_RNDN); 245 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 246 if (!mpfr_equal_p (r, x) || inexact != 0) 247 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 248 mpfr_neg (x, x, MPFR_RNDN); 249 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 250 if (!mpfr_equal_p (r, x) || inexact != 0) 251 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 252 253 /* fmod (+0, -inf) = +0 */ 254 mpfr_set_ui (x, 0, MPFR_RNDN); 255 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 256 if (!mpfr_equal_p (r, x) || inexact != 0) 257 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 258 259 /* fmod (-0, -inf) = -0 */ 260 mpfr_neg (x, x, MPFR_RNDN); 261 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 262 if (!mpfr_equal_p (r, x) || inexact != 0) 263 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 264 265 /* fmod (+0, +0) is NaN */ 266 mpfr_set_ui (x, 0, MPFR_RNDN); 267 mpfr_set_ui (y, 0, MPFR_RNDN); 268 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 269 if (!mpfr_nan_p (r) || inexact != 0) 270 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 271 272 /* fmod (+0, -0) is NaN */ 273 mpfr_neg (y, y, MPFR_RNDN); 274 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 275 if (!mpfr_nan_p (r) || inexact != 0) 276 test_failed (r, t, 0, inexact, x, y, MPFR_RNDN); 277 278 /* fmod (+0, +1) = +0 */ 279 mpfr_set_ui (y, 1, MPFR_RNDN); 280 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 281 if (!mpfr_equal_p (r, x) || inexact != 0) 282 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 283 284 /* fmod (+0, -1) = +0 */ 285 mpfr_neg (y, y, MPFR_RNDN); 286 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 287 if (!mpfr_equal_p (r, x) || inexact != 0) 288 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 289 290 /* fmod (-0, -1) = -0 */ 291 mpfr_neg (x, x, MPFR_RNDN); 292 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 293 if (!mpfr_equal_p (r, x) || inexact != 0) 294 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 295 296 /* fmod (-0, +1) = -0 */ 297 mpfr_neg (y, y, MPFR_RNDN); 298 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 299 if (!mpfr_equal_p (r, x) || inexact != 0) 300 test_failed (r, x, 0, inexact, x, y, MPFR_RNDN); 301 302 mpfr_set_prec (x, 380); 303 mpfr_set_prec (y, 385); 304 mpfr_set_str_binary (x, "0.11011010010110011101011000100100101100101011010001011100110001100101111001010100001011111110111100101110101010110011010101000100000100011101101100001011101110100111101111111010001001000010000110010110011100111000001110111010000100101001010111100100010001101001110100011110010000000001110001111001101100111011001000110110011100100011111110010100011001000001001011010111010000000000E0"); 305 mpfr_set_str_binary (y, "0.1100011000011101011010001100010111001110110111001101010010111100111100011010010011011101111101111001010111111110001001100001111101001000000010100101111001001110010110000111001000101010111001001000100101011111000010100110001111000110011011010101111101100110010101011010011101100001011101001000101111110110110110000001001101110111110110111110111111001001011110001110011111100000000000000E-1"); 306 mpfr_set_prec (r, 2); 307 inexact = mpfr_fmod (r, x, y, MPFR_RNDA); 308 mpfr_set_prec (t, 2); 309 mpfr_set_ui_2exp (t, 3, -5, MPFR_RNDN); 310 if (mpfr_cmp_ui_2exp (r, 3, -5) || inexact <= 0) 311 test_failed (r, t, 1, inexact, x, y, MPFR_RNDA); 312 313 mpfr_clears (x, y, r, t, (mpfr_ptr) 0); 314 return; 315 } 316 317 /* bug reported by Eric Veach */ 318 static void 319 bug20090519 (void) 320 { 321 mpfr_t x, y, r; 322 int inexact; 323 324 mpfr_inits2 (100, x, y, r, (mpfr_ptr) 0); 325 326 mpfr_set_prec (x, 3); 327 mpfr_set_prec (y, 3); 328 mpfr_set_prec (r, 3); 329 mpfr_set_si (x, 8, MPFR_RNDN); 330 mpfr_set_si (y, 7, MPFR_RNDN); 331 check (r, x, y, MPFR_RNDN); 332 333 mpfr_set_prec (x, 10); 334 mpfr_set_prec (y, 10); 335 mpfr_set_prec (r, 10); 336 mpfr_set_ui_2exp (x, 3, 26, MPFR_RNDN); 337 mpfr_set_si (y, (1 << 9) - 1, MPFR_RNDN); 338 check (r, x, y, MPFR_RNDN); 339 340 mpfr_set_prec (x, 100); 341 mpfr_set_prec (y, 100); 342 mpfr_set_prec (r, 100); 343 mpfr_set_str (x, "3.5", 10, MPFR_RNDN); 344 mpfr_set_str (y, "1.1", 10, MPFR_RNDN); 345 check (r, x, y, MPFR_RNDN); 346 /* double check, with a pre-computed value */ 347 { 348 mpfr_t er; 349 mpfr_init2 (er, 100); 350 mpfr_set_str (er, "CCCCCCCCCCCCCCCCCCCCCCCC8p-102", 16, MPFR_RNDN); 351 352 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 353 if (!mpfr_equal_p (r, er) || inexact != 0) 354 test_failed (er, r, 0, inexact, x, y, MPFR_RNDN); 355 356 mpfr_clear (er); 357 } 358 359 mpfr_set_si (x, 20, MPFR_RNDN); 360 mpfr_set_ui_2exp (y, 1, 1, MPFR_RNDN); /* exact */ 361 mpfr_sin (y, y, MPFR_RNDN); 362 check (r, x, y, MPFR_RNDN); 363 364 mpfr_clears (x, y, r, (mpfr_ptr) 0); 365 } 366 367 static void 368 bug20160217 (void) 369 { 370 mpfr_t x, y, r; 371 int inexact, i; 372 mpfr_exp_t emin, emax; 373 374 mpfr_inits2 (53, x, y, r, (mpfr_ptr) 0); 375 376 emin = mpfr_get_emin (); 377 emax = mpfr_get_emax (); 378 379 for (i = 0; i <= 1; i++) 380 { 381 mpfr_set_zero (x, 1); 382 mpfr_nextabove (x); 383 mpfr_set_inf (y, 1); 384 mpfr_nextbelow (y); 385 inexact = mpfr_fmod (r, x, y, MPFR_RNDN); 386 if (!mpfr_equal_p (r, x) || inexact != 0) 387 { 388 printf ("Error for mpfr_fmod (r, nextabove(0), nextbelow(+inf)," 389 " MPFR_RNDN)%s\n", i ? "extended exponent range" : ""); 390 printf ("Expected inex = 0, r = "); 391 mpfr_dump (x); 392 printf ("Got inex = %d, r = ", inexact); 393 mpfr_dump (r); 394 exit (1); 395 } 396 set_emin (MPFR_EMIN_MIN); 397 set_emax (MPFR_EMAX_MAX); 398 } 399 400 set_emin (emin); 401 set_emax (emax); 402 403 mpfr_clears (x, y, r, (mpfr_ptr) 0); 404 } 405 406 int 407 main (int argc, char *argv[]) 408 { 409 tests_start_mpfr (); 410 411 bug20090519 (); 412 bug20160217 (); 413 414 test_generic (MPFR_PREC_MIN, 100, 100); 415 416 special (); 417 regular (); 418 419 tests_end_mpfr (); 420 return 0; 421 } 422