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