/src/dropbear/libtommath/bn_s_mp_sqr_fast.c

Source
#include "tommath_private.h"
#ifdef BN_S_MP_SQR_FAST_C
/* LibTomMath, multiple-precision integer library -- Tom St Denis */
/* SPDX-License-Identifier: Unlicense */

/* the jist of squaring...
 * you do like mult except the offset of the tmpx [one that
 * starts closer to zero] can't equal the offset of tmpy.
 * So basically you set up iy like before then you min it with
 * (ty-tx) so that it never happens.  You double all those
 * you add in the inner loop

After that loop you do the squares and add them in.
*/

mp_err s_mp_sqr_fast(const mp_int *a, mp_int *b)
{
   int       olduse, pa, ix, iz;
   mp_digit  W[MP_WARRAY], *tmpx;
   mp_word   W1;
   mp_err    err;

   /* grow the destination as required */
   pa = a->used + a->used;
   if (b->alloc < pa) {
      if ((err = mp_grow(b, pa)) != MP_OKAY) {
         return err;
      }
   }

   /* number of output digits to produce */
   W1 = 0;
   for (ix = 0; ix < pa; ix++) {
      int      tx, ty, iy;
      mp_word  _W;
      mp_digit *tmpy;

      /* clear counter */
      _W = 0;

      /* get offsets into the two bignums */
      ty = MP_MIN(a->used-1, ix);
      tx = ix - ty;

      /* setup temp aliases */
      tmpx = a->dp + tx;
      tmpy = a->dp + ty;

      /* this is the number of times the loop will iterrate, essentially
         while (tx++ < a->used && ty-- >= 0) { ... }
       */
      iy = MP_MIN(a->used-tx, ty+1);

      /* now for squaring tx can never equal ty
       * we halve the distance since they approach at a rate of 2x
       * and we have to round because odd cases need to be executed
       */
      iy = MP_MIN(iy, ((ty-tx)+1)>>1);

      /* execute loop */
      for (iz = 0; iz < iy; iz++) {
         _W += (mp_word)*tmpx++ * (mp_word)*tmpy--;
      }

      /* double the inner product and add carry */
      _W = _W + _W + W1;

      /* even columns have the square term in them */
      if (((unsigned)ix & 1u) == 0u) {
         _W += (mp_word)a->dp[ix>>1] * (mp_word)a->dp[ix>>1];
      }

      /* store it */
      W[ix] = (mp_digit)_W & MP_MASK;

      /* make next carry */
      W1 = _W >> (mp_word)MP_DIGIT_BIT;
   }

   /* setup dest */
   olduse  = b->used;
   b->used = a->used+a->used;

   {
      mp_digit *tmpb;
      tmpb = b->dp;
      for (ix = 0; ix < pa; ix++) {
         *tmpb++ = W[ix] & MP_MASK;
      }

      /* clear unused digits [that existed in the old copy of c] */
      MP_ZERO_DIGITS(tmpb, olduse - ix);
   }
   mp_clamp(b);
   return MP_OKAY;
}
#endif

Line	Count	Source
1		#include "tommath_private.h"
2		#ifdef BN_S_MP_SQR_FAST_C
3		/* LibTomMath, multiple-precision integer library -- Tom St Denis */
4		/* SPDX-License-Identifier: Unlicense */
5
6		/* the jist of squaring...
7		* you do like mult except the offset of the tmpx [one that
8		* starts closer to zero] can't equal the offset of tmpy.
9		* So basically you set up iy like before then you min it with
10		* (ty-tx) so that it never happens. You double all those
11		* you add in the inner loop
12
13		After that loop you do the squares and add them in.
14		*/
15
16		mp_err s_mp_sqr_fast(const mp_int a, mp_int b)
17	65.0M	{
18	65.0M	int olduse, pa, ix, iz;
19	65.0M	mp_digit W[MP_WARRAY], *tmpx;
20	65.0M	mp_word W1;
21	65.0M	mp_err err;
22
23		/* grow the destination as required */
24	65.0M	pa = a->used + a->used;
25	65.0M	if (b->alloc < pa) {
26	87.3k	if ((err = mp_grow(b, pa)) != MP_OKAY) {
27	0	return err;
28	0	}
29	87.3k	}
30
31		/* number of output digits to produce */
32	65.0M	W1 = 0;
33	3.76G	for (ix = 0; ix < pa; ix++) {
34	3.70G	int tx, ty, iy;
35	3.70G	mp_word _W;
36	3.70G	mp_digit *tmpy;
37
38		/* clear counter */
39	3.70G	_W = 0;
40
41		/* get offsets into the two bignums */
42	3.70G	ty = MP_MIN(a->used-1, ix);
43	3.70G	tx = ix - ty;
44
45		/* setup temp aliases */
46	3.70G	tmpx = a->dp + tx;
47	3.70G	tmpy = a->dp + ty;
48
49		/* this is the number of times the loop will iterrate, essentially
50		while (tx++ < a->used && ty-- >= 0) { ... }
51		*/
52	3.70G	iy = MP_MIN(a->used-tx, ty+1);
53
54		/* now for squaring tx can never equal ty
55		* we halve the distance since they approach at a rate of 2x
56		* and we have to round because odd cases need to be executed
57		*/
58	3.70G	iy = MP_MIN(iy, ((ty-tx)+1)>>1);
59
60		/* execute loop */
61	33.7G	for (iz = 0; iz < iy; iz++) {
62	30.0G	_W += (mp_word)tmpx++ (mp_word)*tmpy--;
63	30.0G	}
64
65		/* double the inner product and add carry */
66	3.70G	_W = _W + _W + W1;
67
68		/* even columns have the square term in them */
69	3.70G	if (((unsigned)ix & 1u) == 0u) {
70	1.85G	_W += (mp_word)a->dp[ix>>1] * (mp_word)a->dp[ix>>1];
71	1.85G	}
72
73		/* store it */
74	3.70G	W[ix] = (mp_digit)_W & MP_MASK;
75
76		/* make next carry */
77	3.70G	W1 = _W >> (mp_word)MP_DIGIT_BIT;
78	3.70G	}
79
80		/* setup dest */
81	65.0M	olduse = b->used;
82	65.0M	b->used = a->used+a->used;
83
84	65.0M	{
85	65.0M	mp_digit *tmpb;
86	65.0M	tmpb = b->dp;
87	3.76G	for (ix = 0; ix < pa; ix++) {
88	3.70G	*tmpb++ = W[ix] & MP_MASK;
89	3.70G	}
90
91		/* clear unused digits [that existed in the old copy of c] */
92	65.0M	MP_ZERO_DIGITS(tmpb, olduse - ix);
93	65.0M	}
94	65.0M	mp_clamp(b);
95	65.0M	return MP_OKAY;
96	65.0M	}
97		#endif

Coverage Report

Created: 2025-11-20 06:33