sec-modinv.c

   1 /* sec-modinv.c */
   2
   3 /* nettle, low-level cryptographics library
   4  *
   5  * Copyright (C) 2013 Niels Möller
   6  *
   7  * The nettle library is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU Lesser General Public License as published by
   9  * the Free Software Foundation; either version 2.1 of the License, or (at your
  10  * option) any later version.
  11  *
  12  * The nettle library is distributed in the hope that it will be useful, but
  13  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15  * License for more details.
  16  *
  17  * You should have received a copy of the GNU Lesser General Public License
  18  * along with the nettle library; see the file COPYING.LIB.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
  20  * MA 02111-1301, USA.
  21  */
  22
  23 /* Development of Nettle's ECC support was funded by the .SE Internet Fund. */
  24
  25 #if HAVE_CONFIG_H
  26 # include "config.h"
  27 #endif
  28
  29 #include <assert.h>
  30
  31 #include "ecc-internal.h"
  32
  33 static void
  34 cnd_neg (int cnd, mp_limb_t *rp, const mp_limb_t *ap, mp_size_t n)
  35 {
  36   mp_limb_t cy = (cnd != 0);
  37   mp_limb_t mask = -cy;
  38   mp_size_t i;
  39
  40   for (i = 0; i < n; i++)
  41     {
  42       mp_limb_t r = (ap[i] ^ mask) + cy;
  43       cy = r < cy;
  44       rp[i] = r;
  45     }
  46 }
  47
  48 static void
  49 cnd_swap (int cnd, mp_limb_t *ap, mp_limb_t *bp, mp_size_t n)
  50 {
  51   mp_limb_t mask = - (mp_limb_t) (cnd != 0);
  52   mp_size_t i;
  53   for (i = 0; i < n; i++)
  54     {
  55       mp_limb_t a, b, t;
  56       a = ap[i];
  57       b = bp[i];
  58       t = (a ^ b) & mask;
  59       ap[i] = a ^ t;
  60       bp[i] = b ^ t;
  61     }
  62 }
  63
  64 /* Compute a^{-1} mod m, with running time depending only on the size.
  65    Also needs (m+1)/2, and m must be odd. */
  66 void
  67 sec_modinv (mp_limb_t *vp, mp_limb_t *ap, mp_size_t n,
  68             const mp_limb_t *mp, const mp_limb_t *mp1h, mp_size_t bit_size,
  69             mp_limb_t *scratch)
  70 {
  71 #define bp scratch
  72 #define dp (scratch + n)
  73 #define up (scratch + 2*n)
  74
  75   /* Avoid the mp_bitcnt_t type for compatibility with older GMP
  76      versions. */
  77   unsigned i;
  78
  79   /* Maintain
  80
  81        a = u * orig_a (mod m)
  82        b = v * orig_a (mod m)
  83
  84      and b odd at all times. Initially,
  85
  86        a = a_orig, u = 1
  87        b = m,      v = 0
  88      */
  89
  90   assert (ap != vp);
  91
  92   up[0] = 1;
  93   mpn_zero (up+1, n - 1);
  94   mpn_copyi (bp, mp, n);
  95   mpn_zero (vp, n);
  96
  97   for (i = bit_size + GMP_NUMB_BITS * n; i-- > 0; )
  98     {
  99       mp_limb_t odd, swap, cy;
 100
 101       /* Always maintain b odd. The logic of the iteration is as
 102          follows. For a, b:
 103
 104            odd = a & 1
 105            a -= odd * b
 106            if (underflow from a-b)
 107              {
 108                b += a, assigns old a
 109                a = B^n-a
 110              }
 111
 112            a /= 2
 113
 114          For u, v:
 115
 116            if (underflow from a - b)
 117              swap u, v
 118            u -= odd * v
 119            if (underflow from u - v)
 120              u += m
 121
 122            u /= 2
 123            if (a one bit was shifted out)
 124              u += (m+1)/2
 125
 126          As long as a > 0, the quantity
 127
 128            (bitsize of a) + (bitsize of b)
 129
 130          is reduced by at least one bit per iteration, hence after
 131          (bit_size of orig_a) + (bit_size of m) - 1 iterations we
 132          surely have a = 0. Then b = gcd(orig_a, m) and if b = 1 then
 133          also v = orig_a^{-1} (mod m)
 134       */
 135
 136       assert (bp[0] & 1);
 137       odd = ap[0] & 1;
 138
 139       /* Which variant is fastest depends on the speed of the various
 140          cnd_* functions. Assembly implementation would help. */
 141 #if 1
 142       swap = cnd_sub_n (odd, ap, bp, n);
 143       cnd_add_n (swap, bp, ap, n);
 144       cnd_neg (swap, ap, ap, n);
 145 #else
 146       swap = odd & mpn_sub_n (dp, ap, bp, n);
 147       cnd_copy (swap, bp, ap, n);
 148       cnd_neg (swap, dp, dp, n);
 149       cnd_copy (odd, ap, dp, n);
 150 #endif
 151
 152 #if 1
 153       cnd_swap (swap, up, vp, n);
 154       cy = cnd_sub_n (odd, up, vp, n);
 155       cy -= cnd_add_n (cy, up, mp, n);
 156 #else
 157       cy = cnd_sub_n (odd, up, vp, n);
 158       cnd_add_n (swap, vp, up, n);
 159       cnd_neg (swap, up, up, n);
 160       cnd_add_n (cy ^ swap, up, mp, n);
 161 #endif
 162       cy = mpn_rshift (ap, ap, n, 1);
 163       assert (cy == 0);
 164       cy = mpn_rshift (up, up, n, 1);
 165       cy = cnd_add_n (cy, up, mp1h, n);
 166       assert (cy == 0);
 167     }
 168   assert ( (ap[0] | ap[n-1]) == 0);
 169 #undef bp
 170 #undef dp
 171 #undef up
 172 }