Mercurial > dropbear

comparison libtommath/bn_s_mp_sqr_fast.c @ 1692:1051e4eea25a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update LibTomMath to 1.2.0 (#84) * update C files * update other files * update headers * update makefiles * remove mp_set/get_double() * use ltm 1.2.0 API * update ltm_desc * use bundled tommath if system-tommath is too old * XMALLOC etc. were changed to MP_MALLOC etc.
author Steffen Jaeckel <s@jaeckel.eu>
date Tue, 26 May 2020 17:36:47 +0200
parents
children
comparison
equal deleted inserted replaced
1691:2d3745d58843 1692:1051e4eea25a
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 {
18 int olduse, pa, ix, iz;
19 mp_digit W[MP_WARRAY], *tmpx;
20 mp_word W1;
21 mp_err err;
22
23 /* grow the destination as required */
24 pa = a->used + a->used;
25 if (b->alloc < pa) {
26 if ((err = mp_grow(b, pa)) != MP_OKAY) {
27 return err;
28 }
29 }
30
31 /* number of output digits to produce */
32 W1 = 0;
33 for (ix = 0; ix < pa; ix++) {
34 int tx, ty, iy;
35 mp_word _W;
36 mp_digit *tmpy;
37
38 /* clear counter */
39 _W = 0;
40
41 /* get offsets into the two bignums */
42 ty = MP_MIN(a->used-1, ix);
43 tx = ix - ty;
44
45 /* setup temp aliases */
46 tmpx = a->dp + tx;
47 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 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 iy = MP_MIN(iy, ((ty-tx)+1)>>1);
59
60 /* execute loop */
61 for (iz = 0; iz < iy; iz++) {
62 _W += (mp_word)*tmpx++ * (mp_word)*tmpy--;
63 }
64
65 /* double the inner product and add carry */
66 _W = _W + _W + W1;
67
68 /* even columns have the square term in them */
69 if (((unsigned)ix & 1u) == 0u) {
70 _W += (mp_word)a->dp[ix>>1] * (mp_word)a->dp[ix>>1];
71 }
72
73 /* store it */
74 W[ix] = (mp_digit)_W & MP_MASK;
75
76 /* make next carry */
77 W1 = _W >> (mp_word)MP_DIGIT_BIT;
78 }
79
80 /* setup dest */
81 olduse = b->used;
82 b->used = a->used+a->used;
83
84 {
85 mp_digit *tmpb;
86 tmpb = b->dp;
87 for (ix = 0; ix < pa; ix++) {
88 *tmpb++ = W[ix] & MP_MASK;
89 }
90
91 /* clear unused digits [that existed in the old copy of c] */
92 MP_ZERO_DIGITS(tmpb, olduse - ix);
93 }
94 mp_clamp(b);
95 return MP_OKAY;
96 }
97 #endif