Home | History | Annotate | Line # | Download | only in tests
tmul.c revision 1.1.1.4
      1      1.1  mrg /* Test file for mpfr_mul.
      2      1.1  mrg 
      3  1.1.1.4  mrg Copyright 1999, 2001-2018 Free Software Foundation, Inc.
      4  1.1.1.3  mrg Contributed by the AriC and Caramba projects, INRIA.
      5      1.1  mrg 
      6      1.1  mrg This file is part of the GNU MPFR Library.
      7      1.1  mrg 
      8      1.1  mrg The GNU MPFR Library is free software; you can redistribute it and/or modify
      9      1.1  mrg it under the terms of the GNU Lesser General Public License as published by
     10      1.1  mrg the Free Software Foundation; either version 3 of the License, or (at your
     11      1.1  mrg option) any later version.
     12      1.1  mrg 
     13      1.1  mrg The GNU MPFR Library is distributed in the hope that it will be useful, but
     14      1.1  mrg WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
     15      1.1  mrg or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
     16      1.1  mrg License for more details.
     17      1.1  mrg 
     18      1.1  mrg You should have received a copy of the GNU Lesser General Public License
     19      1.1  mrg along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
     20      1.1  mrg http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
     21      1.1  mrg 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
     22      1.1  mrg 
     23      1.1  mrg #include "mpfr-test.h"
     24      1.1  mrg 
     25      1.1  mrg #ifdef CHECK_EXTERNAL
     26      1.1  mrg static int
     27      1.1  mrg test_mul (mpfr_ptr a, mpfr_srcptr b, mpfr_srcptr c, mpfr_rnd_t rnd_mode)
     28      1.1  mrg {
     29      1.1  mrg   int res;
     30      1.1  mrg   int ok = rnd_mode == MPFR_RNDN && mpfr_number_p (b) && mpfr_number_p (c);
     31      1.1  mrg   if (ok)
     32      1.1  mrg     {
     33      1.1  mrg       mpfr_print_raw (b);
     34      1.1  mrg       printf (" ");
     35      1.1  mrg       mpfr_print_raw (c);
     36      1.1  mrg     }
     37      1.1  mrg   res = mpfr_mul (a, b, c, rnd_mode);
     38      1.1  mrg   if (ok)
     39      1.1  mrg     {
     40      1.1  mrg       printf (" ");
     41      1.1  mrg       mpfr_print_raw (a);
     42      1.1  mrg       printf ("\n");
     43      1.1  mrg     }
     44      1.1  mrg   return res;
     45      1.1  mrg }
     46      1.1  mrg #else
     47      1.1  mrg #define test_mul mpfr_mul
     48      1.1  mrg #endif
     49      1.1  mrg 
     50      1.1  mrg /* checks that xs * ys gives the expected result res */
     51      1.1  mrg static void
     52      1.1  mrg check (const char *xs, const char *ys, mpfr_rnd_t rnd_mode,
     53      1.1  mrg         unsigned int px, unsigned int py, unsigned int pz, const char *res)
     54      1.1  mrg {
     55      1.1  mrg   mpfr_t xx, yy, zz;
     56      1.1  mrg 
     57      1.1  mrg   mpfr_init2 (xx, px);
     58      1.1  mrg   mpfr_init2 (yy, py);
     59      1.1  mrg   mpfr_init2 (zz, pz);
     60      1.1  mrg   mpfr_set_str1 (xx, xs);
     61      1.1  mrg   mpfr_set_str1 (yy, ys);
     62      1.1  mrg   test_mul(zz, xx, yy, rnd_mode);
     63      1.1  mrg   if (mpfr_cmp_str1 (zz, res) )
     64      1.1  mrg     {
     65  1.1.1.4  mrg       printf ("(1) mpfr_mul failed for x=%s y=%s with rnd=%s\n",
     66      1.1  mrg               xs, ys, mpfr_print_rnd_mode (rnd_mode));
     67      1.1  mrg       printf ("correct is %s, mpfr_mul gives ", res);
     68  1.1.1.4  mrg       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
     69  1.1.1.4  mrg       putchar ('\n');
     70      1.1  mrg       exit (1);
     71      1.1  mrg     }
     72  1.1.1.4  mrg   mpfr_clear (xx);
     73  1.1.1.4  mrg   mpfr_clear (yy);
     74  1.1.1.4  mrg   mpfr_clear (zz);
     75      1.1  mrg }
     76      1.1  mrg 
     77      1.1  mrg static void
     78      1.1  mrg check53 (const char *xs, const char *ys, mpfr_rnd_t rnd_mode, const char *zs)
     79      1.1  mrg {
     80      1.1  mrg   mpfr_t xx, yy, zz;
     81      1.1  mrg 
     82      1.1  mrg   mpfr_inits2 (53, xx, yy, zz, (mpfr_ptr) 0);
     83      1.1  mrg   mpfr_set_str1 (xx, xs);
     84      1.1  mrg   mpfr_set_str1 (yy, ys);
     85      1.1  mrg   test_mul (zz, xx, yy, rnd_mode);
     86      1.1  mrg   if (mpfr_cmp_str1 (zz, zs) )
     87      1.1  mrg     {
     88      1.1  mrg       printf ("(2) mpfr_mul failed for x=%s y=%s with rnd=%s\n",
     89      1.1  mrg               xs, ys, mpfr_print_rnd_mode(rnd_mode));
     90      1.1  mrg       printf ("correct result is %s,\n mpfr_mul gives ", zs);
     91  1.1.1.4  mrg       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
     92  1.1.1.4  mrg       putchar ('\n');
     93      1.1  mrg       exit (1);
     94      1.1  mrg     }
     95      1.1  mrg   mpfr_clears (xx, yy, zz, (mpfr_ptr) 0);
     96      1.1  mrg }
     97      1.1  mrg 
     98      1.1  mrg /* checks that x*y gives the right result with 24 bits of precision */
     99      1.1  mrg static void
    100      1.1  mrg check24 (const char *xs, const char *ys, mpfr_rnd_t rnd_mode, const char *zs)
    101      1.1  mrg {
    102      1.1  mrg   mpfr_t xx, yy, zz;
    103      1.1  mrg 
    104      1.1  mrg   mpfr_inits2 (24, xx, yy, zz, (mpfr_ptr) 0);
    105      1.1  mrg   mpfr_set_str1 (xx, xs);
    106      1.1  mrg   mpfr_set_str1 (yy, ys);
    107      1.1  mrg   test_mul (zz, xx, yy, rnd_mode);
    108      1.1  mrg   if (mpfr_cmp_str1 (zz, zs) )
    109      1.1  mrg     {
    110      1.1  mrg       printf ("(3) mpfr_mul failed for x=%s y=%s with "
    111      1.1  mrg               "rnd=%s\n", xs, ys, mpfr_print_rnd_mode(rnd_mode));
    112      1.1  mrg       printf ("correct result is gives %s, mpfr_mul gives ", zs);
    113  1.1.1.4  mrg       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
    114  1.1.1.4  mrg       putchar ('\n');
    115      1.1  mrg       exit (1);
    116      1.1  mrg     }
    117      1.1  mrg   mpfr_clears (xx, yy, zz, (mpfr_ptr) 0);
    118      1.1  mrg }
    119      1.1  mrg 
    120      1.1  mrg /* the following examples come from the paper "Number-theoretic Test
    121      1.1  mrg    Generation for Directed Rounding" from Michael Parks, Table 1 */
    122      1.1  mrg static void
    123      1.1  mrg check_float (void)
    124      1.1  mrg {
    125      1.1  mrg   check24("8388609.0",  "8388609.0", MPFR_RNDN, "70368760954880.0");
    126      1.1  mrg   check24("16777213.0", "8388609.0", MPFR_RNDN, "140737479966720.0");
    127      1.1  mrg   check24("8388611.0",  "8388609.0", MPFR_RNDN, "70368777732096.0");
    128      1.1  mrg   check24("12582911.0", "8388610.0", MPFR_RNDN, "105553133043712.0");
    129      1.1  mrg   check24("12582914.0", "8388610.0", MPFR_RNDN, "105553158209536.0");
    130      1.1  mrg   check24("13981013.0", "8388611.0", MPFR_RNDN, "117281279442944.0");
    131      1.1  mrg   check24("11184811.0", "8388611.0", MPFR_RNDN, "93825028587520.0");
    132      1.1  mrg   check24("11184810.0", "8388611.0", MPFR_RNDN, "93825020198912.0");
    133      1.1  mrg   check24("13981014.0", "8388611.0", MPFR_RNDN, "117281287831552.0");
    134      1.1  mrg 
    135      1.1  mrg   check24("8388609.0",  "8388609.0", MPFR_RNDZ, "70368760954880.0");
    136      1.1  mrg   check24("16777213.0", "8388609.0", MPFR_RNDZ, "140737471578112.0");
    137      1.1  mrg   check24("8388611.0",  "8388609.0", MPFR_RNDZ, "70368777732096.0");
    138      1.1  mrg   check24("12582911.0", "8388610.0", MPFR_RNDZ, "105553124655104.0");
    139      1.1  mrg   check24("12582914.0", "8388610.0", MPFR_RNDZ, "105553158209536.0");
    140      1.1  mrg   check24("13981013.0", "8388611.0", MPFR_RNDZ, "117281271054336.0");
    141      1.1  mrg   check24("11184811.0", "8388611.0", MPFR_RNDZ, "93825028587520.0");
    142      1.1  mrg   check24("11184810.0", "8388611.0", MPFR_RNDZ, "93825011810304.0");
    143      1.1  mrg   check24("13981014.0", "8388611.0", MPFR_RNDZ, "117281287831552.0");
    144      1.1  mrg 
    145      1.1  mrg   check24("8388609.0",  "8388609.0", MPFR_RNDU, "70368769343488.0");
    146      1.1  mrg   check24("16777213.0", "8388609.0", MPFR_RNDU, "140737479966720.0");
    147      1.1  mrg   check24("8388611.0",  "8388609.0", MPFR_RNDU, "70368786120704.0");
    148      1.1  mrg   check24("12582911.0", "8388610.0", MPFR_RNDU, "105553133043712.0");
    149      1.1  mrg   check24("12582914.0", "8388610.0", MPFR_RNDU, "105553166598144.0");
    150      1.1  mrg   check24("13981013.0", "8388611.0", MPFR_RNDU, "117281279442944.0");
    151      1.1  mrg   check24("11184811.0", "8388611.0", MPFR_RNDU, "93825036976128.0");
    152      1.1  mrg   check24("11184810.0", "8388611.0", MPFR_RNDU, "93825020198912.0");
    153      1.1  mrg   check24("13981014.0", "8388611.0", MPFR_RNDU, "117281296220160.0");
    154      1.1  mrg 
    155      1.1  mrg   check24("8388609.0",  "8388609.0", MPFR_RNDD, "70368760954880.0");
    156      1.1  mrg   check24("16777213.0", "8388609.0", MPFR_RNDD, "140737471578112.0");
    157      1.1  mrg   check24("8388611.0",  "8388609.0", MPFR_RNDD, "70368777732096.0");
    158      1.1  mrg   check24("12582911.0", "8388610.0", MPFR_RNDD, "105553124655104.0");
    159      1.1  mrg   check24("12582914.0", "8388610.0", MPFR_RNDD, "105553158209536.0");
    160      1.1  mrg   check24("13981013.0", "8388611.0", MPFR_RNDD, "117281271054336.0");
    161      1.1  mrg   check24("11184811.0", "8388611.0", MPFR_RNDD, "93825028587520.0");
    162      1.1  mrg   check24("11184810.0", "8388611.0", MPFR_RNDD, "93825011810304.0");
    163      1.1  mrg   check24("13981014.0", "8388611.0", MPFR_RNDD, "117281287831552.0");
    164      1.1  mrg }
    165      1.1  mrg 
    166      1.1  mrg /* check sign of result */
    167      1.1  mrg static void
    168      1.1  mrg check_sign (void)
    169      1.1  mrg {
    170      1.1  mrg   mpfr_t a, b;
    171      1.1  mrg 
    172      1.1  mrg   mpfr_init2 (a, 53);
    173      1.1  mrg   mpfr_init2 (b, 53);
    174      1.1  mrg   mpfr_set_si (a, -1, MPFR_RNDN);
    175      1.1  mrg   mpfr_set_ui (b, 2, MPFR_RNDN);
    176      1.1  mrg   test_mul(a, b, b, MPFR_RNDN);
    177      1.1  mrg   if (mpfr_cmp_ui (a, 4) )
    178      1.1  mrg     {
    179      1.1  mrg       printf ("2.0*2.0 gives \n");
    180  1.1.1.4  mrg       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
    181  1.1.1.4  mrg       putchar ('\n');
    182      1.1  mrg       exit (1);
    183      1.1  mrg     }
    184      1.1  mrg   mpfr_clear(a); mpfr_clear(b);
    185      1.1  mrg }
    186      1.1  mrg 
    187      1.1  mrg /* checks that the inexact return value is correct */
    188      1.1  mrg static void
    189      1.1  mrg check_exact (void)
    190      1.1  mrg {
    191      1.1  mrg   mpfr_t a, b, c, d;
    192      1.1  mrg   mpfr_prec_t prec;
    193      1.1  mrg   int i, inexact;
    194      1.1  mrg   mpfr_rnd_t rnd;
    195      1.1  mrg 
    196      1.1  mrg   mpfr_init (a);
    197      1.1  mrg   mpfr_init (b);
    198      1.1  mrg   mpfr_init (c);
    199      1.1  mrg   mpfr_init (d);
    200      1.1  mrg 
    201      1.1  mrg   mpfr_set_prec (a, 17);
    202      1.1  mrg   mpfr_set_prec (b, 17);
    203      1.1  mrg   mpfr_set_prec (c, 32);
    204      1.1  mrg   mpfr_set_str_binary (a, "1.1000111011000100e-1");
    205      1.1  mrg   mpfr_set_str_binary (b, "1.0010001111100111e-1");
    206      1.1  mrg   if (test_mul (c, a, b, MPFR_RNDZ))
    207      1.1  mrg     {
    208      1.1  mrg       printf ("wrong return value (1)\n");
    209      1.1  mrg       exit (1);
    210      1.1  mrg     }
    211      1.1  mrg 
    212  1.1.1.4  mrg   for (prec = MPFR_PREC_MIN; prec < 100; prec++)
    213      1.1  mrg     {
    214      1.1  mrg       mpfr_set_prec (a, prec);
    215      1.1  mrg       mpfr_set_prec (b, prec);
    216  1.1.1.4  mrg       /* for prec=1, ensure PREC(c) >= 1 */
    217  1.1.1.4  mrg       mpfr_set_prec (c, 2 * prec - 2 + (prec == 1));
    218      1.1  mrg       mpfr_set_prec (d, 2 * prec);
    219      1.1  mrg       for (i = 0; i < 1000; i++)
    220      1.1  mrg         {
    221      1.1  mrg           mpfr_urandomb (a, RANDS);
    222      1.1  mrg           mpfr_urandomb (b, RANDS);
    223      1.1  mrg           rnd = RND_RAND ();
    224      1.1  mrg           inexact = test_mul (c, a, b, rnd);
    225      1.1  mrg           if (test_mul (d, a, b, rnd)) /* should be always exact */
    226      1.1  mrg             {
    227      1.1  mrg               printf ("unexpected inexact return value\n");
    228      1.1  mrg               exit (1);
    229      1.1  mrg             }
    230  1.1.1.4  mrg           if ((inexact == 0) && mpfr_cmp (c, d) && rnd != MPFR_RNDF)
    231      1.1  mrg             {
    232  1.1.1.4  mrg               printf ("rnd=%s: inexact=0 but results differ\n",
    233  1.1.1.4  mrg                       mpfr_print_rnd_mode (rnd));
    234  1.1.1.4  mrg               printf ("a=");
    235  1.1.1.4  mrg               mpfr_out_str (stdout, 2, 0, a, rnd);
    236  1.1.1.4  mrg               printf ("\nb=");
    237  1.1.1.4  mrg               mpfr_out_str (stdout, 2, 0, b, rnd);
    238  1.1.1.4  mrg               printf ("\nc=");
    239  1.1.1.4  mrg               mpfr_out_str (stdout, 2, 0, c, rnd);
    240  1.1.1.4  mrg               printf ("\nd=");
    241  1.1.1.4  mrg               mpfr_out_str (stdout, 2, 0, d, rnd);
    242  1.1.1.4  mrg               printf ("\n");
    243      1.1  mrg               exit (1);
    244      1.1  mrg             }
    245  1.1.1.4  mrg           else if (inexact && (mpfr_cmp (c, d) == 0) && rnd != MPFR_RNDF)
    246      1.1  mrg             {
    247      1.1  mrg               printf ("inexact!=0 but results agree\n");
    248      1.1  mrg               printf ("prec=%u rnd=%s a=", (unsigned int) prec,
    249      1.1  mrg                       mpfr_print_rnd_mode (rnd));
    250      1.1  mrg               mpfr_out_str (stdout, 2, 0, a, rnd);
    251      1.1  mrg               printf ("\nb=");
    252      1.1  mrg               mpfr_out_str (stdout, 2, 0, b, rnd);
    253      1.1  mrg               printf ("\nc=");
    254      1.1  mrg               mpfr_out_str (stdout, 2, 0, c, rnd);
    255      1.1  mrg               printf ("\nd=");
    256      1.1  mrg               mpfr_out_str (stdout, 2, 0, d, rnd);
    257      1.1  mrg               printf ("\n");
    258      1.1  mrg               exit (1);
    259      1.1  mrg             }
    260      1.1  mrg         }
    261      1.1  mrg     }
    262      1.1  mrg 
    263      1.1  mrg   mpfr_clear (a);
    264      1.1  mrg   mpfr_clear (b);
    265      1.1  mrg   mpfr_clear (c);
    266      1.1  mrg   mpfr_clear (d);
    267      1.1  mrg }
    268      1.1  mrg 
    269      1.1  mrg static void
    270      1.1  mrg check_max(void)
    271      1.1  mrg {
    272      1.1  mrg   mpfr_t xx, yy, zz;
    273      1.1  mrg   mpfr_exp_t emin;
    274      1.1  mrg 
    275      1.1  mrg   mpfr_init2(xx, 4);
    276      1.1  mrg   mpfr_init2(yy, 4);
    277      1.1  mrg   mpfr_init2(zz, 4);
    278      1.1  mrg   mpfr_set_str1 (xx, "0.68750");
    279      1.1  mrg   mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT/2, MPFR_RNDN);
    280      1.1  mrg   mpfr_set_str1 (yy, "0.68750");
    281      1.1  mrg   mpfr_mul_2si(yy, yy, MPFR_EMAX_DEFAULT - MPFR_EMAX_DEFAULT/2 + 1, MPFR_RNDN);
    282      1.1  mrg   mpfr_clear_flags();
    283      1.1  mrg   test_mul(zz, xx, yy, MPFR_RNDU);
    284      1.1  mrg   if (!(mpfr_overflow_p() && MPFR_IS_INF(zz)))
    285      1.1  mrg     {
    286      1.1  mrg       printf("check_max failed (should be an overflow)\n");
    287      1.1  mrg       exit(1);
    288      1.1  mrg     }
    289      1.1  mrg 
    290      1.1  mrg   mpfr_clear_flags();
    291      1.1  mrg   test_mul(zz, xx, yy, MPFR_RNDD);
    292      1.1  mrg   if (mpfr_overflow_p() || MPFR_IS_INF(zz))
    293      1.1  mrg     {
    294      1.1  mrg       printf("check_max failed (should NOT be an overflow)\n");
    295      1.1  mrg       exit(1);
    296      1.1  mrg     }
    297      1.1  mrg   mpfr_set_str1 (xx, "0.93750");
    298      1.1  mrg   mpfr_mul_2si(xx, xx, MPFR_EMAX_DEFAULT, MPFR_RNDN);
    299      1.1  mrg   if (!(MPFR_IS_FP(xx) && MPFR_IS_FP(zz)))
    300      1.1  mrg     {
    301      1.1  mrg       printf("check_max failed (internal error)\n");
    302      1.1  mrg       exit(1);
    303      1.1  mrg     }
    304      1.1  mrg   if (mpfr_cmp(xx, zz) != 0)
    305      1.1  mrg     {
    306      1.1  mrg       printf("check_max failed: got ");
    307      1.1  mrg       mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
    308      1.1  mrg       printf(" instead of ");
    309      1.1  mrg       mpfr_out_str(stdout, 2, 0, xx, MPFR_RNDZ);
    310      1.1  mrg       printf("\n");
    311      1.1  mrg       exit(1);
    312      1.1  mrg     }
    313      1.1  mrg 
    314      1.1  mrg   /* check underflow */
    315      1.1  mrg   emin = mpfr_get_emin ();
    316      1.1  mrg   set_emin (0);
    317      1.1  mrg   mpfr_set_str_binary (xx, "0.1E0");
    318      1.1  mrg   mpfr_set_str_binary (yy, "0.1E0");
    319      1.1  mrg   test_mul (zz, xx, yy, MPFR_RNDN);
    320      1.1  mrg   /* exact result is 0.1E-1, which should round to 0 */
    321      1.1  mrg   MPFR_ASSERTN(mpfr_cmp_ui (zz, 0) == 0 && MPFR_IS_POS(zz));
    322      1.1  mrg   set_emin (emin);
    323      1.1  mrg 
    324      1.1  mrg   /* coverage test for mpfr_powerof2_raw */
    325      1.1  mrg   emin = mpfr_get_emin ();
    326      1.1  mrg   set_emin (0);
    327      1.1  mrg   mpfr_set_prec (xx, mp_bits_per_limb + 1);
    328      1.1  mrg   mpfr_set_str_binary (xx, "0.1E0");
    329      1.1  mrg   mpfr_nextabove (xx);
    330      1.1  mrg   mpfr_set_str_binary (yy, "0.1E0");
    331      1.1  mrg   test_mul (zz, xx, yy, MPFR_RNDN);
    332      1.1  mrg   /* exact result is just above 0.1E-1, which should round to minfloat */
    333      1.1  mrg   MPFR_ASSERTN(mpfr_cmp (zz, yy) == 0);
    334      1.1  mrg   set_emin (emin);
    335      1.1  mrg 
    336      1.1  mrg   mpfr_clear(xx);
    337      1.1  mrg   mpfr_clear(yy);
    338      1.1  mrg   mpfr_clear(zz);
    339      1.1  mrg }
    340      1.1  mrg 
    341      1.1  mrg static void
    342      1.1  mrg check_min(void)
    343      1.1  mrg {
    344      1.1  mrg   mpfr_t xx, yy, zz;
    345      1.1  mrg 
    346      1.1  mrg   mpfr_init2(xx, 4);
    347      1.1  mrg   mpfr_init2(yy, 4);
    348      1.1  mrg   mpfr_init2(zz, 3);
    349      1.1  mrg   mpfr_set_str1(xx, "0.9375");
    350      1.1  mrg   mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT/2, MPFR_RNDN);
    351      1.1  mrg   mpfr_set_str1(yy, "0.9375");
    352      1.1  mrg   mpfr_mul_2si(yy, yy, MPFR_EMIN_DEFAULT - MPFR_EMIN_DEFAULT/2 - 1, MPFR_RNDN);
    353      1.1  mrg   test_mul(zz, xx, yy, MPFR_RNDD);
    354  1.1.1.4  mrg   if (MPFR_NOTZERO (zz))
    355      1.1  mrg     {
    356      1.1  mrg       printf("check_min failed: got ");
    357      1.1  mrg       mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
    358      1.1  mrg       printf(" instead of 0\n");
    359      1.1  mrg       exit(1);
    360      1.1  mrg     }
    361      1.1  mrg 
    362      1.1  mrg   test_mul(zz, xx, yy, MPFR_RNDU);
    363      1.1  mrg   mpfr_set_str1 (xx, "0.5");
    364      1.1  mrg   mpfr_mul_2si(xx, xx, MPFR_EMIN_DEFAULT, MPFR_RNDN);
    365      1.1  mrg   if (mpfr_sgn(xx) <= 0)
    366      1.1  mrg     {
    367      1.1  mrg       printf("check_min failed (internal error)\n");
    368      1.1  mrg       exit(1);
    369      1.1  mrg     }
    370      1.1  mrg   if (mpfr_cmp(xx, zz) != 0)
    371      1.1  mrg     {
    372      1.1  mrg       printf("check_min failed: got ");
    373      1.1  mrg       mpfr_out_str(stdout, 2, 0, zz, MPFR_RNDZ);
    374      1.1  mrg       printf(" instead of ");
    375      1.1  mrg       mpfr_out_str(stdout, 2, 0, xx, MPFR_RNDZ);
    376      1.1  mrg       printf("\n");
    377      1.1  mrg       exit(1);
    378      1.1  mrg     }
    379      1.1  mrg 
    380      1.1  mrg   mpfr_clear(xx);
    381      1.1  mrg   mpfr_clear(yy);
    382      1.1  mrg   mpfr_clear(zz);
    383      1.1  mrg }
    384      1.1  mrg 
    385      1.1  mrg static void
    386      1.1  mrg check_nans (void)
    387      1.1  mrg {
    388      1.1  mrg   mpfr_t  p, x, y;
    389      1.1  mrg 
    390      1.1  mrg   mpfr_init2 (x, 123L);
    391      1.1  mrg   mpfr_init2 (y, 123L);
    392      1.1  mrg   mpfr_init2 (p, 123L);
    393      1.1  mrg 
    394      1.1  mrg   /* nan * 0 == nan */
    395      1.1  mrg   mpfr_set_nan (x);
    396      1.1  mrg   mpfr_set_ui (y, 0L, MPFR_RNDN);
    397      1.1  mrg   test_mul (p, x, y, MPFR_RNDN);
    398      1.1  mrg   MPFR_ASSERTN (mpfr_nan_p (p));
    399      1.1  mrg 
    400      1.1  mrg   /* 1 * nan == nan */
    401      1.1  mrg   mpfr_set_ui (x, 1L, MPFR_RNDN);
    402      1.1  mrg   mpfr_set_nan (y);
    403      1.1  mrg   test_mul (p, x, y, MPFR_RNDN);
    404      1.1  mrg   MPFR_ASSERTN (mpfr_nan_p (p));
    405      1.1  mrg 
    406      1.1  mrg   /* 0 * +inf == nan */
    407      1.1  mrg   mpfr_set_ui (x, 0L, MPFR_RNDN);
    408      1.1  mrg   mpfr_set_nan (y);
    409      1.1  mrg   test_mul (p, x, y, MPFR_RNDN);
    410      1.1  mrg   MPFR_ASSERTN (mpfr_nan_p (p));
    411      1.1  mrg 
    412      1.1  mrg   /* +1 * +inf == +inf */
    413      1.1  mrg   mpfr_set_ui (x, 1L, MPFR_RNDN);
    414      1.1  mrg   mpfr_set_inf (y, 1);
    415      1.1  mrg   test_mul (p, x, y, MPFR_RNDN);
    416      1.1  mrg   MPFR_ASSERTN (mpfr_inf_p (p));
    417      1.1  mrg   MPFR_ASSERTN (mpfr_sgn (p) > 0);
    418      1.1  mrg 
    419      1.1  mrg   /* -1 * +inf == -inf */
    420      1.1  mrg   mpfr_set_si (x, -1L, MPFR_RNDN);
    421      1.1  mrg   mpfr_set_inf (y, 1);
    422      1.1  mrg   test_mul (p, x, y, MPFR_RNDN);
    423      1.1  mrg   MPFR_ASSERTN (mpfr_inf_p (p));
    424      1.1  mrg   MPFR_ASSERTN (mpfr_sgn (p) < 0);
    425      1.1  mrg 
    426      1.1  mrg   mpfr_clear (x);
    427      1.1  mrg   mpfr_clear (y);
    428      1.1  mrg   mpfr_clear (p);
    429      1.1  mrg }
    430      1.1  mrg 
    431      1.1  mrg #define BUFSIZE 1552
    432      1.1  mrg 
    433      1.1  mrg static void
    434      1.1  mrg get_string (char *s, FILE *fp)
    435      1.1  mrg {
    436      1.1  mrg   int c, n = BUFSIZE;
    437      1.1  mrg 
    438      1.1  mrg   while ((c = getc (fp)) != '\n')
    439      1.1  mrg     {
    440      1.1  mrg       if (c == EOF)
    441      1.1  mrg         {
    442      1.1  mrg           printf ("Error in get_string: end of file\n");
    443      1.1  mrg           exit (1);
    444      1.1  mrg         }
    445      1.1  mrg       *(unsigned char *)s++ = c;
    446      1.1  mrg       if (--n == 0)
    447      1.1  mrg         {
    448      1.1  mrg           printf ("Error in get_string: buffer is too small\n");
    449      1.1  mrg           exit (1);
    450      1.1  mrg         }
    451      1.1  mrg     }
    452      1.1  mrg   *s = '\0';
    453      1.1  mrg }
    454      1.1  mrg 
    455      1.1  mrg static void
    456      1.1  mrg check_regression (void)
    457      1.1  mrg {
    458      1.1  mrg   mpfr_t x, y, z;
    459      1.1  mrg   int i;
    460      1.1  mrg   FILE *fp;
    461      1.1  mrg   char s[BUFSIZE];
    462      1.1  mrg 
    463      1.1  mrg   mpfr_inits2 (6177, x, y, z, (mpfr_ptr) 0);
    464      1.1  mrg   /* we read long strings from a file since ISO C90 does not support strings of
    465      1.1  mrg      length > 509 */
    466      1.1  mrg   fp = src_fopen ("tmul.dat", "r");
    467      1.1  mrg   if (fp == NULL)
    468      1.1  mrg     {
    469      1.1  mrg       fprintf (stderr, "Error, cannot open tmul.dat in srcdir\n");
    470      1.1  mrg       exit (1);
    471      1.1  mrg     }
    472      1.1  mrg   get_string (s, fp);
    473      1.1  mrg   mpfr_set_str (y, s, 16, MPFR_RNDN);
    474      1.1  mrg   get_string (s, fp);
    475      1.1  mrg   mpfr_set_str (z, s, 16, MPFR_RNDN);
    476      1.1  mrg   i = mpfr_mul (x, y, z, MPFR_RNDN);
    477      1.1  mrg   get_string (s, fp);
    478      1.1  mrg   if (mpfr_cmp_str (x, s, 16, MPFR_RNDN) != 0 || i != -1)
    479      1.1  mrg     {
    480      1.1  mrg       printf ("Regression test 1 failed (i=%d, expected -1)\nx=", i);
    481      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
    482      1.1  mrg       exit (1);
    483      1.1  mrg     }
    484      1.1  mrg   fclose (fp);
    485      1.1  mrg 
    486      1.1  mrg   mpfr_set_prec (x, 606);
    487      1.1  mrg   mpfr_set_prec (y, 606);
    488      1.1  mrg   mpfr_set_prec (z, 606);
    489      1.1  mrg 
    490      1.1  mrg   mpfr_set_str (y, "-f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff92daefc3f8052ca9f58736564d9e93e62d324@-1", 16, MPFR_RNDN);
    491      1.1  mrg   mpfr_set_str (z, "-f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff92daefc3f8052ca9f58736564d9e93e62d324@-1", 16, MPFR_RNDN);
    492      1.1  mrg   i = mpfr_mul (x, y, z, MPFR_RNDU);
    493      1.1  mrg   mpfr_set_str (y, "f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff25b5df87f00a5953eb0e6cac9b3d27cc5a64c@-1", 16, MPFR_RNDN);
    494      1.1  mrg   if (mpfr_cmp (x, y) || i <= 0)
    495      1.1  mrg     {
    496      1.1  mrg       printf ("Regression test (2) failed! (i=%d - Expected 1)\n", i);
    497      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN); putchar ('\n');
    498      1.1  mrg       exit (1);
    499      1.1  mrg     }
    500      1.1  mrg 
    501      1.1  mrg   mpfr_set_prec (x, 184);
    502      1.1  mrg   mpfr_set_prec (y, 92);
    503      1.1  mrg   mpfr_set_prec (z, 1023);
    504      1.1  mrg 
    505      1.1  mrg   mpfr_set_str (y, "6.9b8c8498882770d8038c3b0@-1", 16, MPFR_RNDN);
    506      1.1  mrg   mpfr_set_str (z, "7.44e24b986e7fb296f1e936ce749fec3504cbf0d5ba769466b1c9f1578115efd5d29b4c79271191a920a99280c714d3a657ad6e3afbab77ffce9d697e9bb9110e26d676069afcea8b69f1d1541f2365042d80a97c21dcccd8ace4f1bb58b49922003e738e6f37bb82ef653cb2e87f763974e6ae50ae54e7724c38b80653e3289@255", 16, MPFR_RNDN);
    507      1.1  mrg   i = mpfr_mul (x, y, z, MPFR_RNDU);
    508      1.1  mrg   mpfr_set_prec (y, 184);
    509      1.1  mrg   mpfr_set_str (y, "3.0080038f2ac5054e3e71ccbb95f76aaab2221715025a28@255",
    510      1.1  mrg                 16, MPFR_RNDN);
    511      1.1  mrg   if (mpfr_cmp (x, y) || i <= 0)
    512      1.1  mrg     {
    513      1.1  mrg       printf ("Regression test (4) failed! (i=%d - expected 1)\n", i);
    514      1.1  mrg       printf ("Ref: 3.0080038f2ac5054e3e71ccbb95f76aaab2221715025a28@255\n"
    515      1.1  mrg               "Got: ");
    516      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
    517      1.1  mrg       printf ("\n");
    518      1.1  mrg       exit (1);
    519      1.1  mrg     }
    520      1.1  mrg 
    521      1.1  mrg   mpfr_set_prec (x, 908);
    522      1.1  mrg   mpfr_set_prec (y, 908);
    523      1.1  mrg   mpfr_set_prec (z, 908);
    524      1.1  mrg   mpfr_set_str (y, "-f.fffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    525      1.1  mrg "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    526      1.1  mrg "ffffffffffffffffffffffffffffffffffffffffffffffffffffff99be91f83ec6f0ed28a3d42"
    527      1.1  mrg "e6e9a327230345ea6@-1", 16, MPFR_RNDN);
    528      1.1  mrg   mpfr_set_str (z, "-f.fffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    529      1.1  mrg "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    530      1.1  mrg "ffffffffffffffffffffffffffffffffffffffffffffffffffffff99be91f83ec6f0ed28a3d42"
    531      1.1  mrg                 "e6e9a327230345ea6@-1", 16, MPFR_RNDN);
    532      1.1  mrg   i = mpfr_mul (x, y, z, MPFR_RNDU);
    533      1.1  mrg   mpfr_set_str (y, "f.ffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    534      1.1  mrg "fffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff"
    535      1.1  mrg "fffffffffffffffffffffffffffffffffffffffffffffffffffff337d23f07d8de1da5147a85c"
    536      1.1  mrg "dd3464e46068bd4d@-1", 16, MPFR_RNDN);
    537      1.1  mrg   if (mpfr_cmp (x, y) || i <= 0)
    538      1.1  mrg     {
    539      1.1  mrg       printf ("Regression test (5) failed! (i=%d - expected 1)\n", i);
    540      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
    541      1.1  mrg       printf ("\n");
    542      1.1  mrg       exit (1);
    543      1.1  mrg     }
    544      1.1  mrg 
    545      1.1  mrg 
    546      1.1  mrg   mpfr_set_prec (x, 50);
    547      1.1  mrg   mpfr_set_prec (y, 40);
    548      1.1  mrg   mpfr_set_prec (z, 53);
    549      1.1  mrg   mpfr_set_str (y, "4.1ffffffff8", 16, MPFR_RNDN);
    550      1.1  mrg   mpfr_set_str (z, "4.2000000ffe0000@-4", 16, MPFR_RNDN);
    551      1.1  mrg   i = mpfr_mul (x, y, z, MPFR_RNDN);
    552      1.1  mrg   if (mpfr_cmp_str (x, "1.104000041d6c0@-3", 16, MPFR_RNDN) != 0
    553      1.1  mrg       || i <= 0)
    554      1.1  mrg     {
    555      1.1  mrg       printf ("Regression test (6) failed! (i=%d - expected 1)\nx=", i);
    556      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
    557      1.1  mrg       printf ("\nMore prec=");
    558      1.1  mrg       mpfr_set_prec (x, 93);
    559      1.1  mrg       mpfr_mul (x, y, z, MPFR_RNDN);
    560      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
    561      1.1  mrg       printf ("\n");
    562      1.1  mrg       exit (1);
    563      1.1  mrg     }
    564      1.1  mrg 
    565      1.1  mrg   mpfr_set_prec (x, 439);
    566      1.1  mrg   mpfr_set_prec (y, 393);
    567      1.1  mrg   mpfr_set_str (y, "-1.921fb54442d18469898cc51701b839a252049c1114cf98e804177d"
    568      1.1  mrg                 "4c76273644a29410f31c6809bbdf2a33679a748636600",
    569      1.1  mrg                 16, MPFR_RNDN);
    570      1.1  mrg   i = mpfr_mul (x, y, y, MPFR_RNDU);
    571      1.1  mrg   if (mpfr_cmp_str (x, "2.77a79937c8bbcb495b89b36602306b1c2159a8ff834288a19a08"
    572      1.1  mrg     "84094f1cda3dc426da61174c4544a173de83c2500f8bfea2e0569e3698",
    573      1.1  mrg                     16, MPFR_RNDN) != 0
    574      1.1  mrg       || i <= 0)
    575      1.1  mrg     {
    576      1.1  mrg       printf ("Regression test (7) failed! (i=%d - expected 1)\nx=", i);
    577      1.1  mrg       mpfr_out_str (stdout, 16, 0, x, MPFR_RNDN);
    578      1.1  mrg       printf ("\n");
    579      1.1  mrg       exit (1);
    580      1.1  mrg     }
    581      1.1  mrg 
    582      1.1  mrg   mpfr_set_prec (x, 1023);
    583      1.1  mrg   mpfr_set_prec (y, 1023);
    584      1.1  mrg   mpfr_set_prec (z, 511);
    585      1.1  mrg   mpfr_set_ui (x, 17, MPFR_RNDN);
    586      1.1  mrg   mpfr_set_ui (y, 42, MPFR_RNDN);
    587      1.1  mrg   i = mpfr_mul (z, x, y, MPFR_RNDN);
    588      1.1  mrg   if (mpfr_cmp_ui (z, 17*42) != 0 || i != 0)
    589      1.1  mrg     {
    590      1.1  mrg       printf ("Regression test (8) failed! (i=%d - expected 0)\nz=", i);
    591      1.1  mrg       mpfr_out_str (stdout, 16, 0, z, MPFR_RNDN);
    592      1.1  mrg       printf ("\n");
    593      1.1  mrg       exit (1);
    594      1.1  mrg     }
    595      1.1  mrg 
    596      1.1  mrg   mpfr_clears (x, y, z, (mpfr_ptr) 0);
    597      1.1  mrg }
    598      1.1  mrg 
    599      1.1  mrg #define TEST_FUNCTION test_mul
    600      1.1  mrg #define TWO_ARGS
    601      1.1  mrg #define RAND_FUNCTION(x) mpfr_random2(x, MPFR_LIMB_SIZE (x), randlimb () % 100, RANDS)
    602      1.1  mrg #include "tgeneric.c"
    603      1.1  mrg 
    604      1.1  mrg /* multiplies x by 53-bit approximation of Pi */
    605      1.1  mrg static int
    606      1.1  mrg mpfr_mulpi (mpfr_ptr y, mpfr_srcptr x, mpfr_rnd_t r)
    607      1.1  mrg {
    608      1.1  mrg   mpfr_t z;
    609      1.1  mrg   int inex;
    610      1.1  mrg 
    611      1.1  mrg   mpfr_init2 (z, 53);
    612      1.1  mrg   mpfr_set_str_binary (z, "11.001001000011111101101010100010001000010110100011");
    613      1.1  mrg   inex = mpfr_mul (y, x, z, r);
    614      1.1  mrg   mpfr_clear (z);
    615      1.1  mrg   return inex;
    616      1.1  mrg }
    617      1.1  mrg 
    618  1.1.1.2  mrg static void
    619  1.1.1.2  mrg valgrind20110503 (void)
    620  1.1.1.2  mrg {
    621  1.1.1.2  mrg   mpfr_t a, b, c;
    622  1.1.1.2  mrg 
    623  1.1.1.2  mrg   mpfr_init2 (a, 2);
    624  1.1.1.2  mrg   mpfr_init2 (b, 2005);
    625  1.1.1.2  mrg   mpfr_init2 (c, 2);
    626  1.1.1.2  mrg 
    627  1.1.1.2  mrg   mpfr_set_ui (b, 5, MPFR_RNDN);
    628  1.1.1.2  mrg   mpfr_nextabove (b);
    629  1.1.1.2  mrg   mpfr_set_ui (c, 1, MPFR_RNDN);
    630  1.1.1.2  mrg   mpfr_mul (a, b, c, MPFR_RNDZ);
    631  1.1.1.2  mrg   /* After the call to mpfr_mulhigh_n, valgrind complains:
    632  1.1.1.2  mrg      Conditional jump or move depends on uninitialised value(s) */
    633  1.1.1.2  mrg 
    634  1.1.1.2  mrg   mpfr_clears (a, b, c, (mpfr_ptr) 0);
    635  1.1.1.2  mrg }
    636  1.1.1.2  mrg 
    637  1.1.1.4  mrg static void
    638  1.1.1.4  mrg testall_rndf (mpfr_prec_t pmax)
    639  1.1.1.4  mrg {
    640  1.1.1.4  mrg   mpfr_t a, b, c, d;
    641  1.1.1.4  mrg   mpfr_prec_t pa, pb, pc;
    642  1.1.1.4  mrg 
    643  1.1.1.4  mrg   for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
    644  1.1.1.4  mrg     {
    645  1.1.1.4  mrg       mpfr_init2 (a, pa);
    646  1.1.1.4  mrg       mpfr_init2 (d, pa);
    647  1.1.1.4  mrg       for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
    648  1.1.1.4  mrg         {
    649  1.1.1.4  mrg           mpfr_init2 (b, pb);
    650  1.1.1.4  mrg           mpfr_set_ui (b, 1, MPFR_RNDN);
    651  1.1.1.4  mrg           while (mpfr_cmp_ui (b, 2) < 0)
    652  1.1.1.4  mrg             {
    653  1.1.1.4  mrg               for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
    654  1.1.1.4  mrg                 {
    655  1.1.1.4  mrg                   mpfr_init2 (c, pc);
    656  1.1.1.4  mrg                   mpfr_set_ui (c, 1, MPFR_RNDN);
    657  1.1.1.4  mrg                   while (mpfr_cmp_ui (c, 2) < 0)
    658  1.1.1.4  mrg                     {
    659  1.1.1.4  mrg                       mpfr_mul (a, b, c, MPFR_RNDF);
    660  1.1.1.4  mrg                       mpfr_mul (d, b, c, MPFR_RNDD);
    661  1.1.1.4  mrg                       if (!mpfr_equal_p (a, d))
    662  1.1.1.4  mrg                         {
    663  1.1.1.4  mrg                           mpfr_mul (d, b, c, MPFR_RNDU);
    664  1.1.1.4  mrg                           if (!mpfr_equal_p (a, d))
    665  1.1.1.4  mrg                             {
    666  1.1.1.4  mrg                               printf ("Error: mpfr_mul(a,b,c,RNDF) does not "
    667  1.1.1.4  mrg                                       "match RNDD/RNDU\n");
    668  1.1.1.4  mrg                               printf ("b="); mpfr_dump (b);
    669  1.1.1.4  mrg                               printf ("c="); mpfr_dump (c);
    670  1.1.1.4  mrg                               printf ("a="); mpfr_dump (a);
    671  1.1.1.4  mrg                               exit (1);
    672  1.1.1.4  mrg                             }
    673  1.1.1.4  mrg                         }
    674  1.1.1.4  mrg                       mpfr_nextabove (c);
    675  1.1.1.4  mrg                     }
    676  1.1.1.4  mrg                   mpfr_clear (c);
    677  1.1.1.4  mrg                 }
    678  1.1.1.4  mrg               mpfr_nextabove (b);
    679  1.1.1.4  mrg             }
    680  1.1.1.4  mrg           mpfr_clear (b);
    681  1.1.1.4  mrg         }
    682  1.1.1.4  mrg       mpfr_clear (a);
    683  1.1.1.4  mrg       mpfr_clear (d);
    684  1.1.1.4  mrg     }
    685  1.1.1.4  mrg }
    686  1.1.1.4  mrg 
    687  1.1.1.3  mrg /* Check underflow flag corresponds to *after* rounding.
    688  1.1.1.3  mrg  *
    689  1.1.1.3  mrg  * More precisely, we want to test mpfr_mul on inputs b and c such that
    690  1.1.1.3  mrg  * EXP(b*c) < emin but EXP(round(b*c, p, rnd)) = emin. Thus an underflow
    691  1.1.1.3  mrg  * must not be generated.
    692  1.1.1.3  mrg  */
    693  1.1.1.3  mrg static void
    694  1.1.1.3  mrg test_underflow (mpfr_prec_t pmax)
    695  1.1.1.3  mrg {
    696  1.1.1.3  mrg   mpfr_exp_t emin;
    697  1.1.1.3  mrg   mpfr_prec_t p;
    698  1.1.1.3  mrg   mpfr_t a0, a, b, c;
    699  1.1.1.3  mrg   int inex;
    700  1.1.1.3  mrg 
    701  1.1.1.3  mrg   mpfr_init2 (a0, MPFR_PREC_MIN);
    702  1.1.1.3  mrg   emin = mpfr_get_emin ();
    703  1.1.1.3  mrg   mpfr_setmin (a0, emin);  /* 0.5 * 2^emin */
    704  1.1.1.3  mrg 
    705  1.1.1.3  mrg   /* for RNDN, we want b*c < 0.5 * 2^emin but RNDN(b*c, p) = 0.5 * 2^emin,
    706  1.1.1.3  mrg      thus b*c >= (0.5 - 1/4 * ulp_p(0.5)) * 2^emin */
    707  1.1.1.3  mrg   for (p = MPFR_PREC_MIN; p <= pmax; p++)
    708  1.1.1.3  mrg     {
    709  1.1.1.3  mrg       mpfr_init2 (a, p + 1);
    710  1.1.1.3  mrg       mpfr_init2 (b, p + 10);
    711  1.1.1.3  mrg       mpfr_init2 (c, p + 10);
    712  1.1.1.3  mrg       do mpfr_urandomb (b, RANDS); while (MPFR_IS_ZERO (b));
    713  1.1.1.3  mrg       inex = mpfr_set_ui_2exp (a, 1, -1, MPFR_RNDZ); /* a = 0.5 */
    714  1.1.1.3  mrg       MPFR_ASSERTN (inex == 0);
    715  1.1.1.3  mrg       mpfr_nextbelow (a); /* 0.5 - 1/2*ulp_{p+1}(0.5) = 0.5 - 1/4*ulp_p(0.5) */
    716  1.1.1.3  mrg       inex = mpfr_div (c, a, b, MPFR_RNDU);
    717  1.1.1.3  mrg       /* 0.5 - 1/4 * ulp_p(0.5) = a <= b*c < 0.5 */
    718  1.1.1.3  mrg       mpfr_mul_2si (b, b, emin / 2, MPFR_RNDZ);
    719  1.1.1.3  mrg       mpfr_mul_2si (c, c, (emin - 1) / 2, MPFR_RNDZ);
    720  1.1.1.3  mrg       /* now (0.5 - 1/4 * ulp_p(0.5)) * 2^emin <= b*c < 0.5 * 2^emin,
    721  1.1.1.3  mrg          thus b*c should be rounded to 0.5 * 2^emin */
    722  1.1.1.3  mrg       mpfr_set_prec (a, p);
    723  1.1.1.3  mrg       mpfr_clear_underflow ();
    724  1.1.1.3  mrg       mpfr_mul (a, b, c, MPFR_RNDN);
    725  1.1.1.3  mrg       if (mpfr_underflow_p () || ! mpfr_equal_p (a, a0))
    726  1.1.1.3  mrg         {
    727  1.1.1.3  mrg           printf ("Error, b*c incorrect or underflow flag incorrectly set"
    728  1.1.1.3  mrg                   " for emin=%" MPFR_EXP_FSPEC "d, rnd=%s\n",
    729  1.1.1.3  mrg                   (mpfr_eexp_t) emin, mpfr_print_rnd_mode (MPFR_RNDN));
    730  1.1.1.3  mrg           printf ("b="); mpfr_dump (b);
    731  1.1.1.3  mrg           printf ("c="); mpfr_dump (c);
    732  1.1.1.3  mrg           printf ("a="); mpfr_dump (a);
    733  1.1.1.3  mrg           mpfr_set_prec (a, mpfr_get_prec (b) + mpfr_get_prec (c));
    734  1.1.1.3  mrg           mpfr_mul_2exp (b, b, 1, MPFR_RNDN);
    735  1.1.1.3  mrg           inex = mpfr_mul (a, b, c, MPFR_RNDN);
    736  1.1.1.3  mrg           MPFR_ASSERTN (inex == 0);
    737  1.1.1.3  mrg           printf ("Exact 2*a="); mpfr_dump (a);
    738  1.1.1.3  mrg           exit (1);
    739  1.1.1.3  mrg         }
    740  1.1.1.3  mrg       mpfr_clear (a);
    741  1.1.1.3  mrg       mpfr_clear (b);
    742  1.1.1.3  mrg       mpfr_clear (c);
    743  1.1.1.3  mrg     }
    744  1.1.1.3  mrg 
    745  1.1.1.3  mrg   /* for RNDU, we want b*c < 0.5*2^emin but RNDU(b*c, p) = 0.5*2^emin thus
    746  1.1.1.3  mrg      b*c > (0.5 - 1/2 * ulp_p(0.5)) * 2^emin */
    747  1.1.1.3  mrg   for (p = MPFR_PREC_MIN; p <= pmax; p++)
    748  1.1.1.3  mrg     {
    749  1.1.1.3  mrg       mpfr_init2 (a, p);
    750  1.1.1.3  mrg       mpfr_init2 (b, p + 10);
    751  1.1.1.3  mrg       mpfr_init2 (c, p + 10);
    752  1.1.1.3  mrg       do mpfr_urandomb (b, RANDS); while (MPFR_IS_ZERO (b));
    753  1.1.1.3  mrg       inex = mpfr_set_ui_2exp (a, 1, -1, MPFR_RNDZ); /* a = 0.5 */
    754  1.1.1.3  mrg       MPFR_ASSERTN (inex == 0);
    755  1.1.1.3  mrg       mpfr_nextbelow (a); /* 0.5 - 1/2 * ulp_p(0.5) */
    756  1.1.1.3  mrg       inex = mpfr_div (c, a, b, MPFR_RNDU);
    757  1.1.1.3  mrg       /* 0.5 - 1/2 * ulp_p(0.5) <= b*c < 0.5 */
    758  1.1.1.3  mrg       mpfr_mul_2si (b, b, emin / 2, MPFR_RNDZ);
    759  1.1.1.3  mrg       mpfr_mul_2si (c, c, (emin - 1) / 2, MPFR_RNDZ);
    760  1.1.1.3  mrg       if (inex == 0)
    761  1.1.1.3  mrg         mpfr_nextabove (c); /* ensures b*c > (0.5 - 1/2 * ulp_p(0.5)) * 2^emin.
    762  1.1.1.3  mrg                                Warning: for p=1, 0.5 - 1/2 * ulp_p(0.5)
    763  1.1.1.3  mrg                                = 0.25, thus b*c > 2^(emin-2), which should
    764  1.1.1.3  mrg                                also be rounded up with p=1 to 0.5 * 2^emin
    765  1.1.1.3  mrg                                with an unbounded exponent range. */
    766  1.1.1.3  mrg       mpfr_clear_underflow ();
    767  1.1.1.3  mrg       mpfr_mul (a, b, c, MPFR_RNDU);
    768  1.1.1.3  mrg       if (mpfr_underflow_p () || ! mpfr_equal_p (a, a0))
    769  1.1.1.3  mrg         {
    770  1.1.1.3  mrg           printf ("Error, b*c incorrect or underflow flag incorrectly set"
    771  1.1.1.3  mrg                   " for emin=%" MPFR_EXP_FSPEC "d, rnd=%s\n",
    772  1.1.1.3  mrg                   (mpfr_eexp_t) emin, mpfr_print_rnd_mode (MPFR_RNDU));
    773  1.1.1.3  mrg           printf ("b="); mpfr_dump (b);
    774  1.1.1.3  mrg           printf ("c="); mpfr_dump (c);
    775  1.1.1.3  mrg           printf ("a="); mpfr_dump (a);
    776  1.1.1.3  mrg           mpfr_set_prec (a, mpfr_get_prec (b) + mpfr_get_prec (c));
    777  1.1.1.3  mrg           mpfr_mul_2exp (b, b, 1, MPFR_RNDN);
    778  1.1.1.3  mrg           inex = mpfr_mul (a, b, c, MPFR_RNDN);
    779  1.1.1.3  mrg           MPFR_ASSERTN (inex == 0);
    780  1.1.1.3  mrg           printf ("Exact 2*a="); mpfr_dump (a);
    781  1.1.1.3  mrg           exit (1);
    782  1.1.1.3  mrg         }
    783  1.1.1.3  mrg       mpfr_clear (a);
    784  1.1.1.3  mrg       mpfr_clear (b);
    785  1.1.1.3  mrg       mpfr_clear (c);
    786  1.1.1.3  mrg     }
    787  1.1.1.3  mrg 
    788  1.1.1.3  mrg   mpfr_clear (a0);
    789  1.1.1.3  mrg }
    790  1.1.1.3  mrg 
    791  1.1.1.4  mrg /* checks special case where no underflow should occur */
    792  1.1.1.4  mrg static void
    793  1.1.1.4  mrg bug20161209 (void)
    794  1.1.1.4  mrg {
    795  1.1.1.4  mrg   mpfr_exp_t emin;
    796  1.1.1.4  mrg   mpfr_t x, y, z;
    797  1.1.1.4  mrg 
    798  1.1.1.4  mrg   emin = mpfr_get_emin ();
    799  1.1.1.4  mrg   set_emin (-1);
    800  1.1.1.4  mrg 
    801  1.1.1.4  mrg   /* test for mpfr_mul_1 for 64-bit limb, mpfr_mul_2 for 32-bit limb */
    802  1.1.1.4  mrg   mpfr_init2 (x, 53);
    803  1.1.1.4  mrg   mpfr_init2 (y, 53);
    804  1.1.1.4  mrg   mpfr_init2 (z, 53);
    805  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.101000001E-1"); /* x = 321/2^10 */
    806  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
    807  1.1.1.4  mrg   /* y = 28059810762433/2^46 */
    808  1.1.1.4  mrg   /* x * y = (2^53+1)/2^56 = 0.001...000[1]000..., and should round to 0.25 */
    809  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDN);
    810  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    811  1.1.1.4  mrg 
    812  1.1.1.4  mrg   /* test for mpfr_mul_2 for 64-bit limb */
    813  1.1.1.4  mrg   mpfr_set_prec (x, 65);
    814  1.1.1.4  mrg   mpfr_set_prec (y, 65);
    815  1.1.1.4  mrg   mpfr_set_prec (z, 65);
    816  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.101101000010010110100001E-1"); /* 11806113/2^25 */
    817  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
    818  1.1.1.4  mrg   /* y = 3124947910241/2^43 */
    819  1.1.1.4  mrg   /* x * y = (2^65+1)/2^68 = 0.001...000[1]000..., and should round to 0.25 */
    820  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDN);
    821  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    822  1.1.1.4  mrg 
    823  1.1.1.4  mrg   /* test for the generic code */
    824  1.1.1.4  mrg   mpfr_set_prec (x, 54);
    825  1.1.1.4  mrg   mpfr_set_prec (y, 55);
    826  1.1.1.4  mrg   mpfr_set_prec (z, 53);
    827  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.101000001E-1");
    828  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
    829  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDN);
    830  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    831  1.1.1.4  mrg 
    832  1.1.1.4  mrg   /* another test for the generic code */
    833  1.1.1.4  mrg   mpfr_set_prec (x, 66);
    834  1.1.1.4  mrg   mpfr_set_prec (y, 67);
    835  1.1.1.4  mrg   mpfr_set_prec (z, 65);
    836  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.101101000010010110100001E-1");
    837  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
    838  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDN);
    839  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    840  1.1.1.4  mrg 
    841  1.1.1.4  mrg   mpfr_clear (x);
    842  1.1.1.4  mrg   mpfr_clear (y);
    843  1.1.1.4  mrg   mpfr_clear (z);
    844  1.1.1.4  mrg   set_emin (emin);
    845  1.1.1.4  mrg }
    846  1.1.1.4  mrg 
    847  1.1.1.4  mrg /* test for case in mpfr_mul_1() where:
    848  1.1.1.4  mrg    ax = __gmpfr_emin - 1
    849  1.1.1.4  mrg    ap[0] == ~mask
    850  1.1.1.4  mrg    rnd_mode = MPFR_RNDZ.
    851  1.1.1.4  mrg    Whatever the values of rb and sb, we should round to zero (underflow). */
    852  1.1.1.4  mrg static void
    853  1.1.1.4  mrg bug20161209a (void)
    854  1.1.1.4  mrg {
    855  1.1.1.4  mrg   mpfr_exp_t emin;
    856  1.1.1.4  mrg   mpfr_t x, y, z;
    857  1.1.1.4  mrg 
    858  1.1.1.4  mrg   emin = mpfr_get_emin ();
    859  1.1.1.4  mrg   set_emin (-1);
    860  1.1.1.4  mrg 
    861  1.1.1.4  mrg   mpfr_init2 (x, 53);
    862  1.1.1.4  mrg   mpfr_init2 (y, 53);
    863  1.1.1.4  mrg   mpfr_init2 (z, 53);
    864  1.1.1.4  mrg 
    865  1.1.1.4  mrg   /* case rb = sb = 0 */
    866  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.11010010100110000110110011111E-1");
    867  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.1001101110011000110100001");
    868  1.1.1.4  mrg   /* x = 441650591/2^30, y = 20394401/2^25 */
    869  1.1.1.4  mrg   /* x * y = (2^53-1)/2^55 = 0.00111...111[0]000..., and should round to 0 */
    870  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDZ);
    871  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_zero_p (z));
    872  1.1.1.4  mrg 
    873  1.1.1.4  mrg   /* case rb = 1, sb = 0 */
    874  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.111111111000000000000000000111111111E-1");
    875  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.1000000001000000001");
    876  1.1.1.4  mrg   /* x = 68585259519/2^37, y = 262657/2^19 */
    877  1.1.1.4  mrg   /* x * y = (2^54-1)/2^56 = 0.00111...111[1]000..., and should round to 0 */
    878  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDZ);
    879  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_zero_p (z));
    880  1.1.1.4  mrg 
    881  1.1.1.4  mrg   /* case rb = 0, sb = 1 */
    882  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.110010011001011110001100100001000001E-1");
    883  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.10100010100010111101");
    884  1.1.1.4  mrg   /* x = 541144371852^37, y = 665789/2^20 */
    885  1.1.1.4  mrg   /* x * y = (2^55-3)/2^57 = 0.00111...111[0]100..., and should round to 0 */
    886  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDZ);
    887  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_zero_p (z));
    888  1.1.1.4  mrg 
    889  1.1.1.4  mrg   /* case rb = sb = 1 */
    890  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.10100110001001001010001111110010100111E-1");
    891  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.110001010011101001");
    892  1.1.1.4  mrg   /* x = 178394823847/2^39, y = 201961/2^18 */
    893  1.1.1.4  mrg   /* x * y = (2^55-1)/2^57 = 0.00111...111[1]100..., and should round to 0 */
    894  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDZ);
    895  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_zero_p (z));
    896  1.1.1.4  mrg 
    897  1.1.1.4  mrg   /* similar test for mpfr_mul_2 (we only check rb = sb = 1 here) */
    898  1.1.1.4  mrg   mpfr_set_prec (x, 65);
    899  1.1.1.4  mrg   mpfr_set_prec (y, 65);
    900  1.1.1.4  mrg   mpfr_set_prec (z, 65);
    901  1.1.1.4  mrg   /* 2^67-1 = 193707721 * 761838257287 */
    902  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.1011100010111011111011001001E-1");
    903  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.1011000101100001000110010100010010000111");
    904  1.1.1.4  mrg   mpfr_mul (z, x, y, MPFR_RNDZ);
    905  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_zero_p (z));
    906  1.1.1.4  mrg 
    907  1.1.1.4  mrg   mpfr_clear (x);
    908  1.1.1.4  mrg   mpfr_clear (y);
    909  1.1.1.4  mrg   mpfr_clear (z);
    910  1.1.1.4  mrg   set_emin (emin);
    911  1.1.1.4  mrg }
    912  1.1.1.4  mrg 
    913  1.1.1.4  mrg /* bug for RNDF */
    914  1.1.1.4  mrg static void
    915  1.1.1.4  mrg bug20170602 (void)
    916  1.1.1.4  mrg {
    917  1.1.1.4  mrg   mpfr_t x, u, y, yd, yu;
    918  1.1.1.4  mrg 
    919  1.1.1.4  mrg   mpfr_init2 (x, 493);
    920  1.1.1.4  mrg   mpfr_init2 (u, 493);
    921  1.1.1.4  mrg   mpfr_init2 (y, 503);
    922  1.1.1.4  mrg   mpfr_init2 (yd, 503);
    923  1.1.1.4  mrg   mpfr_init2 (yu, 503);
    924  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.1111100000000000001111111110000000001111111111111000000000000000000011111111111111111111111000000000000000000001111111111111111111111111111111111111111000000000011111111111111111111000000000011111111111111000000000000001110000000000000000000000000000000000000000011111111111110011111111111100000000000000011111111111111111110000000011111111111111111110011111111111110000000000001111111111111111000000000000000000000000000000000000111111111111111111111111111111011111111111111111111111111111100E44");
    925  1.1.1.4  mrg   mpfr_set_str_binary (u, "0.1110000000000000001111111111111111111111111111111111111000000000111111111111111111111111111111000000000000000000001111111000000000000000011111111111111111111111111111111111111111111111111111111000000000000000011111111111111000000011111111111111111110000000000000001111111111111111111111111111111111111110000000000001111111111111111111111111111111111111000000000000000000000000000000000001111111111111000000000000000000001111111111100000000000000011111111111111111111111111111111111111111111111E35");
    926  1.1.1.4  mrg   mpfr_mul (y, x, u, MPFR_RNDF);
    927  1.1.1.4  mrg   mpfr_mul (yd, x, u, MPFR_RNDD);
    928  1.1.1.4  mrg   mpfr_mul (yu, x, u, MPFR_RNDU);
    929  1.1.1.4  mrg   if (mpfr_cmp (y, yd) != 0 && mpfr_cmp (y, yu) != 0)
    930  1.1.1.4  mrg     {
    931  1.1.1.4  mrg       printf ("RNDF is neither RNDD nor RNDU\n");
    932  1.1.1.4  mrg       printf ("x"); mpfr_dump (x);
    933  1.1.1.4  mrg       printf ("u"); mpfr_dump (u);
    934  1.1.1.4  mrg       printf ("y(RNDF)="); mpfr_dump (y);
    935  1.1.1.4  mrg       printf ("y(RNDD)="); mpfr_dump (yd);
    936  1.1.1.4  mrg       printf ("y(RNDU)="); mpfr_dump (yu);
    937  1.1.1.4  mrg       exit (1);
    938  1.1.1.4  mrg     }
    939  1.1.1.4  mrg   mpfr_clear (x);
    940  1.1.1.4  mrg   mpfr_clear (u);
    941  1.1.1.4  mrg   mpfr_clear (y);
    942  1.1.1.4  mrg   mpfr_clear (yd);
    943  1.1.1.4  mrg   mpfr_clear (yu);
    944  1.1.1.4  mrg }
    945  1.1.1.4  mrg 
    946  1.1.1.4  mrg /* Test for 1 to 3 limbs. */
    947  1.1.1.4  mrg static void
    948  1.1.1.4  mrg small_prec (void)
    949  1.1.1.4  mrg {
    950  1.1.1.4  mrg   mpfr_exp_t emin, emax;
    951  1.1.1.4  mrg   mpfr_t x, y, z1, z2, zz;
    952  1.1.1.4  mrg   int xq, yq, zq;
    953  1.1.1.4  mrg   mpfr_rnd_t rnd;
    954  1.1.1.4  mrg   mpfr_flags_t flags1, flags2;
    955  1.1.1.4  mrg   int inex1, inex2;
    956  1.1.1.4  mrg   int i, j, r, s, ediff;
    957  1.1.1.4  mrg 
    958  1.1.1.4  mrg   emin = mpfr_get_emin ();
    959  1.1.1.4  mrg   emax = mpfr_get_emax ();
    960  1.1.1.4  mrg 
    961  1.1.1.4  mrg   /* The mpfr_mul implementation doesn't extend the exponent range,
    962  1.1.1.4  mrg      so that it is OK to reduce it here for the test to make sure that
    963  1.1.1.4  mrg      mpfr_mul_2si can be used. */
    964  1.1.1.4  mrg   set_emin (-1000);
    965  1.1.1.4  mrg   set_emax (1000);
    966  1.1.1.4  mrg 
    967  1.1.1.4  mrg   mpfr_inits2 (3 * GMP_NUMB_BITS, x, y, z1, z2, (mpfr_ptr) 0);
    968  1.1.1.4  mrg   mpfr_init2 (zz, 6 * GMP_NUMB_BITS);
    969  1.1.1.4  mrg   for (i = 0; i < 3; i++)
    970  1.1.1.4  mrg     for (j = 0; j < 10000; j++)
    971  1.1.1.4  mrg       {
    972  1.1.1.4  mrg         xq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
    973  1.1.1.4  mrg         mpfr_set_prec (x, xq);
    974  1.1.1.4  mrg         yq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
    975  1.1.1.4  mrg         mpfr_set_prec (y, yq);
    976  1.1.1.4  mrg         zq = i * GMP_NUMB_BITS + 1 + randlimb () % (GMP_NUMB_BITS-1);
    977  1.1.1.4  mrg         mpfr_set_prec (z1, zq);
    978  1.1.1.4  mrg         mpfr_set_prec (z2, zq);
    979  1.1.1.4  mrg         s = r = randlimb () & 0x7f;
    980  1.1.1.4  mrg         do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
    981  1.1.1.4  mrg         if (s & 1)
    982  1.1.1.4  mrg           mpfr_neg (x, x, MPFR_RNDN);
    983  1.1.1.4  mrg         s >>= 1;
    984  1.1.1.4  mrg         if (s & 1)
    985  1.1.1.4  mrg           {
    986  1.1.1.4  mrg             do mpfr_urandomb (y, RANDS); while (mpfr_zero_p (y));
    987  1.1.1.4  mrg           }
    988  1.1.1.4  mrg         else
    989  1.1.1.4  mrg           {
    990  1.1.1.4  mrg             mpfr_ui_div (y, 1, x, MPFR_RNDN);
    991  1.1.1.4  mrg             mpfr_set_exp (y, 0);
    992  1.1.1.4  mrg           }
    993  1.1.1.4  mrg         s >>= 1;
    994  1.1.1.4  mrg         if (s & 1)
    995  1.1.1.4  mrg           mpfr_neg (y, y, MPFR_RNDN);
    996  1.1.1.4  mrg         s >>= 1;
    997  1.1.1.4  mrg         rnd = RND_RAND_NO_RNDF ();
    998  1.1.1.4  mrg         inex1 = mpfr_mul (zz, x, y, MPFR_RNDN);
    999  1.1.1.4  mrg         MPFR_ASSERTN (inex1 == 0);
   1000  1.1.1.4  mrg         if (s == 0)
   1001  1.1.1.4  mrg           {
   1002  1.1.1.4  mrg             ediff = __gmpfr_emin - MPFR_EXP (x);
   1003  1.1.1.4  mrg             mpfr_set_exp (x, __gmpfr_emin);
   1004  1.1.1.4  mrg           }
   1005  1.1.1.4  mrg         else if (s == 1)
   1006  1.1.1.4  mrg           {
   1007  1.1.1.4  mrg             ediff = __gmpfr_emax - MPFR_EXP (x) + 1;
   1008  1.1.1.4  mrg             mpfr_set_exp (x, __gmpfr_emax);
   1009  1.1.1.4  mrg             mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
   1010  1.1.1.4  mrg           }
   1011  1.1.1.4  mrg         else
   1012  1.1.1.4  mrg           ediff = 0;
   1013  1.1.1.4  mrg         mpfr_clear_flags ();
   1014  1.1.1.4  mrg         inex1 = mpfr_mul_2si (z1, zz, ediff, rnd);
   1015  1.1.1.4  mrg         flags1 = __gmpfr_flags;
   1016  1.1.1.4  mrg         mpfr_clear_flags ();
   1017  1.1.1.4  mrg         inex2 = mpfr_mul (z2, x, y, rnd);
   1018  1.1.1.4  mrg         flags2 = __gmpfr_flags;
   1019  1.1.1.4  mrg         if (!(mpfr_equal_p (z1, z2) &&
   1020  1.1.1.4  mrg               SAME_SIGN (inex1, inex2) &&
   1021  1.1.1.4  mrg               flags1 == flags2))
   1022  1.1.1.4  mrg           {
   1023  1.1.1.4  mrg             printf ("Error in small_prec on i = %d, j = %d\n", i, j);
   1024  1.1.1.4  mrg             printf ("r = 0x%x, xq = %d, yq = %d, zq = %d, rnd = %s\n",
   1025  1.1.1.4  mrg                     r, xq, yq, zq, mpfr_print_rnd_mode (rnd));
   1026  1.1.1.4  mrg             printf ("x = ");
   1027  1.1.1.4  mrg             mpfr_dump (x);
   1028  1.1.1.4  mrg             printf ("y = ");
   1029  1.1.1.4  mrg             mpfr_dump (y);
   1030  1.1.1.4  mrg             printf ("ediff = %d\n", ediff);
   1031  1.1.1.4  mrg             printf ("zz = ");
   1032  1.1.1.4  mrg             mpfr_dump (zz);
   1033  1.1.1.4  mrg             printf ("Expected ");
   1034  1.1.1.4  mrg             mpfr_dump (z1);
   1035  1.1.1.4  mrg             printf ("with inex = %d and flags =", inex1);
   1036  1.1.1.4  mrg             flags_out (flags1);
   1037  1.1.1.4  mrg             printf ("Got      ");
   1038  1.1.1.4  mrg             mpfr_dump (z2);
   1039  1.1.1.4  mrg             printf ("with inex = %d and flags =", inex2);
   1040  1.1.1.4  mrg             flags_out (flags2);
   1041  1.1.1.4  mrg             exit (1);
   1042  1.1.1.4  mrg           }
   1043  1.1.1.4  mrg       }
   1044  1.1.1.4  mrg   mpfr_clears (x, y, z1, z2, zz, (mpfr_ptr) 0);
   1045  1.1.1.4  mrg 
   1046  1.1.1.4  mrg   set_emin (emin);
   1047  1.1.1.4  mrg   set_emax (emax);
   1048  1.1.1.4  mrg }
   1049  1.1.1.4  mrg 
   1050  1.1.1.4  mrg /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */
   1051  1.1.1.4  mrg static void
   1052  1.1.1.4  mrg test_underflow2 (void)
   1053  1.1.1.4  mrg {
   1054  1.1.1.4  mrg   mpfr_t x, y, z;
   1055  1.1.1.4  mrg   mpfr_exp_t emin;
   1056  1.1.1.4  mrg 
   1057  1.1.1.4  mrg   emin = mpfr_get_emin ();
   1058  1.1.1.4  mrg   mpfr_set_emin (0);
   1059  1.1.1.4  mrg 
   1060  1.1.1.4  mrg   mpfr_init2 (x, 24);
   1061  1.1.1.4  mrg   mpfr_init2 (y, 24);
   1062  1.1.1.4  mrg   mpfr_init2 (z, 24);
   1063  1.1.1.4  mrg 
   1064  1.1.1.4  mrg   mpfr_set_ui_2exp (x, 12913, -14, MPFR_RNDN);
   1065  1.1.1.4  mrg   mpfr_set_ui_2exp (y, 5197,  -13, MPFR_RNDN);
   1066  1.1.1.4  mrg   mpfr_clear_underflow ();
   1067  1.1.1.4  mrg   /* x*y = 0.0111111111111111111111111[01] thus underflow */
   1068  1.1.1.4  mrg   mpfr_mul (z, y, x, MPFR_RNDN);
   1069  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
   1070  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_underflow_p ());
   1071  1.1.1.4  mrg 
   1072  1.1.1.4  mrg   mpfr_set_prec (x, 65);
   1073  1.1.1.4  mrg   mpfr_set_prec (y, 65);
   1074  1.1.1.4  mrg   mpfr_set_prec (z, 65);
   1075  1.1.1.4  mrg   mpfr_set_str_binary (x, "0.10011101000110111011111100101001111");
   1076  1.1.1.4  mrg   mpfr_set_str_binary (y, "0.110100001001000111000011011110011");
   1077  1.1.1.4  mrg   mpfr_clear_underflow ();
   1078  1.1.1.4  mrg   /* x*y = 0.011111111111111111111111111111111111111111111111111111111111111111[01] thus underflow */
   1079  1.1.1.4  mrg   mpfr_mul (z, y, x, MPFR_RNDN);
   1080  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
   1081  1.1.1.4  mrg   MPFR_ASSERTN(mpfr_underflow_p ());
   1082  1.1.1.4  mrg 
   1083  1.1.1.4  mrg   mpfr_clear (y);
   1084  1.1.1.4  mrg   mpfr_clear (x);
   1085  1.1.1.4  mrg   mpfr_clear (z);
   1086  1.1.1.4  mrg 
   1087  1.1.1.4  mrg   mpfr_set_emin (emin);
   1088  1.1.1.4  mrg }
   1089  1.1.1.4  mrg 
   1090      1.1  mrg int
   1091      1.1  mrg main (int argc, char *argv[])
   1092      1.1  mrg {
   1093      1.1  mrg   tests_start_mpfr ();
   1094      1.1  mrg 
   1095  1.1.1.4  mrg   testall_rndf (9);
   1096      1.1  mrg   check_nans ();
   1097      1.1  mrg   check_exact ();
   1098      1.1  mrg   check_float ();
   1099      1.1  mrg 
   1100      1.1  mrg   check53("6.9314718055994530941514e-1", "0.0", MPFR_RNDZ, "0.0");
   1101      1.1  mrg   check53("0.0", "6.9314718055994530941514e-1", MPFR_RNDZ, "0.0");
   1102      1.1  mrg   check_sign();
   1103      1.1  mrg   check53("-4.165000000e4", "-0.00004801920768307322868063274915", MPFR_RNDN,
   1104      1.1  mrg           "2.0");
   1105      1.1  mrg   check53("2.71331408349172961467e-08", "-6.72658901114033715233e-165",
   1106      1.1  mrg           MPFR_RNDZ, "-1.8251348697787782844e-172");
   1107      1.1  mrg   check53("2.71331408349172961467e-08", "-6.72658901114033715233e-165",
   1108      1.1  mrg           MPFR_RNDA, "-1.8251348697787786e-172");
   1109      1.1  mrg   check53("0.31869277231188065", "0.88642843322303122", MPFR_RNDZ,
   1110      1.1  mrg           "2.8249833483992453642e-1");
   1111      1.1  mrg   check("8.47622108205396074254e-01", "3.24039313247872939883e-01", MPFR_RNDU,
   1112      1.1  mrg         28, 45, 2, "0.375");
   1113      1.1  mrg   check("8.47622108205396074254e-01", "3.24039313247872939883e-01", MPFR_RNDA,
   1114      1.1  mrg         28, 45, 2, "0.375");
   1115      1.1  mrg   check("2.63978122803639081440e-01", "6.8378615379333496093e-1", MPFR_RNDN,
   1116      1.1  mrg         34, 23, 31, "0.180504585267044603");
   1117      1.1  mrg   check("1.0", "0.11835170935876249132", MPFR_RNDU, 6, 41, 36,
   1118      1.1  mrg         "0.1183517093595583");
   1119      1.1  mrg   check53("67108865.0", "134217729.0", MPFR_RNDN, "9.007199456067584e15");
   1120      1.1  mrg   check("1.37399642157394197284e-01", "2.28877275604219221350e-01", MPFR_RNDN,
   1121      1.1  mrg         49, 15, 32, "0.0314472340833162888");
   1122      1.1  mrg   check("4.03160720978664954828e-01", "5.854828e-1"
   1123      1.1  mrg         /*"5.85483042917246621073e-01"*/, MPFR_RNDZ,
   1124      1.1  mrg         51, 22, 32, "0.2360436821472831");
   1125      1.1  mrg   check("3.90798504668055102229e-14", "9.85394674650308388664e-04", MPFR_RNDN,
   1126      1.1  mrg         46, 22, 12, "0.385027296503914762e-16");
   1127      1.1  mrg   check("4.58687081072827851358e-01", "2.20543551472118792844e-01", MPFR_RNDN,
   1128      1.1  mrg         49, 3, 2, "0.09375");
   1129      1.1  mrg   check_max();
   1130      1.1  mrg   check_min();
   1131  1.1.1.4  mrg   small_prec ();
   1132      1.1  mrg 
   1133      1.1  mrg   check_regression ();
   1134  1.1.1.4  mrg   test_generic (MPFR_PREC_MIN, 500, 100);
   1135      1.1  mrg 
   1136      1.1  mrg   data_check ("data/mulpi", mpfr_mulpi, "mpfr_mulpi");
   1137      1.1  mrg 
   1138  1.1.1.2  mrg   valgrind20110503 ();
   1139  1.1.1.3  mrg   test_underflow (128);
   1140  1.1.1.4  mrg   bug20161209 ();
   1141  1.1.1.4  mrg   bug20161209a ();
   1142  1.1.1.4  mrg   bug20170602 ();
   1143  1.1.1.4  mrg   test_underflow2 ();
   1144  1.1.1.2  mrg 
   1145      1.1  mrg   tests_end_mpfr ();
   1146      1.1  mrg   return 0;
   1147      1.1  mrg }
   1148