Mercurial > dropbear

comparison libtommath/bn_s_mp_toom_sqr.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_TOOM_SQR_C
3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
4 /* SPDX-License-Identifier: Unlicense */
5
6 /* squaring using Toom-Cook 3-way algorithm */
7
8 /*
9 This file contains code from J. Arndt's book "Matters Computational"
10 and the accompanying FXT-library with permission of the author.
11 */
12
13 /* squaring using Toom-Cook 3-way algorithm */
14 /*
15 Setup and interpolation from algorithm SQR_3 in
16
17 Chung, Jaewook, and M. Anwar Hasan. "Asymmetric squaring formulae."
18 18th IEEE Symposium on Computer Arithmetic (ARITH'07). IEEE, 2007.
19
20 */
21 mp_err s_mp_toom_sqr(const mp_int *a, mp_int *b)
22 {
23 mp_int S0, a0, a1, a2;
24 mp_digit *tmpa, *tmpc;
25 int B, count;
26 mp_err err;
27
28
29 /* init temps */
30 if ((err = mp_init(&S0)) != MP_OKAY) {
31 return err;
32 }
33
34 /* B */
35 B = a->used / 3;
36
37 /** a = a2 * x^2 + a1 * x + a0; */
38 if ((err = mp_init_size(&a0, B)) != MP_OKAY) goto LBL_ERRa0;
39
40 a0.used = B;
41 if ((err = mp_init_size(&a1, B)) != MP_OKAY) goto LBL_ERRa1;
42 a1.used = B;
43 if ((err = mp_init_size(&a2, B + (a->used - (3 * B)))) != MP_OKAY) goto LBL_ERRa2;
44
45 tmpa = a->dp;
46 tmpc = a0.dp;
47 for (count = 0; count < B; count++) {
48 *tmpc++ = *tmpa++;
49 }
50 tmpc = a1.dp;
51 for (; count < (2 * B); count++) {
52 *tmpc++ = *tmpa++;
53 }
54 tmpc = a2.dp;
55 for (; count < a->used; count++) {
56 *tmpc++ = *tmpa++;
57 a2.used++;
58 }
59 mp_clamp(&a0);
60 mp_clamp(&a1);
61 mp_clamp(&a2);
62
63 /** S0 = a0^2; */
64 if ((err = mp_sqr(&a0, &S0)) != MP_OKAY) goto LBL_ERR;
65
66 /** \\S1 = (a2 + a1 + a0)^2 */
67 /** \\S2 = (a2 - a1 + a0)^2 */
68 /** \\S1 = a0 + a2; */
69 /** a0 = a0 + a2; */
70 if ((err = mp_add(&a0, &a2, &a0)) != MP_OKAY) goto LBL_ERR;
71 /** \\S2 = S1 - a1; */
72 /** b = a0 - a1; */
73 if ((err = mp_sub(&a0, &a1, b)) != MP_OKAY) goto LBL_ERR;
74 /** \\S1 = S1 + a1; */
75 /** a0 = a0 + a1; */
76 if ((err = mp_add(&a0, &a1, &a0)) != MP_OKAY) goto LBL_ERR;
77 /** \\S1 = S1^2; */
78 /** a0 = a0^2; */
79 if ((err = mp_sqr(&a0, &a0)) != MP_OKAY) goto LBL_ERR;
80 /** \\S2 = S2^2; */
81 /** b = b^2; */
82 if ((err = mp_sqr(b, b)) != MP_OKAY) goto LBL_ERR;
83
84 /** \\ S3 = 2 * a1 * a2 */
85 /** \\S3 = a1 * a2; */
86 /** a1 = a1 * a2; */
87 if ((err = mp_mul(&a1, &a2, &a1)) != MP_OKAY) goto LBL_ERR;
88 /** \\S3 = S3 << 1; */
89 /** a1 = a1 << 1; */
90 if ((err = mp_mul_2(&a1, &a1)) != MP_OKAY) goto LBL_ERR;
91
92 /** \\S4 = a2^2; */
93 /** a2 = a2^2; */
94 if ((err = mp_sqr(&a2, &a2)) != MP_OKAY) goto LBL_ERR;
95
96 /** \\ tmp = (S1 + S2)/2 */
97 /** \\tmp = S1 + S2; */
98 /** b = a0 + b; */
99 if ((err = mp_add(&a0, b, b)) != MP_OKAY) goto LBL_ERR;
100 /** \\tmp = tmp >> 1; */
101 /** b = b >> 1; */
102 if ((err = mp_div_2(b, b)) != MP_OKAY) goto LBL_ERR;
103
104 /** \\ S1 = S1 - tmp - S3 */
105 /** \\S1 = S1 - tmp; */
106 /** a0 = a0 - b; */
107 if ((err = mp_sub(&a0, b, &a0)) != MP_OKAY) goto LBL_ERR;
108 /** \\S1 = S1 - S3; */
109 /** a0 = a0 - a1; */
110 if ((err = mp_sub(&a0, &a1, &a0)) != MP_OKAY) goto LBL_ERR;
111
112 /** \\S2 = tmp - S4 -S0 */
113 /** \\S2 = tmp - S4; */
114 /** b = b - a2; */
115 if ((err = mp_sub(b, &a2, b)) != MP_OKAY) goto LBL_ERR;
116 /** \\S2 = S2 - S0; */
117 /** b = b - S0; */
118 if ((err = mp_sub(b, &S0, b)) != MP_OKAY) goto LBL_ERR;
119
120
121 /** \\P = S4*x^4 + S3*x^3 + S2*x^2 + S1*x + S0; */
122 /** P = a2*x^4 + a1*x^3 + b*x^2 + a0*x + S0; */
123
124 if ((err = mp_lshd(&a2, 4 * B)) != MP_OKAY) goto LBL_ERR;
125 if ((err = mp_lshd(&a1, 3 * B)) != MP_OKAY) goto LBL_ERR;
126 if ((err = mp_lshd(b, 2 * B)) != MP_OKAY) goto LBL_ERR;
127 if ((err = mp_lshd(&a0, 1 * B)) != MP_OKAY) goto LBL_ERR;
128 if ((err = mp_add(&a2, &a1, &a2)) != MP_OKAY) goto LBL_ERR;
129 if ((err = mp_add(&a2, b, b)) != MP_OKAY) goto LBL_ERR;
130 if ((err = mp_add(b, &a0, b)) != MP_OKAY) goto LBL_ERR;
131 if ((err = mp_add(b, &S0, b)) != MP_OKAY) goto LBL_ERR;
132 /** a^2 - P */
133
134
135 LBL_ERR:
136 mp_clear(&a2);
137 LBL_ERRa2:
138 mp_clear(&a1);
139 LBL_ERRa1:
140 mp_clear(&a0);
141 LBL_ERRa0:
142 mp_clear(&S0);
143
144 return err;
145 }
146
147 #endif