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 (at) gmail.com, http://math.libtomcrypt.com 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 281 282 /* $Source: /cvs/libtom/libtommath/bn_mp_toom_mul.c,v $ */ 283 /* $Revision: 1.3 $ */ 284 /* $Date: 2006/03/31 14:18:44 $ */ 285
