s_mp_mul_comba.c

   1 #include "tommath_private.h"
   2 #ifdef S_MP_MUL_COMBA_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* Fast (comba) multiplier
   7  *
   8  * This is the fast column-array [comba] multiplier.  It is
   9  * designed to compute the columns of the product first
  10  * then handle the carries afterwards.  This has the effect
  11  * of making the nested loops that compute the columns very
  12  * simple and schedulable on super-scalar processors.
  13  *
  14  * This has been modified to produce a variable number of
  15  * digits of output so if say only a half-product is required
  16  * you don't have to compute the upper half (a feature
  17  * required for fast Barrett reduction).
  18  *
  19  * Based on Algorithm 14.12 on pp.595 of HAC.
  20  *
  21  */
  22 mp_err s_mp_mul_comba(const mp_int *a, const mp_int *b, mp_int *c, int digs)
  23 {
  24    int      oldused, pa, ix;
  25    mp_err   err;
  26    mp_digit MP_ALLOC_WARRAY(W);
  27    mp_word  _W;
  28
  29    MP_CHECK_WARRAY(W);
  30
  31    if (digs < 0) {
  32       MP_FREE_WARRAY(W);
  33       return MP_VAL;
  34    }
  35
  36    /* grow the destination as required */
  37    if ((err = mp_grow(c, digs)) != MP_OKAY) {
  38       MP_FREE_WARRAY(W);
  39       return err;
  40    }
  41
  42    /* number of output digits to produce */
  43    pa = MP_MIN(digs, a->used + b->used);
  44
  45    /* clear the carry */
  46    _W = 0;
  47    for (ix = 0; ix < pa; ix++) {
  48       int tx, ty, iy, iz;
  49
  50       /* get offsets into the two bignums */
  51       ty = MP_MIN(b->used-1, ix);
  52       tx = ix - ty;
  53
  54       /* this is the number of times the loop will iterate, essentially
  55          while (tx++ < a->used && ty-- >= 0) { ... }
  56        */
  57       iy = MP_MIN(a->used-tx, ty+1);
  58
  59       /* execute loop */
  60       for (iz = 0; iz < iy; ++iz) {
  61          _W += (mp_word)a->dp[tx + iz] * (mp_word)b->dp[ty - iz];
  62       }
  63
  64       /* store term */
  65       W[ix] = (mp_digit)_W & MP_MASK;
  66
  67       /* make next carry */
  68       _W = _W >> (mp_word)MP_DIGIT_BIT;
  69    }
  70
  71    /* setup dest */
  72    oldused  = c->used;
  73    c->used = pa;
  74
  75    for (ix = 0; ix < pa; ix++) {
  76       /* now extract the previous digit [below the carry] */
  77       c->dp[ix] = W[ix];
  78    }
  79
  80    /* clear unused digits [that existed in the old copy of c] */
  81    s_mp_zero_digs(c->dp + c->used, oldused - c->used);
  82
  83    mp_clamp(c);
  84    MP_FREE_WARRAY(W);
  85    return MP_OKAY;
  86 }
  87 #endif