| 1 |
/* Test of the double rounding effect. |
| 2 |
* |
| 3 |
* This example was presented at the CNC'2 summer school on MPFR and MPC |
| 4 |
* at LORIA, Nancy, France. |
| 5 |
* |
| 6 |
* Arguments: max difference of exponents dmax, significand size n. |
| 7 |
* Optional argument: extended precision p (with double rounding). |
| 8 |
* |
| 9 |
* Return all the couples of positive machine numbers (x,y) such that |
| 10 |
* 1/2 <= y < 1, 0 <= Ex - Ey <= dmax, x - y is exactly representable |
| 11 |
* in precision n and the results of floor(x/y) in the rounding modes |
| 12 |
* toward 0 and to nearest are different. |
| 13 |
*/ |
| 14 |
|
| 15 |
/* |
| 16 |
Copyright 2009-2020 Free Software Foundation, Inc. |
| 17 |
Contributed by the AriC and Caramba projects, INRIA. |
| 18 |
|
| 19 |
This file is part of the GNU MPFR Library. |
| 20 |
|
| 21 |
The GNU MPFR Library is free software; you can redistribute it and/or modify |
| 22 |
it under the terms of the GNU Lesser General Public License as published by |
| 23 |
the Free Software Foundation; either version 3 of the License, or (at your |
| 24 |
option) any later version. |
| 25 |
|
| 26 |
The GNU MPFR Library is distributed in the hope that it will be useful, but |
| 27 |
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
| 28 |
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |
| 29 |
License for more details. |
| 30 |
|
| 31 |
You should have received a copy of the GNU Lesser General Public License |
| 32 |
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see |
| 33 |
https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., |
| 34 |
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. |
| 35 |
*/ |
| 36 |
|
| 37 |
#include <stdio.h> |
| 38 |
#include <stdlib.h> |
| 39 |
#include <mpfr.h> |
| 40 |
|
| 41 |
#define PRECN x, y, z |
| 42 |
#define VARS PRECN, t |
| 43 |
|
| 44 |
static unsigned long |
| 45 |
eval (mpfr_t x, mpfr_t y, mpfr_t z, mpfr_t t, mpfr_rnd_t rnd) |
| 46 |
{ |
| 47 |
mpfr_div (t, x, y, rnd); /* the division x/y in precision p */ |
| 48 |
mpfr_set (z, t, rnd); /* the rounding to the precision n */ |
| 49 |
mpfr_rint_floor (z, z, rnd); |
| 50 |
return mpfr_get_ui (z, rnd); |
| 51 |
} |
| 52 |
|
| 53 |
int main (int argc, char *argv[]) |
| 54 |
{ |
| 55 |
int dmax, n, p; |
| 56 |
mpfr_t VARS; |
| 57 |
|
| 58 |
if (argc != 3 && argc != 4) |
| 59 |
{ |
| 60 |
fprintf (stderr, "Usage: divworst <dmax> <n> [ <p> ]\n"); |
| 61 |
exit (EXIT_FAILURE); |
| 62 |
} |
| 63 |
|
| 64 |
dmax = atoi (argv[1]); |
| 65 |
n = atoi (argv[2]); |
| 66 |
p = argc == 3 ? n : atoi (argv[3]); |
| 67 |
if (p < n) |
| 68 |
{ |
| 69 |
fprintf (stderr, "divworst: p must be greater or equal to n\n"); |
| 70 |
exit (EXIT_FAILURE); |
| 71 |
} |
| 72 |
|
| 73 |
mpfr_inits2 (n, PRECN, (mpfr_ptr) 0); |
| 74 |
mpfr_init2 (t, p); |
| 75 |
|
| 76 |
for (mpfr_set_ui_2exp (x, 1, -1, MPFR_RNDN); |
| 77 |
mpfr_get_exp (x) <= dmax; |
| 78 |
mpfr_nextabove (x)) |
| 79 |
for (mpfr_set_ui_2exp (y, 1, -1, MPFR_RNDN); |
| 80 |
mpfr_get_exp (y) == 0; |
| 81 |
mpfr_nextabove (y)) |
| 82 |
{ |
| 83 |
unsigned long rz, rn; |
| 84 |
|
| 85 |
if (mpfr_sub (z, x, y, MPFR_RNDZ) != 0) |
| 86 |
continue; /* x - y is not representable in precision n */ |
| 87 |
rz = eval (x, y, z, t, MPFR_RNDZ); |
| 88 |
rn = eval (x, y, z, t, MPFR_RNDN); |
| 89 |
if (rz == rn) |
| 90 |
continue; |
| 91 |
mpfr_printf ("x = %.*Rb ; y = %.*Rb ; Z: %lu ; N: %lu\n", |
| 92 |
n - 1, x, n - 1, y, rz, rn); |
| 93 |
} |
| 94 |
|
| 95 |
mpfr_clears (VARS, (mpfr_ptr) 0); |
| 96 |
mpfr_free_cache (); |
| 97 |
return 0; |
| 98 |
} |
