"Fossies" - the Fresh Open Source Software Archive

Member "mpfr-4.0.2/tests/tcmp2.c" (7 Jan 2019, 12280 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 "tcmp2.c": 4.0.1_vs_4.0.2.

    1 /* Test file for mpfr_cmp2.
    2 
    3 Copyright 1999-2003, 2006-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 /* set bit n of x to b, where bit 0 is the most significant one */
   26 static void
   27 set_bit (mpfr_t x, unsigned int n, int b)
   28 {
   29   unsigned l;
   30   mp_size_t xn;
   31 
   32   xn = (MPFR_PREC(x) - 1) / mp_bits_per_limb;
   33   l = n / mp_bits_per_limb;
   34   n %= mp_bits_per_limb;
   35   n = mp_bits_per_limb - 1 - n;
   36   if (b)
   37     MPFR_MANT(x)[xn - l] |= MPFR_LIMB_ONE << n;
   38   else
   39     MPFR_MANT(x)[xn - l] &= ~(MPFR_LIMB_ONE << n);
   40 }
   41 
   42 /* check that for x = 1.u 1 v 0^k low(x)
   43                   y = 1.u 0 v 1^k low(y)
   44    mpfr_cmp2 (x, y) returns 1 + |u| + |v| + k for low(x) >= low(y),
   45                         and 1 + |u| + |v| + k + 1 otherwise */
   46 static void
   47 worst_cases (void)
   48 {
   49   mpfr_t x, y;
   50   unsigned int i, j, k, b, expected;
   51   mpfr_prec_t l;
   52 
   53   mpfr_init2 (x, 200);
   54   mpfr_init2 (y, 200);
   55 
   56   mpfr_set_ui (y, 1, MPFR_RNDN);
   57   for (i = 1; i < MPFR_PREC(x); i++)
   58     {
   59       mpfr_set_ui (x, 1, MPFR_RNDN);
   60       mpfr_div_2exp (y, y, 1, MPFR_RNDN); /* y = 1/2^i */
   61 
   62       l = 0;
   63       if (mpfr_cmp2 (x, y, &l) <= 0 || l != 1)
   64         {
   65           printf ("Error in mpfr_cmp2:\nx=");
   66           mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
   67           printf ("\ny=");
   68           mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
   69           printf ("\ngot %lu instead of 1\n", (unsigned long) l);
   70           exit (1);
   71         }
   72 
   73       mpfr_add (x, x, y, MPFR_RNDN); /* x = 1 + 1/2^i */
   74       l = 0;
   75       if (mpfr_cmp2 (x, y, &l) <= 0 || l != 0)
   76         {
   77           printf ("Error in mpfr_cmp2:\nx=");
   78           mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
   79           printf ("\ny=");
   80           mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
   81           printf ("\ngot %lu instead of 0\n", (unsigned long) l);
   82           exit (1);
   83         }
   84     }
   85 
   86   for (i = 0; i < 64; i++) /* |u| = i */
   87     {
   88       mpfr_urandomb (x, RANDS);
   89       mpfr_set (y, x, MPFR_RNDN);
   90       set_bit (x, i + 1, 1);
   91       set_bit (y, i + 1, 0);
   92       for (j = 0; j < 64; j++) /* |v| = j */
   93         {
   94           b = randlimb () % 2;
   95           set_bit (x, i + j + 2, b);
   96           set_bit (y, i + j + 2, b);
   97           for (k = 0; k < 64; k++)
   98             {
   99               if (k)
  100                 set_bit (x, i + j + k + 1, 0);
  101               if (k)
  102                 set_bit (y, i + j + k + 1, 1);
  103               set_bit (x, i + j + k + 2, 1);
  104               set_bit (y, i + j + k + 2, 0);
  105               l = 0; mpfr_cmp2 (x, y, &l);
  106               expected = i + j + k + 1;
  107               if (l != expected)
  108                 {
  109                   printf ("Error in mpfr_cmp2:\nx=");
  110                   mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
  111                   printf ("\ny=");
  112                   mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
  113                   printf ("\ngot %lu instead of %u\n",
  114                           (unsigned long) l, expected);
  115                   exit (1);
  116                 }
  117               set_bit (x, i + j + k + 2, 0);
  118               set_bit (x, i + j + k + 3, 0);
  119               set_bit (y, i + j + k + 3, 1);
  120               l = 0; mpfr_cmp2 (x, y, &l);
  121               expected = i + j + k + 2;
  122               if (l != expected)
  123                 {
  124                   printf ("Error in mpfr_cmp2:\nx=");
  125                   mpfr_out_str (stdout, 2, 0, x, MPFR_RNDN);
  126                   printf ("\ny=");
  127                   mpfr_out_str (stdout, 2, 0, y, MPFR_RNDN);
  128                   printf ("\ngot %lu instead of %u\n",
  129                           (unsigned long) l, expected);
  130                   exit (1);
  131                 }
  132             }
  133         }
  134     }
  135 
  136   mpfr_clear (x);
  137   mpfr_clear (y);
  138 }
  139 
  140 static void
  141 tcmp2 (double x, double y, int i)
  142 {
  143   mpfr_t xx, yy;
  144   mpfr_prec_t j;
  145 
  146   if (i == -1)
  147     {
  148       if (x == y)
  149         i = 53;
  150       else
  151         i = (int) (__gmpfr_floor_log2 (x) - __gmpfr_floor_log2 (x - y));
  152     }
  153   mpfr_init2(xx, 53); mpfr_init2(yy, 53);
  154   mpfr_set_d (xx, x, MPFR_RNDN);
  155   mpfr_set_d (yy, y, MPFR_RNDN);
  156   j = 0;
  157   if (mpfr_cmp2 (xx, yy, &j) == 0)
  158     {
  159       if (x != y)
  160         {
  161           printf ("Error in mpfr_cmp2 for\nx=");
  162           mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN);
  163           printf ("\ny=");
  164           mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN);
  165           printf ("\ngot sign 0 for x != y\n");
  166           exit (1);
  167         }
  168     }
  169   else if (j != (unsigned) i)
  170     {
  171       printf ("Error in mpfr_cmp2 for\nx=");
  172       mpfr_out_str (stdout, 2, 0, xx, MPFR_RNDN);
  173       printf ("\ny=");
  174       mpfr_out_str (stdout, 2, 0, yy, MPFR_RNDN);
  175       printf ("\ngot %lu instead of %d\n", (unsigned long) j, i);
  176       exit (1);
  177     }
  178   mpfr_clear(xx); mpfr_clear(yy);
  179 }
  180 
  181 static void
  182 special (void)
  183 {
  184   mpfr_t x, y;
  185   mpfr_prec_t j;
  186 
  187   mpfr_init (x); mpfr_init (y);
  188 
  189   /* bug found by Nathalie Revol, 21 March 2001 */
  190   mpfr_set_prec (x, 65);
  191   mpfr_set_prec (y, 65);
  192   mpfr_set_str_binary (x, "0.10000000000000000000000000000000000001110010010110100110011110000E1");
  193   mpfr_set_str_binary (y, "0.11100100101101001100111011111111110001101001000011101001001010010E-35");
  194   j = 0;
  195   if (mpfr_cmp2 (x, y, &j) <= 0 || j != 1)
  196     {
  197       printf ("Error in mpfr_cmp2:\n");
  198       printf ("x=");
  199       mpfr_dump (x);
  200       printf ("y=");
  201       mpfr_dump (y);
  202       printf ("got %lu, expected 1\n", (unsigned long) j);
  203       exit (1);
  204     }
  205 
  206   mpfr_set_prec(x, 127); mpfr_set_prec(y, 127);
  207   mpfr_set_str_binary(x, "0.1011010000110111111000000101011110110001000101101011011110010010011110010000101101000010011001100110010000000010110000101000101E6");
  208   mpfr_set_str_binary(y, "0.1011010000110111111000000101011011111100011101000011001111000010100010100110110100110010011001100110010000110010010110000010110E6");
  209   j = 0;
  210   if (mpfr_cmp2(x, y, &j) <= 0 || j != 32)
  211     {
  212       printf ("Error in mpfr_cmp2:\n");
  213       printf ("x="); mpfr_dump (x);
  214       printf ("y="); mpfr_dump (y);
  215       printf ("got %lu, expected 32\n", (unsigned long) j);
  216       exit (1);
  217     }
  218 
  219   mpfr_set_prec (x, 128); mpfr_set_prec (y, 239);
  220   mpfr_set_str_binary (x, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100100000000000E167");
  221   mpfr_set_str_binary (y, "0.10001000110110000111011000101011111100110010010011001101000011111010010110001000000010100110100111111011011010101100011111111111111111111111111111111111111111111111011111100101011100011001101000100111000010000000000101100110000111111000101E167");
  222   j = 0;
  223   if (mpfr_cmp2(x, y, &j) <= 0 || j != 164)
  224     {
  225       printf ("Error in mpfr_cmp2:\n");
  226       printf ("x="); mpfr_dump (x);
  227       printf ("y="); mpfr_dump (y);
  228       printf ("got %lu, expected 164\n", (unsigned long) j);
  229       exit (1);
  230     }
  231 
  232   /* bug found by Nathalie Revol, 29 March 2001 */
  233   mpfr_set_prec (x, 130); mpfr_set_prec (y, 130);
  234   mpfr_set_str_binary (x, "0.1100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000E2");
  235   mpfr_set_str_binary (y, "0.1011111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111100E2");
  236   j = 0;
  237   if (mpfr_cmp2(x, y, &j) <= 0 || j != 127)
  238     {
  239       printf ("Error in mpfr_cmp2:\n");
  240       printf ("x="); mpfr_dump (x);
  241       printf ("y="); mpfr_dump (y);
  242       printf ("got %lu, expected 127\n", (unsigned long) j);
  243       exit (1);
  244     }
  245 
  246   /* bug found by Nathalie Revol, 2 Apr 2001 */
  247   mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
  248   mpfr_set_ui (x, 5, MPFR_RNDN);
  249   mpfr_set_str_binary (y, "0.10011111111111111111111111111111111111111111111111111111111111101E3");
  250   j = 0;
  251   if (mpfr_cmp2(x, y, &j) <= 0 || j != 63)
  252     {
  253       printf ("Error in mpfr_cmp2:\n");
  254       printf ("x="); mpfr_dump (x);
  255       printf ("y="); mpfr_dump (y);
  256       printf ("got %lu, expected 63\n", (unsigned long) j);
  257       exit (1);
  258     }
  259 
  260   /* bug found by Nathalie Revol, 2 Apr 2001 */
  261   mpfr_set_prec (x, 65); mpfr_set_prec (y, 65);
  262   mpfr_set_str_binary (x, "0.10011011111000101001110000000000000000000000000000000000000000000E-69");
  263   mpfr_set_str_binary (y, "0.10011011111000101001101111111111111111111111111111111111111111101E-69");
  264   j = 0;
  265   if (mpfr_cmp2(x, y, &j) <= 0 || j != 63)
  266     {
  267       printf ("Error in mpfr_cmp2:\n");
  268       printf ("x="); mpfr_dump (x);
  269       printf ("y="); mpfr_dump (y);
  270       printf ("got %lu, expected 63\n", (unsigned long) j);
  271       exit (1);
  272     }
  273 
  274   mpfr_set_prec (x, 65);
  275   mpfr_set_prec (y, 32);
  276   mpfr_set_str_binary (x, "1.1110111011110001110111011111111111101000011001011100101100101101");
  277   mpfr_set_str_binary (y, "0.11101110111100011101110111111111");
  278   if (mpfr_cmp2 (x, y, &j) <= 0 || j != 0)
  279     {
  280       printf ("Error in mpfr_cmp2 (1)\n");
  281       exit (1);
  282     }
  283 
  284   mpfr_set_prec (x, 2 * GMP_NUMB_BITS);
  285   mpfr_set_prec (y, GMP_NUMB_BITS);
  286   mpfr_set_ui (x, 1, MPFR_RNDN); /* x = 1 */
  287   mpfr_set_ui (y, 1, MPFR_RNDN);
  288   mpfr_nextbelow (y);            /* y = 1 - 2^(-GMP_NUMB_BITS) */
  289   mpfr_cmp2 (x, y, &j);
  290   if (mpfr_cmp2 (x, y, &j) <= 0 || j != GMP_NUMB_BITS)
  291     {
  292       printf ("Error in mpfr_cmp2 (2)\n");
  293       exit (1);
  294     }
  295 
  296   mpfr_clear (x);
  297   mpfr_clear (y);
  298 }
  299 
  300 /* Compare (m,kx) and (m,ky), where (m,k) means m fixed limbs followed by
  301    k zero limbs. */
  302 static void
  303 test_equal (void)
  304 {
  305   mpfr_t w, x, y;
  306   int m, kx, ky, inex;
  307   mpfr_prec_t j;
  308 
  309   for (m = 1; m <= 4; m++)
  310     {
  311       mpfr_init2 (w, m * GMP_NUMB_BITS);
  312       for (kx = 0; kx <= 4; kx++)
  313         for (ky = 0; ky <= 4; ky++)
  314           {
  315             do mpfr_urandomb (w, RANDS); while (mpfr_zero_p (w));
  316             mpfr_init2 (x, (m + kx) * GMP_NUMB_BITS
  317                         - (kx == 0 ? 0 : randlimb () % GMP_NUMB_BITS));
  318             mpfr_init2 (y, (m + ky) * GMP_NUMB_BITS
  319                         - (ky == 0 ? 0 : randlimb () % GMP_NUMB_BITS));
  320             inex = mpfr_set (x, w, MPFR_RNDN);
  321             MPFR_ASSERTN (inex == 0);
  322             inex = mpfr_set (y, w, MPFR_RNDN);
  323             MPFR_ASSERTN (inex == 0);
  324             MPFR_ASSERTN (mpfr_equal_p (x, y));
  325             if (randlimb () & 1)
  326               mpfr_neg (x, x, MPFR_RNDN);
  327             if (randlimb () & 1)
  328               mpfr_neg (y, y, MPFR_RNDN);
  329             if (mpfr_cmp2 (x, y, &j) != 0)
  330               {
  331                 printf ("Error in test_equal for m = %d, kx = %d, ky = %d\n",
  332                         m, kx, ky);
  333                 printf ("  x = ");
  334                 mpfr_dump (x);
  335                 printf ("  y = ");
  336                 mpfr_dump (y);
  337                 exit (1);
  338               }
  339             mpfr_clears (x, y, (mpfr_ptr) 0);
  340           }
  341       mpfr_clear (w);
  342     }
  343 }
  344 
  345 int
  346 main (void)
  347 {
  348   int i;
  349   long j;
  350   double x, y, z;
  351 
  352   tests_start_mpfr ();
  353   mpfr_test_init ();
  354 
  355   worst_cases ();
  356   special ();
  357   tcmp2 (5.43885304644369510000e+185, -1.87427265794105340000e-57, 1);
  358   tcmp2 (1.06022698059744327881e+71, 1.05824655795525779205e+71, -1);
  359   tcmp2 (1.0, 1.0, 53);
  360   /* warning: cmp2 does not allow 0 as input */
  361   for (x = 0.5, i = 1; i < 54; i++)
  362     {
  363       tcmp2 (1.0, 1.0-x, i);
  364       x /= 2.0;
  365     }
  366   for (x = 0.5, i = 1; i < 100; i++)
  367     {
  368       tcmp2 (1.0, x, 1);
  369       x /= 2.0;
  370     }
  371   for (j = 0; j < 100000; j++)
  372     {
  373       x = DBL_RAND ();
  374       y = DBL_RAND ();
  375       if (x < y)
  376         {
  377           z = x;
  378           x = y;
  379           y = z;
  380         }
  381       if (y != 0.0)
  382         tcmp2 (x, y, -1);
  383     }
  384 
  385   test_equal ();
  386 
  387   tests_end_mpfr ();
  388 
  389   return 0;
  390 }