github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mini-gmp/tests/t-invert.c (about)

     1  /*
     2  
     3  Copyright 2012, 2016 Free Software Foundation, Inc.
     4  
     5  This file is part of the GNU MP Library test suite.
     6  
     7  The GNU MP Library test suite is free software; you can redistribute it
     8  and/or modify it under the terms of the GNU General Public License as
     9  published by the Free Software Foundation; either version 3 of the License,
    10  or (at your option) any later version.
    11  
    12  The GNU MP Library test suite is distributed in the hope that it will be
    13  useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
    14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General
    15  Public License for more details.
    16  
    17  You should have received a copy of the GNU General Public License along with
    18  the GNU MP Library test suite.  If not, see https://www.gnu.org/licenses/.  */
    19  
    20  #include <assert.h>
    21  #include <limits.h>
    22  #include <stdlib.h>
    23  #include <stdio.h>
    24  
    25  #include "testutils.h"
    26  
    27  #define GMP_LIMB_BITS (sizeof(mp_limb_t) * CHAR_BIT)
    28  
    29  #define COUNT 10000
    30  
    31  static void
    32  test_2by1(const mpz_t u)
    33  {
    34    mpz_t m, p, t;
    35  
    36    mpz_init (m);
    37    mpz_init (p);
    38    mpz_init (t);
    39  
    40    assert (mpz_size (u) == 1);
    41  
    42    mpz_set_ui (m, mpn_invert_limb (u->_mp_d[0]));
    43    mpz_setbit (m, GMP_LIMB_BITS);
    44  
    45    mpz_mul (p, m, u);
    46  
    47    mpz_set_ui (t, 0);
    48    mpz_setbit (t, 2* GMP_LIMB_BITS);
    49    mpz_sub (t, t, p);
    50  
    51    /* Should have 0 < B^2 - m u <= u */
    52    if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
    53      {
    54        fprintf (stderr, "mpn_invert_limb failed:\n");
    55        dump ("u", u);
    56        dump ("m", m);
    57        dump ("p", p);
    58        dump ("t", t);
    59        abort ();
    60      }
    61    mpz_clear (m);
    62    mpz_clear (p);
    63    mpz_clear (t);
    64  }
    65  
    66  static void
    67  test_3by2(const mpz_t u)
    68  {
    69    mpz_t m, p, t;
    70  
    71    mpz_init (m);
    72    mpz_init (p);
    73    mpz_init (t);
    74  
    75    assert (mpz_size (u) == 2);
    76  
    77    mpz_set_ui (m, mpn_invert_3by2 (u->_mp_d[1], u[0]._mp_d[0]));
    78  
    79    mpz_setbit (m, GMP_LIMB_BITS);
    80  
    81    mpz_mul (p, m, u);
    82  
    83    mpz_set_ui (t, 0);
    84    mpz_setbit (t, 3 * GMP_LIMB_BITS);
    85    mpz_sub (t, t, p);
    86  
    87    /* Should have 0 < B^3 - m u <= u */
    88    if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
    89      {
    90        fprintf (stderr, "mpn_invert_3by2 failed:\n");
    91        dump ("u", u);
    92        dump ("m", m);
    93        dump ("p", p);
    94        dump ("t", t);
    95        abort ();
    96      }
    97    mpz_clear (m);
    98    mpz_clear (p);
    99    mpz_clear (t);
   100  }
   101  
   102  void
   103  testmain (int argc, char **argv)
   104  {
   105    unsigned i;
   106    mpz_t u, m, p, t;
   107  
   108    mpz_init (u);
   109    mpz_init (m);
   110    mpz_init (p);
   111    mpz_init (t);
   112  
   113    /* These values trigger 32-bit overflow of ql in mpn_invert_3by2. */
   114    if (GMP_LIMB_BITS == 64)
   115      {
   116        mpz_set_str (u, "80007fff3ffe0000", 16);
   117        test_2by1 (u);
   118        mpz_set_str (u, "80007fff3ffe000000000000000003ff", 16);
   119        test_3by2 (u);
   120      }
   121  
   122    for (i = 0; i < COUNT; i++)
   123      {
   124        mini_urandomb (u, GMP_LIMB_BITS);
   125        mpz_setbit (u, GMP_LIMB_BITS -1);
   126  
   127        test_2by1 (u);
   128      }
   129  
   130    for (i = 0; i < COUNT; i++)
   131      {
   132        mini_urandomb (u, 2*GMP_LIMB_BITS);
   133        mpz_setbit (u, 2*GMP_LIMB_BITS -1);
   134  
   135        test_3by2 (u);
   136      }
   137  
   138    mpz_clear (u);
   139  }