release/src/router/dropbear/libtommath/bn_mp_karatsuba_sqr.c

   1 #include <tommath.h>
   2 #ifdef BN_MP_KARATSUBA_SQR_C
   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, tomstdenis@gmail.com, http://math.libtomcrypt.com
  16  */
  17
  18 /* Karatsuba squaring, computes b = a*a using three
  19  * half size squarings
  20  *
  21  * See comments of karatsuba_mul for details.  It
  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 (s_mp_add (&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 (s_mp_sub (&t1, &t2, &t1) != MP_OKAY)
  93     goto X1X1;           /* t1 = (x1+x0)**2 - (x0x0 + x1x1) */
  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 }
 117 #endif
 118
 119 /* $Source: /cvs/libtom/libtommath/bn_mp_karatsuba_sqr.c,v $ */
 120 /* $Revision: 1.5 $ */
 121 /* $Date: 2006/03/31 14:18:44 $ */