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