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