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

   1 /* mpn_sec_pi1_div_qr, mpn_sec_pi1_div_r -- Compute Q = floor(U / V), U = U
   2    mod V.  Side-channel silent under the assumption that the used instructions
   3    are side-channel silent.
   4
   5    Contributed to the GNU project by Torbjörn Granlund.
   6
   7    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
   8    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
   9    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  10
  11 Copyright 2011-2013 Free Software Foundation, Inc.
  12
  13 This file is part of the GNU MP Library.
  14
  15 The GNU MP Library is free software; you can redistribute it and/or modify
  16 it under the terms of either:
  17
  18   * the GNU Lesser General Public License as published by the Free
  19     Software Foundation; either version 3 of the License, or (at your
  20     option) any later version.
  21
  22 or
  23
  24   * the GNU General Public License as published by the Free Software
  25     Foundation; either version 2 of the License, or (at your option) any
  26     later version.
  27
  28 or both in parallel, as here.
  29
  30 The GNU MP Library is distributed in the hope that it will be useful, but
  31 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  32 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  33 for more details.
  34
  35 You should have received copies of the GNU General Public License and the
  36 GNU Lesser General Public License along with the GNU MP Library.  If not,
  37 see https://www.gnu.org/licenses/.  */
  38
  39 #include "gmp.h"
  40 #include "gmp-impl.h"
  41 #include "longlong.h"
  42
  43 /* This side-channel silent division algorithm reduces the partial remainder by
  44    GMP_NUMB_BITS/2 bits at a time, compared to GMP_NUMB_BITS for the main
  45    division algorithm.  We actually do not insist on reducing by exactly
  46    GMP_NUMB_BITS/2, but may leave a partial remainder that is D*B^i to 3D*B^i
  47    too large (B is the limb base, D is the divisor, and i is the induction
  48    variable); the subsequent step will handle the extra partial remainder bits.
  49
  50    With that partial remainder reduction, each step generates a quotient "half
  51    limb".  The outer loop generates two quotient half limbs, an upper (q1h) and
  52    a lower (q0h) which are stored sparsely in separate limb arrays.  These
  53    arrays are added at the end; using separate arrays avoids data-dependent
  54    carry propagation which could else pose a side-channel leakage problem.
  55
  56    The quotient half limbs may be between -3 to 0 from the accurate value
  57    ("accurate" being the one which corresponds to a reduction to a principal
  58    partial remainder).  Too small quotient half limbs correspond to too large
  59    remainders, which we reduce later, as described above.
  60
  61    In order to keep quotients from getting too big, corresponding to a negative
  62    partial remainder, we use an inverse which is slightly smaller than usually.
  63 */
  64
  65 #if OPERATION_sec_pi1_div_qr
  66 /* Needs (dn + 1) + (nn - dn) + (nn - dn) = 2nn - dn + 1 limbs at tp. */
  67 #define FNAME mpn_sec_pi1_div_qr
  68 #define Q(q) q,
  69 #define RETTYPE mp_limb_t
  70 #endif
  71 #if OPERATION_sec_pi1_div_r
  72 /* Needs (dn + 1) limbs at tp.  */
  73 #define FNAME mpn_sec_pi1_div_r
  74 #define Q(q)
  75 #define RETTYPE void
  76 #endif
  77
  78 RETTYPE
  79 FNAME (Q(mp_ptr qp)
  80        mp_ptr np, mp_size_t nn,
  81        mp_srcptr dp, mp_size_t dn,
  82        mp_limb_t dinv,
  83        mp_ptr tp)
  84 {
  85   mp_limb_t nh, cy, q1h, q0h, dummy, cnd;
  86   mp_size_t i;
  87   mp_ptr hp;
  88 #if OPERATION_sec_pi1_div_qr
  89   mp_limb_t qh;
  90   mp_ptr qlp, qhp;
  91 #endif
  92
  93   ASSERT (dn >= 1);
  94   ASSERT (nn >= dn);
  95   ASSERT ((dp[dn - 1] & GMP_NUMB_HIGHBIT) != 0);
  96
  97   if (nn == dn)
  98     {
  99       cy = mpn_sub_n (np, np, dp, dn);
 100       mpn_cnd_add_n (cy, np, np, dp, dn);
 101 #if OPERATION_sec_pi1_div_qr
 102       return 1 - cy;
 103 #else
 104       return;
 105 #endif
 106     }
 107
 108   /* Create a divisor copy shifted half a limb.  */
 109   hp = tp;                                      /* (dn + 1) limbs */
 110   hp[dn] = mpn_lshift (hp, dp, dn, GMP_NUMB_BITS / 2);
 111
 112 #if OPERATION_sec_pi1_div_qr
 113   qlp = tp + (dn + 1);                          /* (nn - dn) limbs */
 114   qhp = tp + (nn + 1);                          /* (nn - dn) limbs */
 115 #endif
 116
 117   np += nn - dn;
 118   nh = 0;
 119
 120   for (i = nn - dn - 1; i >= 0; i--)
 121     {
 122       np--;
 123
 124       nh = (nh << GMP_NUMB_BITS/2) + (np[dn] >> GMP_NUMB_BITS/2);
 125       umul_ppmm (q1h, dummy, nh, dinv);
 126       q1h += nh;
 127 #if OPERATION_sec_pi1_div_qr
 128       qhp[i] = q1h;
 129 #endif
 130       mpn_submul_1 (np, hp, dn + 1, q1h);
 131
 132       nh = np[dn];
 133       umul_ppmm (q0h, dummy, nh, dinv);
 134       q0h += nh;
 135 #if OPERATION_sec_pi1_div_qr
 136       qlp[i] = q0h;
 137 #endif
 138       nh -= mpn_submul_1 (np, dp, dn, q0h);
 139     }
 140
 141   /* 1st adjustment depends on extra high remainder limb.  */
 142   cnd = nh != 0;                                /* FIXME: cmp-to-int */
 143 #if OPERATION_sec_pi1_div_qr
 144   qlp[0] += cnd;
 145 #endif
 146   nh -= mpn_cnd_sub_n (cnd, np, np, dp, dn);
 147
 148   /* 2nd adjustment depends on remainder/divisor comparison as well as whether
 149      extra remainder limb was nullified by previous subtract.  */
 150   cy = mpn_sub_n (np, np, dp, dn);
 151   cy = cy - nh;
 152 #if OPERATION_sec_pi1_div_qr
 153   qlp[0] += 1 - cy;
 154 #endif
 155   mpn_cnd_add_n (cy, np, np, dp, dn);
 156
 157   /* 3rd adjustment depends on remainder/divisor comparison.  */
 158   cy = mpn_sub_n (np, np, dp, dn);
 159 #if OPERATION_sec_pi1_div_qr
 160   qlp[0] += 1 - cy;
 161 #endif
 162   mpn_cnd_add_n (cy, np, np, dp, dn);
 163
 164 #if OPERATION_sec_pi1_div_qr
 165   /* Combine quotient halves into final quotient.  */
 166   qh = mpn_lshift (qhp, qhp, nn - dn, GMP_NUMB_BITS/2);
 167   qh += mpn_add_n (qp, qhp, qlp, nn - dn);
 168
 169   return qh;
 170 #else
 171   return;
 172 #endif
 173 }