github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpn/generic/sec_invert.c (about)

     1  /* mpn_sec_invert
     2  
     3     Contributed to the GNU project by Niels Möller
     4  
     5  Copyright 2013 Free Software Foundation, Inc.
     6  
     7  This file is part of the GNU MP Library.
     8  
     9  The GNU MP Library is free software; you can redistribute it and/or modify
    10  it under the terms of either:
    11  
    12    * the GNU Lesser General Public License as published by the Free
    13      Software Foundation; either version 3 of the License, or (at your
    14      option) any later version.
    15  
    16  or
    17  
    18    * the GNU General Public License as published by the Free Software
    19      Foundation; either version 2 of the License, or (at your option) any
    20      later version.
    21  
    22  or both in parallel, as here.
    23  
    24  The GNU MP Library is distributed in the hope that it will be useful, but
    25  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    26  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    27  for more details.
    28  
    29  You should have received copies of the GNU General Public License and the
    30  GNU Lesser General Public License along with the GNU MP Library.  If not,
    31  see https://www.gnu.org/licenses/.  */
    32  
    33  #include "gmp.h"
    34  #include "gmp-impl.h"
    35  
    36  #if 0
    37  /* Currently unused. Should be resurrected once mpn_cnd_neg is
    38     advertised. */
    39  static mp_size_t
    40  mpn_cnd_neg_itch (mp_size_t n)
    41  {
    42    return n;
    43  }
    44  #endif
    45  
    46  /* FIXME: Ought to return carry */
    47  static void
    48  mpn_cnd_neg (int cnd, mp_limb_t *rp, const mp_limb_t *ap, mp_size_t n,
    49  	     mp_ptr scratch)
    50  {
    51    mpn_lshift (scratch, ap, n, 1);
    52    mpn_cnd_sub_n (cnd, rp, ap, scratch, n);
    53  }
    54  
    55  static int
    56  mpn_sec_eq_ui (mp_srcptr ap, mp_size_t n, mp_limb_t b)
    57  {
    58    mp_limb_t d;
    59    ASSERT (n > 0);
    60  
    61    d = ap[0] ^ b;
    62  
    63    while (--n > 0)
    64      d |= ap[n];
    65  
    66    return d == 0;
    67  }
    68  
    69  mp_size_t
    70  mpn_sec_invert_itch (mp_size_t n)
    71  {
    72    return 4*n;
    73  }
    74  
    75  /* Compute V <-- A^{-1} (mod M), in data-independent time. M must be
    76     odd. Returns 1 on success, and 0 on failure (i.e., if gcd (A, m) !=
    77     1). Inputs and outputs of size n, and no overlap allowed. The {ap,
    78     n} area is destroyed. For arbitrary inputs, bit_size should be
    79     2*n*GMP_NUMB_BITS, but if A or M are known to be smaller, e.g., if
    80     M = 2^521 - 1 and A < M, bit_size can be any bound on the sum of
    81     the bit sizes of A and M. */
    82  int
    83  mpn_sec_invert (mp_ptr vp, mp_ptr ap, mp_srcptr mp,
    84  		mp_size_t n, mp_bitcnt_t bit_size,
    85  		mp_ptr scratch)
    86  {
    87    ASSERT (n > 0);
    88    ASSERT (bit_size > 0);
    89    ASSERT (mp[0] & 1);
    90    ASSERT (! MPN_OVERLAP_P (ap, n, vp, n));
    91  #define bp (scratch + n)
    92  #define up (scratch + 2*n)
    93  #define m1hp (scratch + 3*n)
    94  
    95    /* Maintain
    96  
    97         a = u * orig_a (mod m)
    98         b = v * orig_a (mod m)
    99  
   100       and b odd at all times. Initially,
   101  
   102         a = a_orig, u = 1
   103         b = m,      v = 0
   104       */
   105  
   106  
   107    up[0] = 1;
   108    mpn_zero (up+1, n - 1);
   109    mpn_copyi (bp, mp, n);
   110    mpn_zero (vp, n);
   111  
   112    ASSERT_CARRY (mpn_rshift (m1hp, mp, n, 1));
   113    ASSERT_NOCARRY (mpn_sec_add_1 (m1hp, m1hp, n, 1, scratch));
   114  
   115    while (bit_size-- > 0)
   116      {
   117        mp_limb_t odd, swap, cy;
   118  
   119        /* Always maintain b odd. The logic of the iteration is as
   120  	 follows. For a, b:
   121  
   122  	   odd = a & 1
   123  	   a -= odd * b
   124  	   if (underflow from a-b)
   125  	     {
   126  	       b += a, assigns old a
   127  	       a = B^n-a
   128  	     }
   129  
   130  	   a /= 2
   131  
   132  	 For u, v:
   133  
   134  	   if (underflow from a - b)
   135  	     swap u, v
   136  	   u -= odd * v
   137  	   if (underflow from u - v)
   138  	     u += m
   139  
   140  	   u /= 2
   141  	   if (a one bit was shifted out)
   142  	     u += (m+1)/2
   143  
   144  	 As long as a > 0, the quantity
   145  
   146  	   (bitsize of a) + (bitsize of b)
   147  
   148  	 is reduced by at least one bit per iteration, hence after (bit_size of
   149  	 orig_a) + (bit_size of m) - 1 iterations we surely have a = 0. Then b
   150  	 = gcd(orig_a, m) and if b = 1 then also v = orig_a^{-1} (mod m).
   151        */
   152  
   153        ASSERT (bp[0] & 1);
   154        odd = ap[0] & 1;
   155  
   156        swap = mpn_cnd_sub_n (odd, ap, ap, bp, n);
   157        mpn_cnd_add_n (swap, bp, bp, ap, n);
   158        mpn_cnd_neg (swap, ap, ap, n, scratch);
   159  
   160        mpn_cnd_swap (swap, up, vp, n);
   161        cy = mpn_cnd_sub_n (odd, up, up, vp, n);
   162        cy -= mpn_cnd_add_n (cy, up, up, mp, n);
   163        ASSERT (cy == 0);
   164  
   165        cy = mpn_rshift (ap, ap, n, 1);
   166        ASSERT (cy == 0);
   167        cy = mpn_rshift (up, up, n, 1);
   168        cy = mpn_cnd_add_n (cy, up, up, m1hp, n);
   169        ASSERT (cy == 0);
   170      }
   171    /* Should be all zeros, but check only extreme limbs */
   172    ASSERT ( (ap[0] | ap[n-1]) == 0);
   173    /* Check if indeed gcd == 1. */
   174    return mpn_sec_eq_ui (bp, n, 1);
   175  #undef bp
   176  #undef up
   177  #undef m1hp
   178  }