1 /* $NetBSD: bn_mp_karatsuba_sqr.c,v 1.1.1.1 2011/04/13 18:14:54 elric 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 */ 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