s_mp_mul_high.c

   1 #include "tommath_private.h"
   2 #ifdef S_MP_MUL_HIGH_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* multiplies |a| * |b| and does not compute the lower digs digits
   7  * [meant to get the higher part of the product]
   8  */
   9 mp_err s_mp_mul_high(const mp_int *a, const mp_int *b, mp_int *c, int digs)
  10 {
  11    mp_int   t;
  12    int      pa, pb, ix;
  13    mp_err   err;
  14
  15    if (digs < 0) {
  16       return MP_VAL;
  17    }
  18
  19    /* can we use the fast multiplier? */
  20    if (MP_HAS(S_MP_MUL_HIGH_COMBA)
  21        && ((a->used + b->used + 1) < MP_WARRAY)
  22        && (MP_MIN(a->used, b->used) < MP_MAX_COMBA)) {
  23       return s_mp_mul_high_comba(a, b, c, digs);
  24    }
  25
  26    if ((err = mp_init_size(&t, a->used + b->used + 1)) != MP_OKAY) {
  27       return err;
  28    }
  29    t.used = a->used + b->used + 1;
  30
  31    pa = a->used;
  32    pb = b->used;
  33    for (ix = 0; ix < pa; ix++) {
  34       int iy;
  35       mp_digit u = 0;
  36
  37       for (iy = digs - ix; iy < pb; iy++) {
  38          /* calculate the double precision result */
  39          mp_word r = (mp_word)t.dp[ix + iy] +
  40                      ((mp_word)a->dp[ix] * (mp_word)b->dp[iy]) +
  41                      (mp_word)u;
  42
  43          /* get the lower part */
  44          t.dp[ix + iy] = (mp_digit)(r & (mp_word)MP_MASK);
  45
  46          /* carry the carry */
  47          u       = (mp_digit)(r >> (mp_word)MP_DIGIT_BIT);
  48       }
  49       t.dp[ix + pb] = u;
  50    }
  51    mp_clamp(&t);
  52    mp_exch(&t, c);
  53    mp_clear(&t);
  54    return MP_OKAY;
  55 }
  56 #endif