Mercurial > dropbear
annotate bn_mp_toom_mul.c @ 190:d8254fc979e9 libtommath-orig LTM_0.35
Initial import of libtommath 0.35
| author | Matt Johnston <matt@ucc.asn.au> |
|---|---|
| date | Fri, 06 May 2005 08:59:30 +0000 |
| parents | d29b64170cf0 |
| children |
| rev | line source |
|---|---|
| 142 | 1 #include <tommath.h> |
| 2 #ifdef BN_MP_TOOM_MUL_C | |
| 2 | 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, [email protected], http://math.libtomcrypt.org | |
| 16 */ | |
| 17 | |
| 142 | 18 /* multiplication using the Toom-Cook 3-way algorithm |
| 19 * | |
|
190
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
20 * Much more complicated than Karatsuba but has a lower |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
21 * asymptotic running time of O(N**1.464). This algorithm is |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
22 * only particularly useful on VERY large inputs |
|
d8254fc979e9
Initial import of libtommath 0.35
Matt Johnston <matt@ucc.asn.au>
parents:
142
diff
changeset
|
23 * (we're talking 1000s of digits here...). |
| 142 | 24 */ |
| 2 | 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 | |
| 142 | 280 #endif |