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    mp_limb_t tl;
      36  
      37    mpz_init (p);
      38  
      39    assert (mpz_size (u) == 1);
      40  
      41    tl = mpn_invert_limb (u->_mp_d[0]);
      42    mpz_roinit_n (t, &tl, 1);
      43    mpz_init_set (m, t);
      44    mpz_setbit (m, GMP_LIMB_BITS);
      45  
      46    mpz_mul (p, m, u);
      47  
      48    mpz_init (t);
      49    mpz_setbit (t, 2* GMP_LIMB_BITS);
      50    mpz_sub (t, t, p);
      51  
      52    /* Should have 0 < B^2 - m u <= u */
      53    if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
      54      {
      55        fprintf (stderr, "mpn_invert_limb failed:\n");
      56        dump ("u", u);
      57        dump ("m", m);
      58        dump ("p", p);
      59        dump ("t", t);
      60        abort ();
      61      }
      62    mpz_clear (m);
      63    mpz_clear (p);
      64    mpz_clear (t);
      65  }
      66  
      67  static void
      68  test_3by2(const mpz_t u)
      69  {
      70    mpz_t m, p, t;
      71    mp_limb_t tl;
      72  
      73    mpz_init (p);
      74  
      75    assert (mpz_size (u) == 2);
      76  
      77    tl = mpn_invert_3by2 (u->_mp_d[1], u->_mp_d[0]);
      78    mpz_roinit_n (t, &tl, 1);
      79    mpz_init_set (m, t);
      80  
      81    mpz_setbit (m, GMP_LIMB_BITS);
      82  
      83    mpz_mul (p, m, u);
      84  
      85    mpz_init (t);
      86    mpz_setbit (t, 3 * GMP_LIMB_BITS);
      87    mpz_sub (t, t, p);
      88  
      89    /* Should have 0 < B^3 - m u <= u */
      90    if (mpz_sgn (t) <= 0 || mpz_cmp (t, u) > 0)
      91      {
      92        fprintf (stderr, "mpn_invert_3by2 failed:\n");
      93        dump ("u", u);
      94        dump ("m", m);
      95        dump ("p", p);
      96        dump ("t", t);
      97        abort ();
      98      }
      99    mpz_clear (m);
     100    mpz_clear (p);
     101    mpz_clear (t);
     102  }
     103  
     104  void
     105  testmain (int argc, char **argv)
     106  {
     107    unsigned i;
     108    mpz_t u, m, p, t;
     109  
     110    mpz_init (u);
     111    mpz_init (m);
     112    mpz_init (p);
     113    mpz_init (t);
     114  
     115    /* These values trigger 32-bit overflow of ql in mpn_invert_3by2. */
     116    if (GMP_LIMB_BITS == 64)
     117      {
     118        mpz_set_str (u, "80007fff3ffe0000", 16);
     119        test_2by1 (u);
     120        mpz_set_str (u, "80007fff3ffe000000000000000003ff", 16);
     121        test_3by2 (u);
     122      }
     123  
     124    for (i = 0; i < COUNT; i++)
     125      {
     126        mini_urandomb (u, GMP_LIMB_BITS);
     127        mpz_setbit (u, GMP_LIMB_BITS -1);
     128  
     129        test_2by1 (u);
     130      }
     131  
     132    for (i = 0; i < COUNT; i++)
     133      {
     134        mini_urandomb (u, 2*GMP_LIMB_BITS);
     135        mpz_setbit (u, 2*GMP_LIMB_BITS -1);
     136  
     137        test_3by2 (u);
     138      }
     139  
     140    mpz_clear (u);
     141  }