source/libs/gmp/gmp-src/mpn/generic/sec_invert.c

   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 }