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