"Fossies" - the Fresh Open Source Software Archive

Member "mpfr-4.0.2/tests/tmul.c" (7 Jan 2019, 44549 Bytes) of package /linux/misc/mpfr-4.0.2.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. See also the latest Fossies "Diffs" side-by-side code changes report for "tmul.c": 4.0.1_vs_4.0.2.

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