1 /*  2 *******************************************************************************  3 \file zz_mul.c  4 \brief Multiple-precision unsigned integers: multiplicative operations  5 \project bee2 [cryptographic library]  6 \author (C) Sergey Agievich [agievich@{bsu.by|gmail.com}]  7 \created 2012.04.22  8 \version 2019.06.26  9 \license This program is released under the GNU General Public License  10 version 3. See Copyright Notices in bee2/info.h.  11 *******************************************************************************  12 */  13 14 #include "bee2/core/mem.h"  15 #include "bee2/core/util.h"  16 #include "bee2/core/word.h"  17 #include "bee2/math/ww.h"  18 #include "bee2/math/zz.h"  19 #include "zz_lcl.h"  20 21 /*  22 *******************************************************************************  23 Умножение / возведение в квадрат  24 25 \todo Возведение в квадрат за один проход (?), сначала с квадратов (?).  26 \todo Умножение Карацубы.  27 *******************************************************************************  28 */  29 30 1 word zzMulW(word b[], const word a[], size_t n, register word w)  31 {  32 1  register word carry = 0;  33  register dword prod;  34  size_t i;  35 1  ASSERT(wwIsSameOrDisjoint(a, b, n));  36 1  for (i = 0; i < n; ++i)  37  {  38 1  _MUL(prod, w, a[i]);  39 1  prod += carry;  40 1  b[i] = (word)prod;  41 1  carry = (word)(prod >> B_PER_W);  42  }  43 1  prod = 0, w = 0;  44 1  return carry;  45 }  46 47 1 word zzAddMulW(word b[], const word a[], size_t n, register word w)  48 {  49 1  register word carry = 0;  50  register dword prod;  51  size_t i;  52 1  ASSERT(wwIsSameOrDisjoint(a, b, n));  53 1  for (i = 0; i < n; ++i)  54  {  55 1  _MUL(prod, w, a[i]);  56 1  prod += carry;  57 1  prod += b[i];  58 1  b[i] = (word)prod;  59 1  carry = (word)(prod >> B_PER_W);  60  }  61 1  prod = 0, w = 0;  62 1  return carry;  63 }  64 65 1 word zzSubMulW(word b[], const word a[], size_t n, register word w)  66 {  67 1  register word borrow = 0;  68  register dword prod;  69  size_t i;  70 1  ASSERT(wwIsSameOrDisjoint(a, b, n));  71 1  for (i = 0; i < n; ++i)  72  {  73 1  _MUL(prod, w, a[i]);  74 1  prod = (dword)0 - prod;  75 1  prod += b[i];  76 1  prod -= borrow;  77 1  b[i] = (word)prod;  78 1  borrow = WORD_0 - (word)(prod >> B_PER_W);  79  }  80 1  prod = 0, w = 0;  81 1  return borrow;  82 }  83 84 1 void zzMul(word c[], const word a[], size_t n, const word b[], size_t m,  85  void* stack)  86 {  87 1  register word carry = 0;  88  register dword prod;  89  size_t i, j;  90 1  ASSERT(wwIsDisjoint2(a, n, c, n + m));  91 1  ASSERT(wwIsDisjoint2(b, m, c, n + m));  92 1  wwSetZero(c, n + m);  93 1  for (i = 0; i < n; ++i)  94  {  95 1  for (j = 0; j < m; ++j)  96  {  97 1  _MUL(prod, a[i], b[j]);  98 1  prod += carry;  99 1  prod += c[i + j];  100 1  c[i + j] = (word)prod;  101 1  carry = (word)(prod >> B_PER_W);  102  }  103 1  c[i + j] = carry;  104 1  carry = 0;  105  }  106 1  prod = 0;  107 }  108 109 1 size_t zzMul_deep(size_t n, size_t m)  110 {  111 1  return 0;  112 }  113 114 1 void zzSqr(word b[], const word a[], size_t n, void* stack)  115 {  116 1  register word carry = 0;  117  register word carry1;  118  register dword prod;  119  size_t i, j;  120 1  ASSERT(wwIsDisjoint2(a, n, b, n + n));  121  // b <- \sum_{i < j} a_i a_j B^{i + j}  122 1  wwSetZero(b, n + n);  123 1  for (i = 0; i < n; ++i)  124  {  125 1  for (j = i + 1; j < n; ++j)  126  {  127 1  _MUL(prod, a[i], a[j]);  128 1  prod += carry;  129 1  prod += b[i + j];  130 1  b[i + j] = (word)prod;  131 1  carry = (word)(prod >> B_PER_W);  132  }  133 1  b[i + j] = carry;  134 1  carry = 0;  135  }  136  // b <- 2 b  137 1  for (i = 0; i < n + n; ++i)  138  {  139 1  carry1 = b[i] >> (B_PER_W - 1);  140 1  b[i] = (b[i] << 1) | carry;  141 1  carry = carry1;  142  }  143  // b <- b + \sum_i a_i^2 B^{i + i}  144 1  for (i = 0; i < n; ++i)  145  {  146 1  _MUL(prod, a[i], a[i]);  147 1  prod += carry;  148 1  prod += b[i + i];  149 1  b[i + i] = (word)prod;  150 1  prod >>= B_PER_W;  151 1  prod += b[i + i + 1];  152 1  b[i + i + 1] = (word)prod;  153 1  carry = (word)(prod >> B_PER_W);  154  }  155 1  prod = 0;  156 1  carry = carry1 = 0;  157 }  158 159 1 size_t zzSqr_deep(size_t n)  160 {  161 1  return 0;  162 }  163 164 /*  165 *******************************************************************************  166 Деление на машинное слово  167 168 В функции zzModW2() сначала определяется значение (b = B \mod mod):  169  r = \sum_i a[i] b^i \equiv \sum_i a[i] B^i = a \mod mod,  170 которое затем приводится \mod mod.  171 Используется следующий алгоритм:  172  r = (r1 r0) <- 0  173  for i = n - 1,..., 0:  174  r <- (r1 b + r0)b + a[i] (*)  175  while (r1 != 0)  176  r <- r1 b + (r0 % mod) (**)  177  return r0 % mod  178 После каждой итерации (*):  179  r <= (B - 1)(1 + b + b^2) <= (B - 1)(mod^2 - mod + 1)  180  <= (B - 1)(B + 1) < B^2.  181 По окончании первой итерации (**):  182  r <= (B - 1)(mod - 1) + (mod - 1) = B(mod - 1).  183 По окончании второй итерации (**):  184  r <= (mod - 1)(mod - 1) + (mod - 1) = mod(mod - 1) < B.  185 Таким образом, r \mod mod = r0 \mod mod.  186 *******************************************************************************  187 */  188 189 1 word zzDivW(word q[], const word a[], size_t n, register word w)  190 {  191 1  register word r = 0;  192  register dword divisor;  193 1  ASSERT(w > 0);  194 1  ASSERT(wwIsSameOrDisjoint(a, q, n));  195 1  while (n--)  196  {  197 1  divisor = r;  198 1  divisor <<= B_PER_W;  199 1  divisor |= a[n];  200 1  q[n] = (word)(divisor / w);  201 1  r = (word)(divisor % w);  202  }  203 1  divisor = 0, w = 0;  204 1  return r;  205 }  206 207 1 word zzModW(const word a[], size_t n, register word w)  208 {  209 1  register word r = 0;  210  register dword divisor;  211 1  ASSERT(w > 0);  212 1  ASSERT(wwIsValid(a, n));  213 1  while (n--)  214  {  215 1  divisor = r;  216 1  divisor <<= B_PER_W;  217 1  divisor |= a[n];  218 1  r = (word)(divisor % w);  219  }  220 1  divisor = 0, w = 0;  221 1  return r;  222 }  223 224 1 word zzModW2(const word a[], size_t n, register word w)  225 {  226 1  register word r0 = 0;  227 1  register dword r1 = 0;  228  register word b;  229  // pre  230 1  ASSERT(w > 0);  231 1  ASSERT(w <= WORD_BIT_HALF);  232 1  ASSERT(wwIsValid(a, n));  233  // b <- B \mod mod  234 1  b = (WORD_MAX - w + 1) % w;  235  // (r1 r0) <- \sum_i a[i] b^i  236 1  while (n--)  237  {  238 1  r1 *= b;  239 1  r1 += r0;  240 1  r1 *= b;  241 1  r1 += a[n];  242 1  r0 = (word)r1;  243 1  r1 >>= B_PER_W;  244  }  245  // нормализация  246 #ifdef SAFE_FAST  247  while (r1 != 0)  248  {  249  r1 *= b;  250  r1 += r0 % w;  251  r0 = (word)r1;  252  r1 >>= B_PER_W;  253  }  254  r0 %= w;  255 #else  256 1  r1 *= b;  257 1  r1 += r0 % w;  258 1  r0 = (word)r1;  259 1  r1 >>= B_PER_W;  260 1  r1 *= b;  261 1  r1 += r0 % w;  262 1  r0 = (word)r1 % w;  263 #endif  264  // очистка и возврат  265 1  r1 = 0, b = w = 0;  266 1  return r0;  267 }  268 269 /*  270 *******************************************************************************  271 Общее деление  272 273 \opt При делении слов определять остаток по частному: (/,*) вместо (/, %).  274 275 \todo Убрать ограничение n >= m в zzDiv().  276 277 \todo T. Jabelean. An Algorithm for exact division. J. of Symb. Computations,  278 15 (2): 169-180, 1993.  279 280 \todo: В zzMod() отказаться от divident.  281 282 В функциях zzDiv(), zzMod() делимое a = a[n - 1]...a[0] и делитель  283 b = b[m - 1]...b[0] предварительно нормализуются:  284  a = a[n]...a[0] <- a * 2^shift;  285  b = b[m - 1]...b[0] <- b * 2^shift.  286 Здесь shift --- минимальное натуральное т.ч. старший бит b[m - 1] * 2^shift  287 равняется 1.  288 289 Деление выполняется по алгоритму 14.20 из [Menezes A., van Oorschot P.,  290 Vanstone S. Handbook of Applied Cryptography]:  291  for i = n, n - 1, ...., m:  292  if a[i] == b[m - 1] (#)  293  q[i - m] <- B - 1  294  else  295  q[i - m] <- a[i]a[i - 1] div b[m - 1]  296  while (q[i - m] * b[m - 1]b[m - 2] > a[i]a[i - 1]a[i - 2]) (##)  297  q[i - m]--  298  a <- a - q[i - m] * b * B^{i - m}  299  if (a < 0)  300  a += b * B^{i - m}, q[i - m]--  301  return q = q[n - m]...q[0] --- частное и a --- остаток  302 303 В реализации вместо (#):  304  d <- a[i]a[i - 1] div b[m - 1]  305  if d >= B  306  d <- B - 1  307  q[i - m] <- d  308 309 \opt Если a[i] == b[m - 1] в (#), то цикл (##) можно не выполнять:  310  q[i - m] * b[m - 1]b[m - 2] <=  311  (B - 1) * (a[i] * B + (B - 1)) =  312  B^2 * a[i] + B^2 - 1 - a[i] * B < a[i]a[i - 1]a[i - 2]  313 314 \opt Если известен остаток d = a[i]a[i - 1] mod b[m - 1], то (##) можно  315  заменить на  316  while (q[i - m] * b[m - 2] > d * B + a[i - 2])  317  q[i - m]--, d += b[m - 1]  318 *******************************************************************************  319 */  320 321 1 void zzDiv(word q[], word r[], const word a[], size_t n, const word b[],  322  size_t m, void* stack)  323 {  324  register dword dividentHi;  325  register word borrow;  326  register size_t shift;  327  size_t i;  328  // переменные в stack  329  word* divident; /*< нормализованное делимое (n + 1 слово) */  330  word* divisor; /*< нормализованный делитель (m слов) */  331  word* mul; /*< вспомогательное произведение (3 слова) */  332  // pre  333 1  ASSERT(n >= m);  334 1  ASSERT(wwIsValid(a, n) && wwIsValid(b, m));  335 1  ASSERT(m > 0 && b[m - 1] > 0);  336 1  ASSERT(wwIsDisjoint2(q, n + 1 - m, r, m));  337 1  ASSERT(a == r || wwIsDisjoint2(a, n, r, m));  338  // a < b?  339 1  if (wwCmp2(a, n, b, m) < 0)  340  {  341  // q <- 0, r <- a  342 1  wwSetZero(q, n - m + 1);  343 1  wwCopy(r, a, m);  344 1  return;  345  }  346  // делим на одноразрядное число?  347 1  if (m == 1)  348  {  349 1  r[0] = zzDivW(q, a, n, b[0]);  350 1  return;  351  }  352  // резервируем переменные в stack  353 1  divident = (word*)stack;  354 1  divisor = divident + n + 1;  355 1  mul = divisor + m;  356 1  stack = mul + 3;  357  // divident <- a  358 1  wwCopy(divident, a, n);  359 1  divident[n] = 0;  360  // divisor <- b  361 1  wwCopy(divisor, b, m);  362  // нормализация  363 1  shift = wordCLZ(b[m - 1]);  364 1  wwShHi(divident, n + 1, shift);  365 1  wwShHi(divisor, m, shift);  366  // цикл по разрядам делимого  367 1  for (i = n; i >= m; --i)  368  {  369  // вычислить пробное частное  370 1  dividentHi = divident[i];  371 1  dividentHi <<= B_PER_W;  372 1  dividentHi |= divident[i - 1];  373 1  dividentHi /= divisor[m - 1];  374 1  if (dividentHi > WORD_MAX)  375 0  q[i - m] = WORD_MAX;  376  else  377 1  q[i - m] = (word)dividentHi;  378  // уточнить пробное частное  379 1  wwCopy(mul, divisor + m - 2, 2);  380 1  mul[2] = zzMulW(mul, mul, 2, q[i - m]);  381 1  while (wwCmp2(mul, 3, divident + i - 2, 3) > 0)  382  {  383 1  q[i - m]--;  384 1  mul[2] -= zzSub2(mul, divisor + m - 2, 2);  385  }  386  // учесть пробное частное  387 1  borrow = zzSubMulW(divident + i - m, divisor, m, q[i - m]);  388 1  divident[i] -= borrow;  389 1  if (divident[i] > (word)~borrow)  390  {  391  // окончательно подправить пробное частное  392 0  q[i - m]--;  393  // корректирующее сложение  394 0  divident[i] += zzAdd2(divident + i - m, divisor, m);  395  }  396  }  397  // денормализация  398 1  wwShLo(divident, n + 1, shift);  399  // сохранить остаток  400 1  wwCopy(r, divident, m);  401  // очистить регистровые переменные  402 1  shift = 0;  403 1  borrow = 0;  404 1  dividentHi = 0;  405 }  406 407 1 size_t zzDiv_deep(size_t n, size_t m)  408 {  409 1  return O_OF_W(n + m + 4);  410 }  411 412 1 void zzMod(word r[], const word a[], size_t n, const word b[], size_t m, void* stack)  413 {  414  register dword dividentHi;  415  register word temp;  416  register size_t shift;  417  size_t i;  418  // переменные в stack  419  word* divident; /*< нормализованное делимое (n + 1 слово) */  420  word* divisor; /*< нормализованный делитель (m слов) */  421  word* mul; /*< вспомогательное произведение (3 слова) */  422  // pre  423 1  ASSERT(wwIsValid(a, n) && wwIsValid(b, m));  424 1  ASSERT(m > 0 && b[m - 1] > 0);  425 1  ASSERT(a == r || wwIsDisjoint2(a, n, r, m));  426  // a < b?  427 1  if (wwCmp2(a, n, b, m) < 0)  428  {  429  // r <- a  430 1  if (n < m)  431 0  wwSetZero(r + n, m - n), m = n;  432 1  wwCopy(r, a, m);  433 1  return;  434  }  435  // делим на одноразрядное число?  436 1  if (m == 1)  437  {  438 1  r[0] = zzModW(a, n, b[0]);  439 1  return;  440  }  441  // резервируем переменные в stack  442 1  divident = (word*)stack;  443 1  divisor = divident + n + 1;  444 1  mul = divisor + m;  445 1  stack = mul + 3;  446  // divident <- a  447 1  wwCopy(divident, a, n);  448 1  divident[n] = 0;  449  // divisor <- b  450 1  wwCopy(divisor, b, m);  451  // нормализация  452 1  shift = wordCLZ(b[m - 1]);  453 1  wwShHi(divident, n + 1, shift);  454 1  wwShHi(divisor, m, shift);  455  // цикл по разрядам делимого  456 1  for (i = n; i >= m; --i)  457  {  458  // вычислить пробное частное  459 1  dividentHi = divident[i];  460 1  dividentHi <<= B_PER_W;  461 1  dividentHi |= divident[i - 1];  462 1  dividentHi /= divisor[m - 1];  463 1  if (dividentHi > WORD_MAX)  464 1  temp = WORD_MAX;  465  else  466 1  temp = (word)dividentHi;  467  // уточнить пробное частное  468 1  wwCopy(mul, divisor + m - 2, 2);  469 1  mul[2] = zzMulW(mul, mul, 2, temp);  470 1  while (wwCmp2(mul, 3, divident + i - 2, 3) > 0)  471  {  472 1  temp--;  473 1  mul[2] -= zzSub2(mul, divisor + m - 2, 2);  474  }  475  // учесть пробное частное  476 1  temp = zzSubMulW(divident + i - m, divisor, m, temp);  477 1  divident[i] -= temp;  478 1  if (divident[i] > (word)~temp)  479  // корректирующее сложение  480 1  divident[i] += zzAdd2(divident + i - m, divisor, m);  481  }  482  // денормализация  483 1  wwShLo(divident, n + 1, shift);  484  // сохранить остаток  485 1  wwCopy(r, divident, m);  486  // очистить регистровые переменные  487 1  shift = 0;  488 1  temp = 0;  489 1  dividentHi = 0;  490 }  491 492 1 size_t zzMod_deep(size_t n, size_t m)  493 {  494 1  return O_OF_W(n + m + 4);  495 } 

Read our documentation on viewing source code .