/* Test file for mpfr_fmma and mpfr_fmms.
Copyright 2016-2020 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.
This file is part of the GNU MPFR Library.
The GNU MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.
The GNU MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
#include "mpfr-test.h"
/* check both mpfr_fmma and mpfr_fmms */
static void
random_test (mpfr_t a, mpfr_t b, mpfr_t c, mpfr_t d, mpfr_rnd_t rnd)
{
mpfr_t ref, res, ab, cd;
int inex_ref, inex_res;
mpfr_prec_t p = MPFR_PREC(a);
mpfr_init2 (res, p);
mpfr_init2 (ref, p);
mpfr_init2 (ab, mpfr_get_prec (a) + mpfr_get_prec (b));
mpfr_init2 (cd, mpfr_get_prec (c) + mpfr_get_prec (d));
/* first check fmma */
inex_res = mpfr_fmma (res, a, b, c, d, rnd);
mpfr_mul (ab, a, b, rnd);
mpfr_mul (cd, c, d, rnd);
inex_ref = mpfr_add (ref, ab, cd, rnd);
if (! SAME_SIGN (inex_res, inex_ref) ||
mpfr_nan_p (res) || mpfr_nan_p (ref) ||
! mpfr_equal_p (res, ref))
{
printf ("mpfr_fmma failed for p=%ld rnd=%s\n",
(long int) p, mpfr_print_rnd_mode (rnd));
printf ("a="); mpfr_dump (a);
printf ("b="); mpfr_dump (b);
printf ("c="); mpfr_dump (c);
printf ("d="); mpfr_dump (d);
printf ("Expected\n ");
mpfr_dump (ref);
printf (" with inex = %d\n", inex_ref);
printf ("Got\n ");
mpfr_dump (res);
printf (" with inex = %d\n", inex_res);
exit (1);
}
/* then check fmms */
inex_res = mpfr_fmms (res, a, b, c, d, rnd);
mpfr_mul (ab, a, b, rnd);
mpfr_mul (cd, c, d, rnd);
inex_ref = mpfr_sub (ref, ab, cd, rnd);
if (! SAME_SIGN (inex_res, inex_ref) ||
mpfr_nan_p (res) || mpfr_nan_p (ref) ||
! mpfr_equal_p (res, ref))
{
printf ("mpfr_fmms failed for p=%ld rnd=%s\n",
(long int) p, mpfr_print_rnd_mode (rnd));
printf ("a="); mpfr_dump (a);
printf ("b="); mpfr_dump (b);
printf ("c="); mpfr_dump (c);
printf ("d="); mpfr_dump (d);
printf ("Expected\n ");
mpfr_dump (ref);
printf (" with inex = %d\n", inex_ref);
printf ("Got\n ");
mpfr_dump (res);
printf (" with inex = %d\n", inex_res);
exit (1);
}
mpfr_clear (ab);
mpfr_clear (cd);
mpfr_clear (res);
mpfr_clear (ref);
}
static void
random_tests (void)
{
mpfr_prec_t p;
int r;
mpfr_t a, b, c, d;
for (p = MPFR_PREC_MIN; p <= 4096; p++)
{
mpfr_inits2 (p, a, b, c, d, (mpfr_ptr) 0);
mpfr_urandomb (a, RANDS);
mpfr_urandomb (b, RANDS);
mpfr_urandomb (c, RANDS);
mpfr_urandomb (d, RANDS);
RND_LOOP_NO_RNDF (r)
random_test (a, b, c, d, (mpfr_rnd_t) r);
mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
}
}
static void
zero_tests (void)
{
mpfr_exp_t emin, emax;
mpfr_t a, b, c, d, res;
int i, r;
emin = mpfr_get_emin ();
emax = mpfr_get_emax ();
set_emin (MPFR_EMIN_MIN);
set_emax (MPFR_EMAX_MAX);
mpfr_inits2 (GMP_NUMB_BITS, a, b, c, d, (mpfr_ptr) 0);
for (i = 0; i <= 4; i++)
{
switch (i)
{
case 0: case 1:
mpfr_set_ui (a, i, MPFR_RNDN);
break;
case 2:
mpfr_setmax (a, MPFR_EMAX_MAX);
break;
case 3:
mpfr_setmin (a, MPFR_EMIN_MIN);
break;
case 4:
mpfr_setmin (a, MPFR_EMIN_MIN+1);
break;
}
RND_LOOP (r)
{
int j, inex;
mpfr_flags_t flags;
mpfr_set (b, a, MPFR_RNDN);
mpfr_set (c, a, MPFR_RNDN);
mpfr_neg (d, a, MPFR_RNDN);
/* We also want to test cases where the precision of the
result is twice the precision of each input, in case
add1sp.c and/or sub1sp.c could be involved. */
for (j = 1; j <= 2; j++)
{
mpfr_init2 (res, GMP_NUMB_BITS * j);
mpfr_clear_flags ();
inex = mpfr_fmma (res, a, b, c, d, (mpfr_rnd_t) r);
flags = __gmpfr_flags;
if (flags != 0 || inex != 0 || ! mpfr_zero_p (res) ||
(r == MPFR_RNDD ? MPFR_IS_POS (res) : MPFR_IS_NEG (res)))
{
printf ("Error in zero_tests on i = %d, %s\n",
i, mpfr_print_rnd_mode ((mpfr_rnd_t) r));
printf ("Expected %c0, inex = 0\n",
r == MPFR_RNDD ? '-' : '+');
printf ("Got ");
if (MPFR_IS_POS (res))
printf ("+");
mpfr_out_str (stdout, 16, 0, res, MPFR_RNDN);
printf (", inex = %d\n", inex);
printf ("Expected flags:");
flags_out (0);
printf ("Got flags: ");
flags_out (flags);
exit (1);
}
mpfr_clear (res);
} /* j */
} /* r */
} /* i */
mpfr_clears (a, b, c, d, (mpfr_ptr) 0);
set_emin (emin);
set_emax (emax);
}
static void
max_tests (void)
{
mpfr_exp_t emin, emax;
mpfr_t x, y1, y2;
int r;
int i, inex1, inex2;
mpfr_flags_t flags1, flags2;
emin = mpfr_get_emin ();
emax = mpfr_get_emax ();
set_emin (MPFR_EMIN_MIN);
set_emax (MPFR_EMAX_MAX);
mpfr_init2 (x, GMP_NUMB_BITS);
mpfr_setmax (x, MPFR_EMAX_MAX);
flags1 = MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT;
RND_LOOP (r)
{
/* We also want to test cases where the precision of the
result is twice the precision of each input, in case
add1sp.c and/or sub1sp.c could be involved. */
for (i = 1; i <= 2; i++)
{
mpfr_inits2 (GMP_NUMB_BITS * i, y1, y2, (mpfr_ptr) 0);
inex1 = mpfr_mul (y1, x, x, (mpfr_rnd_t) r);
mpfr_clear_flags ();
inex2 = mpfr_fmma (y2, x, x, x, x, (mpfr_rnd_t) r);
flags2 = __gmpfr_flags;
if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
mpfr_equal_p (y1, y2)))
{
printf ("Error in max_tests for %s\n",
mpfr_print_rnd_mode ((mpfr_rnd_t) r));
printf ("Expected ");
mpfr_dump (y1);
printf (" with inex = %d, flags =", inex1);
flags_out (flags1);
printf ("Got ");
mpfr_dump (y2);
printf (" with inex = %d, flags =", inex2);
flags_out (flags2);
exit (1);
}
mpfr_clears (y1, y2, (mpfr_ptr) 0);
} /* i */
} /* r */
mpfr_clear (x);
set_emin (emin);
set_emax (emax);
}
/* a^2 - (a+k)(a-k) = k^2 where a^2 overflows but k^2 usually doesn't.
* a^2 + cd where a^2 overflows and cd doesn't affect the overflow
* (and cd may even underflow).
*/
static void
overflow_tests (void)
{
mpfr_exp_t old_emax;
int i;
old_emax = mpfr_get_emax ();
for (i = 0; i < 40; i++)
{
mpfr_exp_t emax, exp_a;
mpfr_t a, k, c, d, z1, z2;
mpfr_rnd_t rnd;
mpfr_prec_t prec_a, prec_z;
int add = i & 1, inex, inex1, inex2;
mpfr_flags_t flags1, flags2;
/* In most cases, we do the test with the maximum exponent. */
emax = i % 8 != 0 ? MPFR_EMAX_MAX : 64 + (randlimb () % 1);
set_emax (emax);
exp_a = emax/2 + 32;
rnd = RND_RAND_NO_RNDF ();
prec_a = 64 + randlimb () % 100;
prec_z = MPFR_PREC_MIN + randlimb () % 160;
mpfr_init2 (a, prec_a);
mpfr_urandom (a, RANDS, MPFR_RNDU);
mpfr_set_exp (a, exp_a);
mpfr_init2 (k, prec_a - 32);
mpfr_urandom (k, RANDS, MPFR_RNDU);
mpfr_set_exp (k, exp_a - 32);
mpfr_init2 (c, prec_a + 1);
inex = mpfr_add (c, a, k, MPFR_RNDN);
MPFR_ASSERTN (inex == 0);
mpfr_init2 (d, prec_a);
inex = mpfr_sub (d, a, k, MPFR_RNDN);
MPFR_ASSERTN (inex == 0);
if (add)
mpfr_neg (d, d, MPFR_RNDN);
mpfr_init2 (z1, prec_z);
mpfr_clear_flags ();
inex1 = mpfr_sqr (z1, k, rnd);
flags1 = __gmpfr_flags;
mpfr_init2 (z2, prec_z);
mpfr_clear_flags ();
inex2 = add ?
mpfr_fmma (z2, a, a, c, d, rnd) :
mpfr_fmms (z2, a, a, c, d, rnd);
flags2 = __gmpfr_flags;
if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
mpfr_equal_p (z1, z2)))
{
printf ("Error 1 in overflow_tests for %s\n",
mpfr_print_rnd_mode (rnd));
printf ("Expected ");
mpfr_dump (z1);
printf (" with inex = %d, flags =", inex1);
flags_out (flags1);
printf ("Got ");
mpfr_dump (z2);
printf (" with inex = %d, flags =", inex2);
flags_out (flags2);
exit (1);
}
/* c and d such that a^2 +/- cd ~= a^2 (overflow) */
mpfr_urandom (c, RANDS, MPFR_RNDU);
mpfr_set_exp (c, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin);
if (randlimb () % 1)
mpfr_neg (c, c, MPFR_RNDN);
mpfr_urandom (d, RANDS, MPFR_RNDU);
mpfr_set_exp (d, randlimb () % 1 ? exp_a - 2 : __gmpfr_emin);
if (randlimb () % 1)
mpfr_neg (d, d, MPFR_RNDN);
mpfr_clear_flags ();
inex1 = mpfr_sqr (z1, a, rnd);
flags1 = __gmpfr_flags;
MPFR_ASSERTN (flags1 == (MPFR_FLAGS_OVERFLOW | MPFR_FLAGS_INEXACT));
mpfr_clear_flags ();
inex2 = add ?
mpfr_fmma (z2, a, a, c, d, rnd) :
mpfr_fmms (z2, a, a, c, d, rnd);
flags2 = __gmpfr_flags;
if (! (flags1 == flags2 && SAME_SIGN (inex1, inex2) &&
mpfr_equal_p (z1, z2)))
{
printf ("Error 2 in overflow_tests for %s\n",
mpfr_print_rnd_mode (rnd));
printf ("Expected ");
mpfr_dump (z1);
printf (" with inex = %d, flags =", inex1);
flags_out (flags1);
printf ("Got ");
mpfr_dump (z2);
printf (" with inex = %d, flags =", inex2);
flags_out (flags2);
printf ("a="); mpfr_dump (a);
printf ("c="); mpfr_dump (c);
printf ("d="); mpfr_dump (d);
exit (1);
}
mpfr_clears (a, k, c, d, z1, z2, (mpfr_ptr) 0);
}
set_emax (old_emax);
}
/* (1/2)x + (1/2)x = x tested near underflow */
static void
half_plus_half (void)
{
mpfr_exp_t emin;
mpfr_t h, x1, x2, y;
int neg, r, i, inex;
mpfr_flags_t flags;
emin = mpfr_get_emin ();
set_emin (MPFR_EMIN_MIN);
mpfr_inits2 (4, h, x1, (mpfr_ptr) 0);
mpfr_init2 (x2, GMP_NUMB_BITS);
mpfr_set_ui_2exp (h, 1, -1, MPFR_RNDN);
for (mpfr_setmin (x1, __gmpfr_emin);
MPFR_GET_EXP (x1) < __gmpfr_emin + 2;
mpfr_nextabove (x1))
{
inex = mpfr_set (x2, x1, MPFR_RNDN);
MPFR_ASSERTN (inex == 0);
for (neg = 0; neg < 2; neg++)
{
RND_LOOP (r)
{
/* We also want to test cases where the precision of the
result is twice the precision of each input, in case
add1sp.c and/or sub1sp.c could be involved. */
for (i = 1; i <= 2; i++)
{
mpfr_init2 (y, GMP_NUMB_BITS * i);
mpfr_clear_flags ();
inex = mpfr_fmma (y, h, x2, h, x2, (mpfr_rnd_t) r);
flags = __gmpfr_flags;
if (! (flags == 0 && inex == 0 && mpfr_equal_p (y, x2)))
{
printf ("Error in half_plus_half for %s\n",
mpfr_print_rnd_mode ((mpfr_rnd_t) r));
printf ("Expected ");
mpfr_dump (x2);
printf (" with inex = 0, flags =");
flags_out (0);
printf ("Got ");
mpfr_dump (y);
printf (" with inex = %d, flags =", inex);
flags_out (flags);
exit (1);
}
mpfr_clear (y);
}
}
mpfr_neg (x2, x2, MPFR_RNDN);
}
}
mpfr_clears (h, x1, x2, (mpfr_ptr) 0);
set_emin (emin);
}
/* check that result has exponent in [emin, emax]
(see https://sympa.inria.fr/sympa/arc/mpfr/2017-04/msg00016.html)
Overflow detection in sub1.c was incorrect (only for UBF cases);
fixed in r11414. */
static void
bug20170405 (void)
{
mpfr_t x, y, z;
mpfr_inits2 (866, x, y, z, (mpfr_ptr) 0);
mpfr_set_str_binary (x, "-0.10010101110110001111000110010100011001111011110001010100010000111110001010111110100001000000011010001000010000101110000000001100001011000110010000100111001100000101110000000001001101101101010110000110100010010111011001101101010011111000101100000010001100010000011100000000011110100010111011101011000101101011110110001010011001101110011101100001111000011000000011000010101010000101001001010000111101100001000001011110011000110010001100001101101001001010000111100101000010111001001101010011001110110001000010101001100000101010110000100100100010101011111001001100010001010110011000000001011110011000110001000100101000111010010111111110010111001101110101010010101101010100111001011100101101010111010011001000001010011001010001101000111011010010100110011001111111000011101111001010111001001011011011110101101001100011010001010110011100001101100100001001100111010100010100E768635654");
mpfr_set_str_binary (y, "-0.11010001010111110010110101010011000010010011010011011101100100110000110101100110111010001001110101110000011101100010110100100011001101111010100011111001011100111101110101101001000101011110101101101011010100110010111111010011011100101111110011001001010101011101111100011101100001010010011000110010110101001110010001101111111001100100000101010100110011101101101010011001000110100001001100000010110010101111000110110000111011000110001000100100100101111110001111100101011100100100110111010000010110110001110010001001101000000110111000101000110101111110000110001110100010101111010110001111010111111111010011001001100110011000110010110011000101110001010001101000100010000110011101010010010011110100000111100000101100110001111010000100011111000001101111110100000011011110010100010010011111111000010110000000011010011001100110001110111111010111110000111110010110011001000010E768635576");
/* since emax = 1073741821, x^2-y^2 should overflow */
mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
/* same test when z has a different precision */
mpfr_set_prec (z, 867);
mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
MPFR_ASSERTN(mpfr_inf_p (z) && mpfr_sgn (z) > 0);
mpfr_set_prec (x, 564);
mpfr_set_prec (y, 564);
mpfr_set_prec (z, 2256);
mpfr_set_str_binary (x, "1e-536870913");
mpfr_set_str_binary (y, "-1e-536870913");
mpfr_fmma (z, x, y, x, y, MPFR_RNDN);
/* we should get -0 as result */
MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z));
mpfr_set_prec (x, 564);
mpfr_set_prec (y, 564);
mpfr_set_prec (z, 2256);
mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100e-737194993");
mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010e-737194903");
mpfr_fmma (z, x, y, x, y, MPFR_RNDN);
/* we should get -0 as result */
MPFR_ASSERTN(mpfr_zero_p (z) && mpfr_signbit (z));
mpfr_set_prec (x, 2);
mpfr_set_prec (y, 2);
mpfr_set_prec (z, 2);
/* (a+i*b)*(c+i*d) with:
a=0.10E1
b=0.10E-536870912
c=0.10E-536870912
d=0.10E1 */
mpfr_set_str_binary (x, "0.10E1"); /* x = a = d */
mpfr_set_str_binary (y, "0.10E-536870912"); /* y = b = c */
/* real part is a*c-b*d = x*y-y*x */
mpfr_fmms (z, x, y, y, x, MPFR_RNDN);
MPFR_ASSERTN(mpfr_zero_p (z) && !mpfr_signbit (z));
/* imaginary part is a*d+b*c = x*x+y*y */
mpfr_fmma (z, x, x, y, y, MPFR_RNDN);
MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0);
mpfr_fmma (z, y, y, x, x, MPFR_RNDN);
MPFR_ASSERTN(mpfr_cmp_ui (z, 1) == 0);
mpfr_clears (x, y, z, (mpfr_ptr) 0);
}
static void
underflow_tests (void)
{
mpfr_t x, y, z;
mpfr_prec_t p;
mpfr_exp_t emin;
emin = mpfr_get_emin ();
mpfr_set_emin (-17);
for (p = MPFR_PREC_MIN; p <= 1024; p++)
{
mpfr_inits2 (p, x, y, (mpfr_ptr) 0);
mpfr_init2 (z, p + 1);
mpfr_set_str_binary (x, "1e-10");
mpfr_set_str_binary (y, "1e-11");
mpfr_clear_underflow ();
mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
/* the exact result is 2^-20-2^-22, and should underflow */
MPFR_ASSERTN(mpfr_underflow_p ());
MPFR_ASSERTN(mpfr_zero_p (z));
MPFR_ASSERTN(mpfr_signbit (z) == 0);
mpfr_clears (x, y, z, (mpfr_ptr) 0);
}
mpfr_set_emin (emin);
}
static void
bug20180604 (void)
{
mpfr_t x, y, yneg, z;
mpfr_exp_t emin;
int ret;
emin = mpfr_get_emin ();
mpfr_set_emin (-1073741821);
mpfr_inits2 (564, x, y, yneg, (mpfr_ptr) 0);
mpfr_init2 (z, 2256);
mpfr_set_str_binary (x, "1.10010000111100110011001101111111101000111001011000110100110010000101000100010001000000111100010000101001011011111001111000110101111100101111001100001100011101100100011110000000011000010110111100111000100101010001011111010111011001110010001011101111001011001110110000010000011100010001010001011100100110111110101001001111001011101111110011101110101010110100010010111011111100010101111100011110111001011111101110101101101110100101111010000101011110100000000110111101000001100001000100010110100111010011011010110011100111010000101110010101111001011100110101100001100E-737194993");
mpfr_set_str_binary (y, "-1.00101000100001001101011110100010110011101010011011010111100110101011111100000100101000111010111101100100110010001110011011100100110110000001011001000111101111101111110101100110111000000011000001101001010100010010001110001000011010000100111001001100101111111100010101110101001101101101111010100011011110001000010000010100011000011000010110101100000111111110111001100100100001101011111011100101110111000100101010110100010011101010110010100110100111000000100111101101101000000011110000100110100100011000010011110010001010000110100011111101101101110001110001101101010E-737194903");
mpfr_clear_underflow ();
ret = mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
MPFR_ASSERTN(mpfr_underflow_p ());
MPFR_ASSERTN(mpfr_zero_p (z));
MPFR_ASSERTN(mpfr_signbit (z) == 1);
MPFR_ASSERTN(ret > 0);
mpfr_clear_underflow ();
ret = mpfr_fmms (z, y, y, x, x, MPFR_RNDN);
MPFR_ASSERTN(mpfr_underflow_p ());
MPFR_ASSERTN(mpfr_zero_p (z));
MPFR_ASSERTN(mpfr_signbit (z) == 0);
MPFR_ASSERTN(ret < 0);
mpfr_neg (yneg, y, MPFR_RNDN);
mpfr_clear_underflow ();
ret = mpfr_fmms (z, x, x, y, yneg, MPFR_RNDN);
MPFR_ASSERTN(mpfr_underflow_p ());
MPFR_ASSERTN(mpfr_zero_p (z));
MPFR_ASSERTN(mpfr_signbit (z) == 0);
MPFR_ASSERTN(ret < 0);
mpfr_clear_underflow ();
ret = mpfr_fmms (z, y, yneg, x, x, MPFR_RNDN);
MPFR_ASSERTN(mpfr_underflow_p ());
MPFR_ASSERTN(mpfr_zero_p (z));
MPFR_ASSERTN(mpfr_signbit (z) == 1);
MPFR_ASSERTN(ret > 0);
mpfr_clears (x, y, yneg, z, (mpfr_ptr) 0);
mpfr_set_emin (emin);
}
/* this bug was discovered from mpc_mul */
static void
bug20200206 (void)
{
mpfr_exp_t emin = mpfr_get_emin ();
mpfr_t xre, xim, yre, yim, zre;
mpfr_set_emin (-1073);
mpfr_inits2 (53, xre, xim, yre, yim, zre, (mpfr_ptr) 0);
mpfr_set_str (xre, "-0x8.294611b331c8p-904", 16, MPFR_RNDN);
mpfr_set_str (xim, "-0x1.859278c2992fap-676", 16, MPFR_RNDN);
mpfr_set_str (yre, "0x9.ac54802a95f8p-820", 16, MPFR_RNDN);
mpfr_set_str (yim, "0x3.17e59e7612aap-412", 16, MPFR_RNDN);
mpfr_fmms (zre, xre, yre, xim, yim, MPFR_RNDN);
if (mpfr_regular_p (zre) && mpfr_get_exp (zre) < -1073)
{
printf ("Error, mpfr_fmms returns an out-of-range exponent:\n");
mpfr_dump (zre);
exit (1);
}
mpfr_clears (xre, xim, yre, yim, zre, (mpfr_ptr) 0);
mpfr_set_emin (emin);
}
/* check for integer overflow (see bug fixed in r13798) */
static void
extreme_underflow (void)
{
mpfr_t x, y, z;
mpfr_exp_t emin = mpfr_get_emin ();
set_emin (MPFR_EMIN_MIN);
mpfr_inits2 (64, x, y, z, (mpfr_ptr) 0);
mpfr_set_ui_2exp (x, 1, MPFR_EMIN_MIN - 1, MPFR_RNDN);
mpfr_set (y, x, MPFR_RNDN);
mpfr_nextabove (x);
mpfr_clear_flags ();
mpfr_fmms (z, x, x, y, y, MPFR_RNDN);
MPFR_ASSERTN (__gmpfr_flags == (MPFR_FLAGS_UNDERFLOW | MPFR_FLAGS_INEXACT));
MPFR_ASSERTN (MPFR_IS_ZERO (z) && MPFR_IS_POS (z));
mpfr_clears (x, y, z, (mpfr_ptr) 0);
set_emin (emin);
}
/* Test double-rounding cases in mpfr_set_1_2, which is called when
all the precisions are the same. With set.c before r13347, this
triggers errors for neg=0. */
static void
double_rounding (void)
{
int p;
for (p = 4; p < 4 * GMP_NUMB_BITS; p++)
{
mpfr_t a, b, c, d, z1, z2;
int neg;
mpfr_inits2 (p, a, b, c, d, z1, z2, (mpfr_ptr) 0);
/* Choose a, b, c, d such that a*b = 2^p and c*d = 1 + epsilon */
mpfr_set_ui (a, 2, MPFR_RNDN);
mpfr_set_ui_2exp (b, 1, p - 1, MPFR_RNDN);
mpfr_set_ui (c, 1, MPFR_RNDN);
mpfr_nextabove (c);
/* c = 1 + ulp(1) = 1 + 2 * ulp(1/2) */
for (neg = 0; neg <= 1; neg++)
{
int inex1, inex2;
/* Set d = 1 - (1 + neg) * ulp(1/2), thus
* c*d = 1 + (1 - neg) * ulp(1/2) - 2 * (1 + neg) * ulp(1/2)^2,
* so that a*b + c*d rounds to 2^p + 1 and epsilon has the
* same sign as (-1)^neg.
*/
mpfr_set_ui (d, 1, MPFR_RNDN);
mpfr_nextbelow (d);
mpfr_set_ui_2exp (z1, 1, p, MPFR_RNDN);
if (neg)
{
mpfr_nextbelow (d);
inex1 = -1;
}
else
{
mpfr_nextabove (z1);
inex1 = 1;
}
inex2 = mpfr_fmma (z2, a, b, c, d, MPFR_RNDN);
if (!(mpfr_equal_p (z1, z2) && SAME_SIGN (inex1, inex2)))
{
printf ("Error in double_rounding for precision %d, neg=%d\n",
p, neg);
printf ("Expected ");
mpfr_dump (z1);
printf ("with inex = %d\n", inex1);
printf ("Got ");
mpfr_dump (z2);
printf ("with inex = %d\n", inex2);
exit (1);
}
}
mpfr_clears (a, b, c, d, z1, z2, (mpfr_ptr) 0);
}
}
int
main (int argc, char *argv[])
{
tests_start_mpfr ();
bug20200206 ();
bug20180604 ();
underflow_tests ();
random_tests ();
zero_tests ();
max_tests ();
overflow_tests ();
half_plus_half ();
bug20170405 ();
double_rounding ();
extreme_underflow ();
tests_end_mpfr ();
return 0;
}