/home/ntakagi/work/tcl8.5.1/libtommath/bn_fast_s_mp_sqr.c Source File

00001 #include <tommath.h>
00002 #ifdef BN_FAST_S_MP_SQR_C
00003 /* LibTomMath, multiple-precision integer library -- Tom St Denis
00004  *
00005  * LibTomMath is a library that provides multiple-precision
00006  * integer arithmetic as well as number theoretic functionality.
00007  *
00008  * The library was designed directly after the MPI library by
00009  * Michael Fromberger but has been written from scratch with
00010  * additional optimizations in place.
00011  *
00012  * The library is free for all purposes without any express
00013  * guarantee it works.
00014  *
00015  * Tom St Denis, tomstdenis@gmail.com, http://math.libtomcrypt.com
00016  */
00017 
00018 /* the jist of squaring...
00019  * you do like mult except the offset of the tmpx [one that 
00020  * starts closer to zero] can't equal the offset of tmpy.  
00021  * So basically you set up iy like before then you min it with
00022  * (ty-tx) so that it never happens.  You double all those 
00023  * you add in the inner loop
00024 
00025 After that loop you do the squares and add them in.
00026 */
00027 
00028 int fast_s_mp_sqr (mp_int * a, mp_int * b)
00029 {
00030   int       olduse, res, pa, ix, iz;
00031   mp_digit   W[MP_WARRAY], *tmpx;
00032   mp_word   W1;
00033 
00034   /* grow the destination as required */
00035   pa = a->used + a->used;
00036   if (b->alloc < pa) {
00037     if ((res = mp_grow (b, pa)) != MP_OKAY) {
00038       return res;
00039     }
00040   }
00041 
00042   /* number of output digits to produce */
00043   W1 = 0;
00044   for (ix = 0; ix < pa; ix++) { 
00045       int      tx, ty, iy;
00046       mp_word  _W;
00047       mp_digit *tmpy;
00048 
00049       /* clear counter */
00050       _W = 0;
00051 
00052       /* get offsets into the two bignums */
00053       ty = MIN(a->used-1, ix);
00054       tx = ix - ty;
00055 
00056       /* setup temp aliases */
00057       tmpx = a->dp + tx;
00058       tmpy = a->dp + ty;
00059 
00060       /* this is the number of times the loop will iterrate, essentially
00061          while (tx++ < a->used && ty-- >= 0) { ... }
00062        */
00063       iy = MIN(a->used-tx, ty+1);
00064 
00065       /* now for squaring tx can never equal ty 
00066        * we halve the distance since they approach at a rate of 2x
00067        * and we have to round because odd cases need to be executed
00068        */
00069       iy = MIN(iy, (ty-tx+1)>>1);
00070 
00071       /* execute loop */
00072       for (iz = 0; iz < iy; iz++) {
00073          _W += ((mp_word)*tmpx++)*((mp_word)*tmpy--);
00074       }
00075 
00076       /* double the inner product and add carry */
00077       _W = _W + _W + W1;
00078 
00079       /* even columns have the square term in them */
00080       if ((ix&1) == 0) {
00081          _W += ((mp_word)a->dp[ix>>1])*((mp_word)a->dp[ix>>1]);
00082       }
00083 
00084       /* store it */
00085       W[ix] = (mp_digit)(_W & MP_MASK);
00086 
00087       /* make next carry */
00088       W1 = _W >> ((mp_word)DIGIT_BIT);
00089   }
00090 
00091   /* setup dest */
00092   olduse  = b->used;
00093   b->used = a->used+a->used;
00094 
00095   {
00096     mp_digit *tmpb;
00097     tmpb = b->dp;
00098     for (ix = 0; ix < pa; ix++) {
00099       *tmpb++ = W[ix] & MP_MASK;
00100     }
00101 
00102     /* clear unused digits [that existed in the old copy of c] */
00103     for (; ix < olduse; ix++) {
00104       *tmpb++ = 0;
00105     }
00106   }
00107   mp_clamp (b);
00108   return MP_OKAY;
00109 }
00110 #endif
00111 
00112 /* $Source: /cvsroot/tcl/libtommath/bn_fast_s_mp_sqr.c,v $ */
00113 /* $Revision: 1.1.1.4 $ */
00114 /* $Date: 2006/12/01 00:08:11 $ */