|
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 /* Karatsuba squaring, computes b = a*a using three |
|
|
18 * half size squarings |
|
|
19 * |
|
|
20 * See comments of mp_karatsuba_mul for details. It |
|
|
21 * is essentially the same algorithm but merely |
|
|
22 * tuned to perform recursive squarings. |
|
|
23 */ |
|
|
24 int mp_karatsuba_sqr (mp_int * a, mp_int * b) |
|
|
25 { |
|
|
26 mp_int x0, x1, t1, t2, x0x0, x1x1; |
|
|
27 int B, err; |
|
|
28 |
|
|
29 err = MP_MEM; |
|
|
30 |
|
|
31 /* min # of digits */ |
|
|
32 B = a->used; |
|
|
33 |
|
|
34 /* now divide in two */ |
|
|
35 B = B >> 1; |
|
|
36 |
|
|
37 /* init copy all the temps */ |
|
|
38 if (mp_init_size (&x0, B) != MP_OKAY) |
|
|
39 goto ERR; |
|
|
40 if (mp_init_size (&x1, a->used - B) != MP_OKAY) |
|
|
41 goto X0; |
|
|
42 |
|
|
43 /* init temps */ |
|
|
44 if (mp_init_size (&t1, a->used * 2) != MP_OKAY) |
|
|
45 goto X1; |
|
|
46 if (mp_init_size (&t2, a->used * 2) != MP_OKAY) |
|
|
47 goto T1; |
|
|
48 if (mp_init_size (&x0x0, B * 2) != MP_OKAY) |
|
|
49 goto T2; |
|
|
50 if (mp_init_size (&x1x1, (a->used - B) * 2) != MP_OKAY) |
|
|
51 goto X0X0; |
|
|
52 |
|
|
53 { |
|
|
54 register int x; |
|
|
55 register mp_digit *dst, *src; |
|
|
56 |
|
|
57 src = a->dp; |
|
|
58 |
|
|
59 /* now shift the digits */ |
|
|
60 dst = x0.dp; |
|
|
61 for (x = 0; x < B; x++) { |
|
|
62 *dst++ = *src++; |
|
|
63 } |
|
|
64 |
|
|
65 dst = x1.dp; |
|
|
66 for (x = B; x < a->used; x++) { |
|
|
67 *dst++ = *src++; |
|
|
68 } |
|
|
69 } |
|
|
70 |
|
|
71 x0.used = B; |
|
|
72 x1.used = a->used - B; |
|
|
73 |
|
|
74 mp_clamp (&x0); |
|
|
75 |
|
|
76 /* now calc the products x0*x0 and x1*x1 */ |
|
|
77 if (mp_sqr (&x0, &x0x0) != MP_OKAY) |
|
|
78 goto X1X1; /* x0x0 = x0*x0 */ |
|
|
79 if (mp_sqr (&x1, &x1x1) != MP_OKAY) |
|
|
80 goto X1X1; /* x1x1 = x1*x1 */ |
|
|
81 |
|
|
82 /* now calc (x1-x0)**2 */ |
|
|
83 if (mp_sub (&x1, &x0, &t1) != MP_OKAY) |
|
|
84 goto X1X1; /* t1 = x1 - x0 */ |
|
|
85 if (mp_sqr (&t1, &t1) != MP_OKAY) |
|
|
86 goto X1X1; /* t1 = (x1 - x0) * (x1 - x0) */ |
|
|
87 |
|
|
88 /* add x0y0 */ |
|
|
89 if (s_mp_add (&x0x0, &x1x1, &t2) != MP_OKAY) |
|
|
90 goto X1X1; /* t2 = x0x0 + x1x1 */ |
|
|
91 if (mp_sub (&t2, &t1, &t1) != MP_OKAY) |
|
|
92 goto X1X1; /* t1 = x0x0 + x1x1 - (x1-x0)*(x1-x0) */ |
|
|
93 |
|
|
94 /* shift by B */ |
|
|
95 if (mp_lshd (&t1, B) != MP_OKAY) |
|
|
96 goto X1X1; /* t1 = (x0x0 + x1x1 - (x1-x0)*(x1-x0))<<B */ |
|
|
97 if (mp_lshd (&x1x1, B * 2) != MP_OKAY) |
|
|
98 goto X1X1; /* x1x1 = x1x1 << 2*B */ |
|
|
99 |
|
|
100 if (mp_add (&x0x0, &t1, &t1) != MP_OKAY) |
|
|
101 goto X1X1; /* t1 = x0x0 + t1 */ |
|
|
102 if (mp_add (&t1, &x1x1, b) != MP_OKAY) |
|
|
103 goto X1X1; /* t1 = x0x0 + t1 + x1x1 */ |
|
|
104 |
|
|
105 err = MP_OKAY; |
|
|
106 |
|
|
107 X1X1:mp_clear (&x1x1); |
|
|
108 X0X0:mp_clear (&x0x0); |
|
|
109 T2:mp_clear (&t2); |
|
|
110 T1:mp_clear (&t1); |
|
|
111 X1:mp_clear (&x1); |
|
|
112 X0:mp_clear (&x0); |
|
|
113 ERR: |
|
|
114 return err; |
|
|
115 } |
