|
1
|
1 /* LibTomMath, multiple-precision integer library -- Tom St Denis |
|
|
2 * |
|
|
3 * LibTomMath is a library that provides multiple-precision |
|
|
4 * integer arithmetic as well as number theoretic functionality. |
|
|
5 * |
|
|
6 * The library was designed directly after the MPI library by |
|
|
7 * Michael Fromberger but has been written from scratch with |
|
|
8 * additional optimizations in place. |
|
|
9 * |
|
|
10 * The library is free for all purposes without any express |
|
|
11 * guarantee it works. |
|
|
12 * |
|
|
13 * Tom St Denis, [email protected], http://math.libtomcrypt.org |
|
|
14 */ |
|
|
15 #include <tommath.h> |
|
|
16 |
|
|
17 /* c = |a| * |b| using Karatsuba Multiplication using |
|
|
18 * three half size multiplications |
|
|
19 * |
|
|
20 * Let B represent the radix [e.g. 2**DIGIT_BIT] and |
|
|
21 * let n represent half of the number of digits in |
|
|
22 * the min(a,b) |
|
|
23 * |
|
|
24 * a = a1 * B**n + a0 |
|
|
25 * b = b1 * B**n + b0 |
|
|
26 * |
|
|
27 * Then, a * b => |
|
|
28 a1b1 * B**2n + ((a1 - a0)(b1 - b0) + a0b0 + a1b1) * B + a0b0 |
|
|
29 * |
|
|
30 * Note that a1b1 and a0b0 are used twice and only need to be |
|
|
31 * computed once. So in total three half size (half # of |
|
|
32 * digit) multiplications are performed, a0b0, a1b1 and |
|
|
33 * (a1-b1)(a0-b0) |
|
|
34 * |
|
|
35 * Note that a multiplication of half the digits requires |
|
|
36 * 1/4th the number of single precision multiplications so in |
|
|
37 * total after one call 25% of the single precision multiplications |
|
|
38 * are saved. Note also that the call to mp_mul can end up back |
|
|
39 * in this function if the a0, a1, b0, or b1 are above the threshold. |
|
|
40 * This is known as divide-and-conquer and leads to the famous |
|
|
41 * O(N**lg(3)) or O(N**1.584) work which is asymptopically lower than |
|
|
42 * the standard O(N**2) that the baseline/comba methods use. |
|
|
43 * Generally though the overhead of this method doesn't pay off |
|
|
44 * until a certain size (N ~ 80) is reached. |
|
|
45 */ |
|
|
46 int mp_karatsuba_mul (mp_int * a, mp_int * b, mp_int * c) |
|
|
47 { |
|
|
48 mp_int x0, x1, y0, y1, t1, x0y0, x1y1; |
|
|
49 int B, err; |
|
|
50 |
|
|
51 /* default the return code to an error */ |
|
|
52 err = MP_MEM; |
|
|
53 |
|
|
54 /* min # of digits */ |
|
|
55 B = MIN (a->used, b->used); |
|
|
56 |
|
|
57 /* now divide in two */ |
|
|
58 B = B >> 1; |
|
|
59 |
|
|
60 /* init copy all the temps */ |
|
|
61 if (mp_init_size (&x0, B) != MP_OKAY) |
|
|
62 goto ERR; |
|
|
63 if (mp_init_size (&x1, a->used - B) != MP_OKAY) |
|
|
64 goto X0; |
|
|
65 if (mp_init_size (&y0, B) != MP_OKAY) |
|
|
66 goto X1; |
|
|
67 if (mp_init_size (&y1, b->used - B) != MP_OKAY) |
|
|
68 goto Y0; |
|
|
69 |
|
|
70 /* init temps */ |
|
|
71 if (mp_init_size (&t1, B * 2) != MP_OKAY) |
|
|
72 goto Y1; |
|
|
73 if (mp_init_size (&x0y0, B * 2) != MP_OKAY) |
|
|
74 goto T1; |
|
|
75 if (mp_init_size (&x1y1, B * 2) != MP_OKAY) |
|
|
76 goto X0Y0; |
|
|
77 |
|
|
78 /* now shift the digits */ |
|
|
79 x0.sign = x1.sign = a->sign; |
|
|
80 y0.sign = y1.sign = b->sign; |
|
|
81 |
|
|
82 x0.used = y0.used = B; |
|
|
83 x1.used = a->used - B; |
|
|
84 y1.used = b->used - B; |
|
|
85 |
|
|
86 { |
|
|
87 register int x; |
|
|
88 register mp_digit *tmpa, *tmpb, *tmpx, *tmpy; |
|
|
89 |
|
|
90 /* we copy the digits directly instead of using higher level functions |
|
|
91 * since we also need to shift the digits |
|
|
92 */ |
|
|
93 tmpa = a->dp; |
|
|
94 tmpb = b->dp; |
|
|
95 |
|
|
96 tmpx = x0.dp; |
|
|
97 tmpy = y0.dp; |
|
|
98 for (x = 0; x < B; x++) { |
|
|
99 *tmpx++ = *tmpa++; |
|
|
100 *tmpy++ = *tmpb++; |
|
|
101 } |
|
|
102 |
|
|
103 tmpx = x1.dp; |
|
|
104 for (x = B; x < a->used; x++) { |
|
|
105 *tmpx++ = *tmpa++; |
|
|
106 } |
|
|
107 |
|
|
108 tmpy = y1.dp; |
|
|
109 for (x = B; x < b->used; x++) { |
|
|
110 *tmpy++ = *tmpb++; |
|
|
111 } |
|
|
112 } |
|
|
113 |
|
|
114 /* only need to clamp the lower words since by definition the |
|
|
115 * upper words x1/y1 must have a known number of digits |
|
|
116 */ |
|
|
117 mp_clamp (&x0); |
|
|
118 mp_clamp (&y0); |
|
|
119 |
|
|
120 /* now calc the products x0y0 and x1y1 */ |
|
|
121 /* after this x0 is no longer required, free temp [x0==t2]! */ |
|
|
122 if (mp_mul (&x0, &y0, &x0y0) != MP_OKAY) |
|
|
123 goto X1Y1; /* x0y0 = x0*y0 */ |
|
|
124 if (mp_mul (&x1, &y1, &x1y1) != MP_OKAY) |
|
|
125 goto X1Y1; /* x1y1 = x1*y1 */ |
|
|
126 |
|
|
127 /* now calc x1-x0 and y1-y0 */ |
|
|
128 if (mp_sub (&x1, &x0, &t1) != MP_OKAY) |
|
|
129 goto X1Y1; /* t1 = x1 - x0 */ |
|
|
130 if (mp_sub (&y1, &y0, &x0) != MP_OKAY) |
|
|
131 goto X1Y1; /* t2 = y1 - y0 */ |
|
|
132 if (mp_mul (&t1, &x0, &t1) != MP_OKAY) |
|
|
133 goto X1Y1; /* t1 = (x1 - x0) * (y1 - y0) */ |
|
|
134 |
|
|
135 /* add x0y0 */ |
|
|
136 if (mp_add (&x0y0, &x1y1, &x0) != MP_OKAY) |
|
|
137 goto X1Y1; /* t2 = x0y0 + x1y1 */ |
|
|
138 if (mp_sub (&x0, &t1, &t1) != MP_OKAY) |
|
|
139 goto X1Y1; /* t1 = x0y0 + x1y1 - (x1-x0)*(y1-y0) */ |
|
|
140 |
|
|
141 /* shift by B */ |
|
|
142 if (mp_lshd (&t1, B) != MP_OKAY) |
|
|
143 goto X1Y1; /* t1 = (x0y0 + x1y1 - (x1-x0)*(y1-y0))<<B */ |
|
|
144 if (mp_lshd (&x1y1, B * 2) != MP_OKAY) |
|
|
145 goto X1Y1; /* x1y1 = x1y1 << 2*B */ |
|
|
146 |
|
|
147 if (mp_add (&x0y0, &t1, &t1) != MP_OKAY) |
|
|
148 goto X1Y1; /* t1 = x0y0 + t1 */ |
|
|
149 if (mp_add (&t1, &x1y1, c) != MP_OKAY) |
|
|
150 goto X1Y1; /* t1 = x0y0 + t1 + x1y1 */ |
|
|
151 |
|
|
152 /* Algorithm succeeded set the return code to MP_OKAY */ |
|
|
153 err = MP_OKAY; |
|
|
154 |
|
|
155 X1Y1:mp_clear (&x1y1); |
|
|
156 X0Y0:mp_clear (&x0y0); |
|
|
157 T1:mp_clear (&t1); |
|
|
158 Y1:mp_clear (&y1); |
|
|
159 Y0:mp_clear (&y0); |
|
|
160 X1:mp_clear (&x1); |
|
|
161 X0:mp_clear (&x0); |
|
|
162 ERR: |
|
|
163 return err; |
|
|
164 } |
