Home | History | Annotate | Line # | Download | only in tests
tmul.c revision 1.1.1.3.4.1
      1          1.1       mrg /* Test file for mpfr_mul.
      2          1.1       mrg 
      3  1.1.1.3.4.1  christos 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.3.4.1  christos       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.3.4.1  christos       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
     69  1.1.1.3.4.1  christos       putchar ('\n');
     70          1.1       mrg       exit (1);
     71          1.1       mrg     }
     72  1.1.1.3.4.1  christos   mpfr_clear (xx);
     73  1.1.1.3.4.1  christos   mpfr_clear (yy);
     74  1.1.1.3.4.1  christos   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.3.4.1  christos       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
     92  1.1.1.3.4.1  christos       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.3.4.1  christos       mpfr_out_str (stdout, 10, 0, zz, MPFR_RNDN);
    114  1.1.1.3.4.1  christos       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.3.4.1  christos       mpfr_out_str (stdout, 10, 0, a, MPFR_RNDN);
    181  1.1.1.3.4.1  christos       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.3.4.1  christos   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.3.4.1  christos       /* for prec=1, ensure PREC(c) >= 1 */
    217  1.1.1.3.4.1  christos       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.3.4.1  christos           if ((inexact == 0) && mpfr_cmp (c, d) && rnd != MPFR_RNDF)
    231          1.1       mrg             {
    232  1.1.1.3.4.1  christos               printf ("rnd=%s: inexact=0 but results differ\n",
    233  1.1.1.3.4.1  christos                       mpfr_print_rnd_mode (rnd));
    234  1.1.1.3.4.1  christos               printf ("a=");
    235  1.1.1.3.4.1  christos               mpfr_out_str (stdout, 2, 0, a, rnd);
    236  1.1.1.3.4.1  christos               printf ("\nb=");
    237  1.1.1.3.4.1  christos               mpfr_out_str (stdout, 2, 0, b, rnd);
    238  1.1.1.3.4.1  christos               printf ("\nc=");
    239  1.1.1.3.4.1  christos               mpfr_out_str (stdout, 2, 0, c, rnd);
    240  1.1.1.3.4.1  christos               printf ("\nd=");
    241  1.1.1.3.4.1  christos               mpfr_out_str (stdout, 2, 0, d, rnd);
    242  1.1.1.3.4.1  christos               printf ("\n");
    243          1.1       mrg               exit (1);
    244          1.1       mrg             }
    245  1.1.1.3.4.1  christos           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.3.4.1  christos   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.3.4.1  christos static void
    638  1.1.1.3.4.1  christos testall_rndf (mpfr_prec_t pmax)
    639  1.1.1.3.4.1  christos {
    640  1.1.1.3.4.1  christos   mpfr_t a, b, c, d;
    641  1.1.1.3.4.1  christos   mpfr_prec_t pa, pb, pc;
    642  1.1.1.3.4.1  christos 
    643  1.1.1.3.4.1  christos   for (pa = MPFR_PREC_MIN; pa <= pmax; pa++)
    644  1.1.1.3.4.1  christos     {
    645  1.1.1.3.4.1  christos       mpfr_init2 (a, pa);
    646  1.1.1.3.4.1  christos       mpfr_init2 (d, pa);
    647  1.1.1.3.4.1  christos       for (pb = MPFR_PREC_MIN; pb <= pmax; pb++)
    648  1.1.1.3.4.1  christos         {
    649  1.1.1.3.4.1  christos           mpfr_init2 (b, pb);
    650  1.1.1.3.4.1  christos           mpfr_set_ui (b, 1, MPFR_RNDN);
    651  1.1.1.3.4.1  christos           while (mpfr_cmp_ui (b, 2) < 0)
    652  1.1.1.3.4.1  christos             {
    653  1.1.1.3.4.1  christos               for (pc = MPFR_PREC_MIN; pc <= pmax; pc++)
    654  1.1.1.3.4.1  christos                 {
    655  1.1.1.3.4.1  christos                   mpfr_init2 (c, pc);
    656  1.1.1.3.4.1  christos                   mpfr_set_ui (c, 1, MPFR_RNDN);
    657  1.1.1.3.4.1  christos                   while (mpfr_cmp_ui (c, 2) < 0)
    658  1.1.1.3.4.1  christos                     {
    659  1.1.1.3.4.1  christos                       mpfr_mul (a, b, c, MPFR_RNDF);
    660  1.1.1.3.4.1  christos                       mpfr_mul (d, b, c, MPFR_RNDD);
    661  1.1.1.3.4.1  christos                       if (!mpfr_equal_p (a, d))
    662  1.1.1.3.4.1  christos                         {
    663  1.1.1.3.4.1  christos                           mpfr_mul (d, b, c, MPFR_RNDU);
    664  1.1.1.3.4.1  christos                           if (!mpfr_equal_p (a, d))
    665  1.1.1.3.4.1  christos                             {
    666  1.1.1.3.4.1  christos                               printf ("Error: mpfr_mul(a,b,c,RNDF) does not "
    667  1.1.1.3.4.1  christos                                       "match RNDD/RNDU\n");
    668  1.1.1.3.4.1  christos                               printf ("b="); mpfr_dump (b);
    669  1.1.1.3.4.1  christos                               printf ("c="); mpfr_dump (c);
    670  1.1.1.3.4.1  christos                               printf ("a="); mpfr_dump (a);
    671  1.1.1.3.4.1  christos                               exit (1);
    672  1.1.1.3.4.1  christos                             }
    673  1.1.1.3.4.1  christos                         }
    674  1.1.1.3.4.1  christos                       mpfr_nextabove (c);
    675  1.1.1.3.4.1  christos                     }
    676  1.1.1.3.4.1  christos                   mpfr_clear (c);
    677  1.1.1.3.4.1  christos                 }
    678  1.1.1.3.4.1  christos               mpfr_nextabove (b);
    679  1.1.1.3.4.1  christos             }
    680  1.1.1.3.4.1  christos           mpfr_clear (b);
    681  1.1.1.3.4.1  christos         }
    682  1.1.1.3.4.1  christos       mpfr_clear (a);
    683  1.1.1.3.4.1  christos       mpfr_clear (d);
    684  1.1.1.3.4.1  christos     }
    685  1.1.1.3.4.1  christos }
    686  1.1.1.3.4.1  christos 
    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.3.4.1  christos /* checks special case where no underflow should occur */
    792  1.1.1.3.4.1  christos static void
    793  1.1.1.3.4.1  christos bug20161209 (void)
    794  1.1.1.3.4.1  christos {
    795  1.1.1.3.4.1  christos   mpfr_exp_t emin;
    796  1.1.1.3.4.1  christos   mpfr_t x, y, z;
    797  1.1.1.3.4.1  christos 
    798  1.1.1.3.4.1  christos   emin = mpfr_get_emin ();
    799  1.1.1.3.4.1  christos   set_emin (-1);
    800  1.1.1.3.4.1  christos 
    801  1.1.1.3.4.1  christos   /* test for mpfr_mul_1 for 64-bit limb, mpfr_mul_2 for 32-bit limb */
    802  1.1.1.3.4.1  christos   mpfr_init2 (x, 53);
    803  1.1.1.3.4.1  christos   mpfr_init2 (y, 53);
    804  1.1.1.3.4.1  christos   mpfr_init2 (z, 53);
    805  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.101000001E-1"); /* x = 321/2^10 */
    806  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
    807  1.1.1.3.4.1  christos   /* y = 28059810762433/2^46 */
    808  1.1.1.3.4.1  christos   /* x * y = (2^53+1)/2^56 = 0.001...000[1]000..., and should round to 0.25 */
    809  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDN);
    810  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    811  1.1.1.3.4.1  christos 
    812  1.1.1.3.4.1  christos   /* test for mpfr_mul_2 for 64-bit limb */
    813  1.1.1.3.4.1  christos   mpfr_set_prec (x, 65);
    814  1.1.1.3.4.1  christos   mpfr_set_prec (y, 65);
    815  1.1.1.3.4.1  christos   mpfr_set_prec (z, 65);
    816  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.101101000010010110100001E-1"); /* 11806113/2^25 */
    817  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
    818  1.1.1.3.4.1  christos   /* y = 3124947910241/2^43 */
    819  1.1.1.3.4.1  christos   /* x * y = (2^65+1)/2^68 = 0.001...000[1]000..., and should round to 0.25 */
    820  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDN);
    821  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    822  1.1.1.3.4.1  christos 
    823  1.1.1.3.4.1  christos   /* test for the generic code */
    824  1.1.1.3.4.1  christos   mpfr_set_prec (x, 54);
    825  1.1.1.3.4.1  christos   mpfr_set_prec (y, 55);
    826  1.1.1.3.4.1  christos   mpfr_set_prec (z, 53);
    827  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.101000001E-1");
    828  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.110011000010100101111000011011000111011000001E-1");
    829  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDN);
    830  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    831  1.1.1.3.4.1  christos 
    832  1.1.1.3.4.1  christos   /* another test for the generic code */
    833  1.1.1.3.4.1  christos   mpfr_set_prec (x, 66);
    834  1.1.1.3.4.1  christos   mpfr_set_prec (y, 67);
    835  1.1.1.3.4.1  christos   mpfr_set_prec (z, 65);
    836  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.101101000010010110100001E-1");
    837  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.101101011110010101011010001111111001100001E-1");
    838  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDN);
    839  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -2) == 0);
    840  1.1.1.3.4.1  christos 
    841  1.1.1.3.4.1  christos   mpfr_clear (x);
    842  1.1.1.3.4.1  christos   mpfr_clear (y);
    843  1.1.1.3.4.1  christos   mpfr_clear (z);
    844  1.1.1.3.4.1  christos   set_emin (emin);
    845  1.1.1.3.4.1  christos }
    846  1.1.1.3.4.1  christos 
    847  1.1.1.3.4.1  christos /* test for case in mpfr_mul_1() where:
    848  1.1.1.3.4.1  christos    ax = __gmpfr_emin - 1
    849  1.1.1.3.4.1  christos    ap[0] == ~mask
    850  1.1.1.3.4.1  christos    rnd_mode = MPFR_RNDZ.
    851  1.1.1.3.4.1  christos    Whatever the values of rb and sb, we should round to zero (underflow). */
    852  1.1.1.3.4.1  christos static void
    853  1.1.1.3.4.1  christos bug20161209a (void)
    854  1.1.1.3.4.1  christos {
    855  1.1.1.3.4.1  christos   mpfr_exp_t emin;
    856  1.1.1.3.4.1  christos   mpfr_t x, y, z;
    857  1.1.1.3.4.1  christos 
    858  1.1.1.3.4.1  christos   emin = mpfr_get_emin ();
    859  1.1.1.3.4.1  christos   set_emin (-1);
    860  1.1.1.3.4.1  christos 
    861  1.1.1.3.4.1  christos   mpfr_init2 (x, 53);
    862  1.1.1.3.4.1  christos   mpfr_init2 (y, 53);
    863  1.1.1.3.4.1  christos   mpfr_init2 (z, 53);
    864  1.1.1.3.4.1  christos 
    865  1.1.1.3.4.1  christos   /* case rb = sb = 0 */
    866  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.11010010100110000110110011111E-1");
    867  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.1001101110011000110100001");
    868  1.1.1.3.4.1  christos   /* x = 441650591/2^30, y = 20394401/2^25 */
    869  1.1.1.3.4.1  christos   /* x * y = (2^53-1)/2^55 = 0.00111...111[0]000..., and should round to 0 */
    870  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDZ);
    871  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_zero_p (z));
    872  1.1.1.3.4.1  christos 
    873  1.1.1.3.4.1  christos   /* case rb = 1, sb = 0 */
    874  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.111111111000000000000000000111111111E-1");
    875  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.1000000001000000001");
    876  1.1.1.3.4.1  christos   /* x = 68585259519/2^37, y = 262657/2^19 */
    877  1.1.1.3.4.1  christos   /* x * y = (2^54-1)/2^56 = 0.00111...111[1]000..., and should round to 0 */
    878  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDZ);
    879  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_zero_p (z));
    880  1.1.1.3.4.1  christos 
    881  1.1.1.3.4.1  christos   /* case rb = 0, sb = 1 */
    882  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.110010011001011110001100100001000001E-1");
    883  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.10100010100010111101");
    884  1.1.1.3.4.1  christos   /* x = 541144371852^37, y = 665789/2^20 */
    885  1.1.1.3.4.1  christos   /* x * y = (2^55-3)/2^57 = 0.00111...111[0]100..., and should round to 0 */
    886  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDZ);
    887  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_zero_p (z));
    888  1.1.1.3.4.1  christos 
    889  1.1.1.3.4.1  christos   /* case rb = sb = 1 */
    890  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.10100110001001001010001111110010100111E-1");
    891  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.110001010011101001");
    892  1.1.1.3.4.1  christos   /* x = 178394823847/2^39, y = 201961/2^18 */
    893  1.1.1.3.4.1  christos   /* x * y = (2^55-1)/2^57 = 0.00111...111[1]100..., and should round to 0 */
    894  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDZ);
    895  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_zero_p (z));
    896  1.1.1.3.4.1  christos 
    897  1.1.1.3.4.1  christos   /* similar test for mpfr_mul_2 (we only check rb = sb = 1 here) */
    898  1.1.1.3.4.1  christos   mpfr_set_prec (x, 65);
    899  1.1.1.3.4.1  christos   mpfr_set_prec (y, 65);
    900  1.1.1.3.4.1  christos   mpfr_set_prec (z, 65);
    901  1.1.1.3.4.1  christos   /* 2^67-1 = 193707721 * 761838257287 */
    902  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.1011100010111011111011001001E-1");
    903  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.1011000101100001000110010100010010000111");
    904  1.1.1.3.4.1  christos   mpfr_mul (z, x, y, MPFR_RNDZ);
    905  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_zero_p (z));
    906  1.1.1.3.4.1  christos 
    907  1.1.1.3.4.1  christos   mpfr_clear (x);
    908  1.1.1.3.4.1  christos   mpfr_clear (y);
    909  1.1.1.3.4.1  christos   mpfr_clear (z);
    910  1.1.1.3.4.1  christos   set_emin (emin);
    911  1.1.1.3.4.1  christos }
    912  1.1.1.3.4.1  christos 
    913  1.1.1.3.4.1  christos /* bug for RNDF */
    914  1.1.1.3.4.1  christos static void
    915  1.1.1.3.4.1  christos bug20170602 (void)
    916  1.1.1.3.4.1  christos {
    917  1.1.1.3.4.1  christos   mpfr_t x, u, y, yd, yu;
    918  1.1.1.3.4.1  christos 
    919  1.1.1.3.4.1  christos   mpfr_init2 (x, 493);
    920  1.1.1.3.4.1  christos   mpfr_init2 (u, 493);
    921  1.1.1.3.4.1  christos   mpfr_init2 (y, 503);
    922  1.1.1.3.4.1  christos   mpfr_init2 (yd, 503);
    923  1.1.1.3.4.1  christos   mpfr_init2 (yu, 503);
    924  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.1111100000000000001111111110000000001111111111111000000000000000000011111111111111111111111000000000000000000001111111111111111111111111111111111111111000000000011111111111111111111000000000011111111111111000000000000001110000000000000000000000000000000000000000011111111111110011111111111100000000000000011111111111111111110000000011111111111111111110011111111111110000000000001111111111111111000000000000000000000000000000000000111111111111111111111111111111011111111111111111111111111111100E44");
    925  1.1.1.3.4.1  christos   mpfr_set_str_binary (u, "0.1110000000000000001111111111111111111111111111111111111000000000111111111111111111111111111111000000000000000000001111111000000000000000011111111111111111111111111111111111111111111111111111111000000000000000011111111111111000000011111111111111111110000000000000001111111111111111111111111111111111111110000000000001111111111111111111111111111111111111000000000000000000000000000000000001111111111111000000000000000000001111111111100000000000000011111111111111111111111111111111111111111111111E35");
    926  1.1.1.3.4.1  christos   mpfr_mul (y, x, u, MPFR_RNDF);
    927  1.1.1.3.4.1  christos   mpfr_mul (yd, x, u, MPFR_RNDD);
    928  1.1.1.3.4.1  christos   mpfr_mul (yu, x, u, MPFR_RNDU);
    929  1.1.1.3.4.1  christos   if (mpfr_cmp (y, yd) != 0 && mpfr_cmp (y, yu) != 0)
    930  1.1.1.3.4.1  christos     {
    931  1.1.1.3.4.1  christos       printf ("RNDF is neither RNDD nor RNDU\n");
    932  1.1.1.3.4.1  christos       printf ("x"); mpfr_dump (x);
    933  1.1.1.3.4.1  christos       printf ("u"); mpfr_dump (u);
    934  1.1.1.3.4.1  christos       printf ("y(RNDF)="); mpfr_dump (y);
    935  1.1.1.3.4.1  christos       printf ("y(RNDD)="); mpfr_dump (yd);
    936  1.1.1.3.4.1  christos       printf ("y(RNDU)="); mpfr_dump (yu);
    937  1.1.1.3.4.1  christos       exit (1);
    938  1.1.1.3.4.1  christos     }
    939  1.1.1.3.4.1  christos   mpfr_clear (x);
    940  1.1.1.3.4.1  christos   mpfr_clear (u);
    941  1.1.1.3.4.1  christos   mpfr_clear (y);
    942  1.1.1.3.4.1  christos   mpfr_clear (yd);
    943  1.1.1.3.4.1  christos   mpfr_clear (yu);
    944  1.1.1.3.4.1  christos }
    945  1.1.1.3.4.1  christos 
    946  1.1.1.3.4.1  christos /* Test for 1 to 3 limbs. */
    947  1.1.1.3.4.1  christos static void
    948  1.1.1.3.4.1  christos small_prec (void)
    949  1.1.1.3.4.1  christos {
    950  1.1.1.3.4.1  christos   mpfr_exp_t emin, emax;
    951  1.1.1.3.4.1  christos   mpfr_t x, y, z1, z2, zz;
    952  1.1.1.3.4.1  christos   int xq, yq, zq;
    953  1.1.1.3.4.1  christos   mpfr_rnd_t rnd;
    954  1.1.1.3.4.1  christos   mpfr_flags_t flags1, flags2;
    955  1.1.1.3.4.1  christos   int inex1, inex2;
    956  1.1.1.3.4.1  christos   int i, j, r, s, ediff;
    957  1.1.1.3.4.1  christos 
    958  1.1.1.3.4.1  christos   emin = mpfr_get_emin ();
    959  1.1.1.3.4.1  christos   emax = mpfr_get_emax ();
    960  1.1.1.3.4.1  christos 
    961  1.1.1.3.4.1  christos   /* The mpfr_mul implementation doesn't extend the exponent range,
    962  1.1.1.3.4.1  christos      so that it is OK to reduce it here for the test to make sure that
    963  1.1.1.3.4.1  christos      mpfr_mul_2si can be used. */
    964  1.1.1.3.4.1  christos   set_emin (-1000);
    965  1.1.1.3.4.1  christos   set_emax (1000);
    966  1.1.1.3.4.1  christos 
    967  1.1.1.3.4.1  christos   mpfr_inits2 (3 * GMP_NUMB_BITS, x, y, z1, z2, (mpfr_ptr) 0);
    968  1.1.1.3.4.1  christos   mpfr_init2 (zz, 6 * GMP_NUMB_BITS);
    969  1.1.1.3.4.1  christos   for (i = 0; i < 3; i++)
    970  1.1.1.3.4.1  christos     for (j = 0; j < 10000; j++)
    971  1.1.1.3.4.1  christos       {
    972  1.1.1.3.4.1  christos         xq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
    973  1.1.1.3.4.1  christos         mpfr_set_prec (x, xq);
    974  1.1.1.3.4.1  christos         yq = i * GMP_NUMB_BITS + 1 + randlimb () % GMP_NUMB_BITS;
    975  1.1.1.3.4.1  christos         mpfr_set_prec (y, yq);
    976  1.1.1.3.4.1  christos         zq = i * GMP_NUMB_BITS + 1 + randlimb () % (GMP_NUMB_BITS-1);
    977  1.1.1.3.4.1  christos         mpfr_set_prec (z1, zq);
    978  1.1.1.3.4.1  christos         mpfr_set_prec (z2, zq);
    979  1.1.1.3.4.1  christos         s = r = randlimb () & 0x7f;
    980  1.1.1.3.4.1  christos         do mpfr_urandomb (x, RANDS); while (mpfr_zero_p (x));
    981  1.1.1.3.4.1  christos         if (s & 1)
    982  1.1.1.3.4.1  christos           mpfr_neg (x, x, MPFR_RNDN);
    983  1.1.1.3.4.1  christos         s >>= 1;
    984  1.1.1.3.4.1  christos         if (s & 1)
    985  1.1.1.3.4.1  christos           {
    986  1.1.1.3.4.1  christos             do mpfr_urandomb (y, RANDS); while (mpfr_zero_p (y));
    987  1.1.1.3.4.1  christos           }
    988  1.1.1.3.4.1  christos         else
    989  1.1.1.3.4.1  christos           {
    990  1.1.1.3.4.1  christos             mpfr_ui_div (y, 1, x, MPFR_RNDN);
    991  1.1.1.3.4.1  christos             mpfr_set_exp (y, 0);
    992  1.1.1.3.4.1  christos           }
    993  1.1.1.3.4.1  christos         s >>= 1;
    994  1.1.1.3.4.1  christos         if (s & 1)
    995  1.1.1.3.4.1  christos           mpfr_neg (y, y, MPFR_RNDN);
    996  1.1.1.3.4.1  christos         s >>= 1;
    997  1.1.1.3.4.1  christos         rnd = RND_RAND_NO_RNDF ();
    998  1.1.1.3.4.1  christos         inex1 = mpfr_mul (zz, x, y, MPFR_RNDN);
    999  1.1.1.3.4.1  christos         MPFR_ASSERTN (inex1 == 0);
   1000  1.1.1.3.4.1  christos         if (s == 0)
   1001  1.1.1.3.4.1  christos           {
   1002  1.1.1.3.4.1  christos             ediff = __gmpfr_emin - MPFR_EXP (x);
   1003  1.1.1.3.4.1  christos             mpfr_set_exp (x, __gmpfr_emin);
   1004  1.1.1.3.4.1  christos           }
   1005  1.1.1.3.4.1  christos         else if (s == 1)
   1006  1.1.1.3.4.1  christos           {
   1007  1.1.1.3.4.1  christos             ediff = __gmpfr_emax - MPFR_EXP (x) + 1;
   1008  1.1.1.3.4.1  christos             mpfr_set_exp (x, __gmpfr_emax);
   1009  1.1.1.3.4.1  christos             mpfr_mul_2ui (y, y, 1, MPFR_RNDN);
   1010  1.1.1.3.4.1  christos           }
   1011  1.1.1.3.4.1  christos         else
   1012  1.1.1.3.4.1  christos           ediff = 0;
   1013  1.1.1.3.4.1  christos         mpfr_clear_flags ();
   1014  1.1.1.3.4.1  christos         inex1 = mpfr_mul_2si (z1, zz, ediff, rnd);
   1015  1.1.1.3.4.1  christos         flags1 = __gmpfr_flags;
   1016  1.1.1.3.4.1  christos         mpfr_clear_flags ();
   1017  1.1.1.3.4.1  christos         inex2 = mpfr_mul (z2, x, y, rnd);
   1018  1.1.1.3.4.1  christos         flags2 = __gmpfr_flags;
   1019  1.1.1.3.4.1  christos         if (!(mpfr_equal_p (z1, z2) &&
   1020  1.1.1.3.4.1  christos               SAME_SIGN (inex1, inex2) &&
   1021  1.1.1.3.4.1  christos               flags1 == flags2))
   1022  1.1.1.3.4.1  christos           {
   1023  1.1.1.3.4.1  christos             printf ("Error in small_prec on i = %d, j = %d\n", i, j);
   1024  1.1.1.3.4.1  christos             printf ("r = 0x%x, xq = %d, yq = %d, zq = %d, rnd = %s\n",
   1025  1.1.1.3.4.1  christos                     r, xq, yq, zq, mpfr_print_rnd_mode (rnd));
   1026  1.1.1.3.4.1  christos             printf ("x = ");
   1027  1.1.1.3.4.1  christos             mpfr_dump (x);
   1028  1.1.1.3.4.1  christos             printf ("y = ");
   1029  1.1.1.3.4.1  christos             mpfr_dump (y);
   1030  1.1.1.3.4.1  christos             printf ("ediff = %d\n", ediff);
   1031  1.1.1.3.4.1  christos             printf ("zz = ");
   1032  1.1.1.3.4.1  christos             mpfr_dump (zz);
   1033  1.1.1.3.4.1  christos             printf ("Expected ");
   1034  1.1.1.3.4.1  christos             mpfr_dump (z1);
   1035  1.1.1.3.4.1  christos             printf ("with inex = %d and flags =", inex1);
   1036  1.1.1.3.4.1  christos             flags_out (flags1);
   1037  1.1.1.3.4.1  christos             printf ("Got      ");
   1038  1.1.1.3.4.1  christos             mpfr_dump (z2);
   1039  1.1.1.3.4.1  christos             printf ("with inex = %d and flags =", inex2);
   1040  1.1.1.3.4.1  christos             flags_out (flags2);
   1041  1.1.1.3.4.1  christos             exit (1);
   1042  1.1.1.3.4.1  christos           }
   1043  1.1.1.3.4.1  christos       }
   1044  1.1.1.3.4.1  christos   mpfr_clears (x, y, z1, z2, zz, (mpfr_ptr) 0);
   1045  1.1.1.3.4.1  christos 
   1046  1.1.1.3.4.1  christos   set_emin (emin);
   1047  1.1.1.3.4.1  christos   set_emax (emax);
   1048  1.1.1.3.4.1  christos }
   1049  1.1.1.3.4.1  christos 
   1050  1.1.1.3.4.1  christos /* check ax < __gmpfr_emin with rnd_mode == MPFR_RNDN, rb = 0 and sb <> 0 */
   1051  1.1.1.3.4.1  christos static void
   1052  1.1.1.3.4.1  christos test_underflow2 (void)
   1053  1.1.1.3.4.1  christos {
   1054  1.1.1.3.4.1  christos   mpfr_t x, y, z;
   1055  1.1.1.3.4.1  christos   mpfr_exp_t emin;
   1056  1.1.1.3.4.1  christos 
   1057  1.1.1.3.4.1  christos   emin = mpfr_get_emin ();
   1058  1.1.1.3.4.1  christos   mpfr_set_emin (0);
   1059  1.1.1.3.4.1  christos 
   1060  1.1.1.3.4.1  christos   mpfr_init2 (x, 24);
   1061  1.1.1.3.4.1  christos   mpfr_init2 (y, 24);
   1062  1.1.1.3.4.1  christos   mpfr_init2 (z, 24);
   1063  1.1.1.3.4.1  christos 
   1064  1.1.1.3.4.1  christos   mpfr_set_ui_2exp (x, 12913, -14, MPFR_RNDN);
   1065  1.1.1.3.4.1  christos   mpfr_set_ui_2exp (y, 5197,  -13, MPFR_RNDN);
   1066  1.1.1.3.4.1  christos   mpfr_clear_underflow ();
   1067  1.1.1.3.4.1  christos   /* x*y = 0.0111111111111111111111111[01] thus underflow */
   1068  1.1.1.3.4.1  christos   mpfr_mul (z, y, x, MPFR_RNDN);
   1069  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
   1070  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_underflow_p ());
   1071  1.1.1.3.4.1  christos 
   1072  1.1.1.3.4.1  christos   mpfr_set_prec (x, 65);
   1073  1.1.1.3.4.1  christos   mpfr_set_prec (y, 65);
   1074  1.1.1.3.4.1  christos   mpfr_set_prec (z, 65);
   1075  1.1.1.3.4.1  christos   mpfr_set_str_binary (x, "0.10011101000110111011111100101001111");
   1076  1.1.1.3.4.1  christos   mpfr_set_str_binary (y, "0.110100001001000111000011011110011");
   1077  1.1.1.3.4.1  christos   mpfr_clear_underflow ();
   1078  1.1.1.3.4.1  christos   /* x*y = 0.011111111111111111111111111111111111111111111111111111111111111111[01] thus underflow */
   1079  1.1.1.3.4.1  christos   mpfr_mul (z, y, x, MPFR_RNDN);
   1080  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_cmp_ui_2exp (z, 1, -1) == 0);
   1081  1.1.1.3.4.1  christos   MPFR_ASSERTN(mpfr_underflow_p ());
   1082  1.1.1.3.4.1  christos 
   1083  1.1.1.3.4.1  christos   mpfr_clear (y);
   1084  1.1.1.3.4.1  christos   mpfr_clear (x);
   1085  1.1.1.3.4.1  christos   mpfr_clear (z);
   1086  1.1.1.3.4.1  christos 
   1087  1.1.1.3.4.1  christos   mpfr_set_emin (emin);
   1088  1.1.1.3.4.1  christos }
   1089  1.1.1.3.4.1  christos 
   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.3.4.1  christos   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.3.4.1  christos   small_prec ();
   1132          1.1       mrg 
   1133          1.1       mrg   check_regression ();
   1134  1.1.1.3.4.1  christos   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.3.4.1  christos   bug20161209 ();
   1141  1.1.1.3.4.1  christos   bug20161209a ();
   1142  1.1.1.3.4.1  christos   bug20170602 ();
   1143  1.1.1.3.4.1  christos   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