Mercurial > dropbear

comparison bn_mp_karatsuba_sqr.c @ 2:86e0b50a9b58 libtommath-orig ltm-0.30-orig

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d29b64170cf0
comparison
equal deleted inserted replaced
-1:000000000000 2:86e0b50a9b58
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 }