userapps/opensource/sshd/libtommath/bn_mp_toom_mul.c

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