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