userapps/opensource/sshd/libtommath/bn_mp_toom_mul.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_TOOM_MUL_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   4  *
   5  * LibTomMath is a library that provides multiple-precision
   6  * integer arithmetic as well as number theoretic functionality.
   7  *
   8  * The library was designed directly after the MPI library by
   9  * Michael Fromberger but has been written from scratch with
  10  * additional optimizations in place.
  11  *
  12  * The library is free for all purposes without any express
  13  * guarantee it works.
  14  *
  15  * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
  16  */
  17
  18 /* multiplication using the Toom-Cook 3-way algorithm
  19  *
  20  * Much more complicated than Karatsuba but has a lower
  21  * asymptotic running time of O(N**1.464).  This algorithm is
  22  * only particularly useful on VERY large inputs
  23  * (we're talking 1000s of digits here...).
  24 */
  25 int mp_toom_mul(mp_int *a, mp_int *b, mp_int *c)
  26 {
  27     mp_int w0, w1, w2, w3, w4, tmp1, tmp2, a0, a1, a2, b0, b1, b2;
  28     int res, B;
  29
  30     /* init temps */
  31     if ((res = mp_init_multi(&w0, &w1, &w2, &w3, &w4,
  32                              &a0, &a1, &a2, &b0, &b1,
  33                              &b2, &tmp1, &tmp2, NULL)) != MP_OKAY) {
  34        return res;
  35     }
  36
  37     /* B */
  38     B = MIN(a->used, b->used) / 3;
  39
  40     /* a = a2 * B**2 + a1 * B + a0 */
  41     if ((res = mp_mod_2d(a, DIGIT_BIT * B, &a0)) != MP_OKAY) {
  42        goto ERR;
  43     }
  44
  45     if ((res = mp_copy(a, &a1)) != MP_OKAY) {
  46        goto ERR;
  47     }
  48     mp_rshd(&a1, B);
  49     mp_mod_2d(&a1, DIGIT_BIT * B, &a1);
  50
  51     if ((res = mp_copy(a, &a2)) != MP_OKAY) {
  52        goto ERR;
  53     }
  54     mp_rshd(&a2, B*2);
  55
  56     /* b = b2 * B**2 + b1 * B + b0 */
  57     if ((res = mp_mod_2d(b, DIGIT_BIT * B, &b0)) != MP_OKAY) {
  58        goto ERR;
  59     }
  60
  61     if ((res = mp_copy(b, &b1)) != MP_OKAY) {
  62        goto ERR;
  63     }
  64     mp_rshd(&b1, B);
  65     mp_mod_2d(&b1, DIGIT_BIT * B, &b1);
  66
  67     if ((res = mp_copy(b, &b2)) != MP_OKAY) {
  68        goto ERR;
  69     }
  70     mp_rshd(&b2, B*2);
  71
  72     /* w0 = a0*b0 */
  73     if ((res = mp_mul(&a0, &b0, &w0)) != MP_OKAY) {
  74        goto ERR;
  75     }
  76
  77     /* w4 = a2 * b2 */
  78     if ((res = mp_mul(&a2, &b2, &w4)) != MP_OKAY) {
  79        goto ERR;
  80     }
  81
  82     /* w1 = (a2 + 2(a1 + 2a0))(b2 + 2(b1 + 2b0)) */
  83     if ((res = mp_mul_2(&a0, &tmp1)) != MP_OKAY) {
  84        goto ERR;
  85     }
  86     if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
  87        goto ERR;
  88     }
  89     if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
  90        goto ERR;
  91     }
  92     if ((res = mp_add(&tmp1, &a2, &tmp1)) != MP_OKAY) {
  93        goto ERR;
  94     }
  95
  96     if ((res = mp_mul_2(&b0, &tmp2)) != MP_OKAY) {
  97        goto ERR;
  98     }
  99     if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
 100        goto ERR;
 101     }
 102     if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
 103        goto ERR;
 104     }
 105     if ((res = mp_add(&tmp2, &b2, &tmp2)) != MP_OKAY) {
 106        goto ERR;
 107     }
 108
 109     if ((res = mp_mul(&tmp1, &tmp2, &w1)) != MP_OKAY) {
 110        goto ERR;
 111     }
 112
 113     /* w3 = (a0 + 2(a1 + 2a2))(b0 + 2(b1 + 2b2)) */
 114     if ((res = mp_mul_2(&a2, &tmp1)) != MP_OKAY) {
 115        goto ERR;
 116     }
 117     if ((res = mp_add(&tmp1, &a1, &tmp1)) != MP_OKAY) {
 118        goto ERR;
 119     }
 120     if ((res = mp_mul_2(&tmp1, &tmp1)) != MP_OKAY) {
 121        goto ERR;
 122     }
 123     if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
 124        goto ERR;
 125     }
 126
 127     if ((res = mp_mul_2(&b2, &tmp2)) != MP_OKAY) {
 128        goto ERR;
 129     }
 130     if ((res = mp_add(&tmp2, &b1, &tmp2)) != MP_OKAY) {
 131        goto ERR;
 132     }
 133     if ((res = mp_mul_2(&tmp2, &tmp2)) != MP_OKAY) {
 134        goto ERR;
 135     }
 136     if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
 137        goto ERR;
 138     }
 139
 140     if ((res = mp_mul(&tmp1, &tmp2, &w3)) != MP_OKAY) {
 141        goto ERR;
 142     }
 143
 144
 145     /* w2 = (a2 + a1 + a0)(b2 + b1 + b0) */
 146     if ((res = mp_add(&a2, &a1, &tmp1)) != MP_OKAY) {
 147        goto ERR;
 148     }
 149     if ((res = mp_add(&tmp1, &a0, &tmp1)) != MP_OKAY) {
 150        goto ERR;
 151     }
 152     if ((res = mp_add(&b2, &b1, &tmp2)) != MP_OKAY) {
 153        goto ERR;
 154     }
 155     if ((res = mp_add(&tmp2, &b0, &tmp2)) != MP_OKAY) {
 156        goto ERR;
 157     }
 158     if ((res = mp_mul(&tmp1, &tmp2, &w2)) != MP_OKAY) {
 159        goto ERR;
 160     }
 161
 162     /* now solve the matrix
 163
 164        0  0  0  0  1
 165        1  2  4  8  16
 166        1  1  1  1  1
 167        16 8  4  2  1
 168        1  0  0  0  0
 169
 170        using 12 subtractions, 4 shifts,
 171               2 small divisions and 1 small multiplication
 172      */
 173
 174      /* r1 - r4 */
 175      if ((res = mp_sub(&w1, &w4, &w1)) != MP_OKAY) {
 176         goto ERR;
 177      }
 178      /* r3 - r0 */
 179      if ((res = mp_sub(&w3, &w0, &w3)) != MP_OKAY) {
 180         goto ERR;
 181      }
 182      /* r1/2 */
 183      if ((res = mp_div_2(&w1, &w1)) != MP_OKAY) {
 184         goto ERR;
 185      }
 186      /* r3/2 */
 187      if ((res = mp_div_2(&w3, &w3)) != MP_OKAY) {
 188         goto ERR;
 189      }
 190      /* r2 - r0 - r4 */
 191      if ((res = mp_sub(&w2, &w0, &w2)) != MP_OKAY) {
 192         goto ERR;
 193      }
 194      if ((res = mp_sub(&w2, &w4, &w2)) != MP_OKAY) {
 195         goto ERR;
 196      }
 197      /* r1 - r2 */
 198      if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
 199         goto ERR;
 200      }
 201      /* r3 - r2 */
 202      if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
 203         goto ERR;
 204      }
 205      /* r1 - 8r0 */
 206      if ((res = mp_mul_2d(&w0, 3, &tmp1)) != MP_OKAY) {
 207         goto ERR;
 208      }
 209      if ((res = mp_sub(&w1, &tmp1, &w1)) != MP_OKAY) {
 210         goto ERR;
 211      }
 212      /* r3 - 8r4 */
 213      if ((res = mp_mul_2d(&w4, 3, &tmp1)) != MP_OKAY) {
 214         goto ERR;
 215      }
 216      if ((res = mp_sub(&w3, &tmp1, &w3)) != MP_OKAY) {
 217         goto ERR;
 218      }
 219      /* 3r2 - r1 - r3 */
 220      if ((res = mp_mul_d(&w2, 3, &w2)) != MP_OKAY) {
 221         goto ERR;
 222      }
 223      if ((res = mp_sub(&w2, &w1, &w2)) != MP_OKAY) {
 224         goto ERR;
 225      }
 226      if ((res = mp_sub(&w2, &w3, &w2)) != MP_OKAY) {
 227         goto ERR;
 228      }
 229      /* r1 - r2 */
 230      if ((res = mp_sub(&w1, &w2, &w1)) != MP_OKAY) {
 231         goto ERR;
 232      }
 233      /* r3 - r2 */
 234      if ((res = mp_sub(&w3, &w2, &w3)) != MP_OKAY) {
 235         goto ERR;
 236      }
 237      /* r1/3 */
 238      if ((res = mp_div_3(&w1, &w1, NULL)) != MP_OKAY) {
 239         goto ERR;
 240      }
 241      /* r3/3 */
 242      if ((res = mp_div_3(&w3, &w3, NULL)) != MP_OKAY) {
 243         goto ERR;
 244      }
 245
 246      /* at this point shift W[n] by B*n */
 247      if ((res = mp_lshd(&w1, 1*B)) != MP_OKAY) {
 248         goto ERR;
 249      }
 250      if ((res = mp_lshd(&w2, 2*B)) != MP_OKAY) {
 251         goto ERR;
 252      }
 253      if ((res = mp_lshd(&w3, 3*B)) != MP_OKAY) {
 254         goto ERR;
 255      }
 256      if ((res = mp_lshd(&w4, 4*B)) != MP_OKAY) {
 257         goto ERR;
 258      }
 259
 260      if ((res = mp_add(&w0, &w1, c)) != MP_OKAY) {
 261         goto ERR;
 262      }
 263      if ((res = mp_add(&w2, &w3, &tmp1)) != MP_OKAY) {
 264         goto ERR;
 265      }
 266      if ((res = mp_add(&w4, &tmp1, &tmp1)) != MP_OKAY) {
 267         goto ERR;
 268      }
 269      if ((res = mp_add(&tmp1, c, c)) != MP_OKAY) {
 270         goto ERR;
 271      }
 272
 273 ERR:
 274      mp_clear_multi(&w0, &w1, &w2, &w3, &w4,
 275                     &a0, &a1, &a2, &b0, &b1,
 276                     &b2, &tmp1, &tmp2, NULL);
 277      return res;
 278 }
 279
 280 #endif