1*2a598d0bSHerbert Xu // SPDX-License-Identifier: GPL-2.0-or-later
2*2a598d0bSHerbert Xu /* mpihelp-mul.c - MPI helper functions
3*2a598d0bSHerbert Xu * Copyright (C) 1994, 1996, 1998, 1999,
4*2a598d0bSHerbert Xu * 2000 Free Software Foundation, Inc.
5*2a598d0bSHerbert Xu *
6*2a598d0bSHerbert Xu * This file is part of GnuPG.
7*2a598d0bSHerbert Xu *
8*2a598d0bSHerbert Xu * Note: This code is heavily based on the GNU MP Library.
9*2a598d0bSHerbert Xu * Actually it's the same code with only minor changes in the
10*2a598d0bSHerbert Xu * way the data is stored; this is to support the abstraction
11*2a598d0bSHerbert Xu * of an optional secure memory allocation which may be used
12*2a598d0bSHerbert Xu * to avoid revealing of sensitive data due to paging etc.
13*2a598d0bSHerbert Xu * The GNU MP Library itself is published under the LGPL;
14*2a598d0bSHerbert Xu * however I decided to publish this code under the plain GPL.
15*2a598d0bSHerbert Xu */
16*2a598d0bSHerbert Xu
17*2a598d0bSHerbert Xu #include <linux/string.h>
18*2a598d0bSHerbert Xu #include "mpi-internal.h"
19*2a598d0bSHerbert Xu #include "longlong.h"
20*2a598d0bSHerbert Xu
21*2a598d0bSHerbert Xu #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
22*2a598d0bSHerbert Xu do { \
23*2a598d0bSHerbert Xu if ((size) < KARATSUBA_THRESHOLD) \
24*2a598d0bSHerbert Xu mul_n_basecase(prodp, up, vp, size); \
25*2a598d0bSHerbert Xu else \
26*2a598d0bSHerbert Xu mul_n(prodp, up, vp, size, tspace); \
27*2a598d0bSHerbert Xu } while (0);
28*2a598d0bSHerbert Xu
29*2a598d0bSHerbert Xu #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
30*2a598d0bSHerbert Xu do { \
31*2a598d0bSHerbert Xu if ((size) < KARATSUBA_THRESHOLD) \
32*2a598d0bSHerbert Xu mpih_sqr_n_basecase(prodp, up, size); \
33*2a598d0bSHerbert Xu else \
34*2a598d0bSHerbert Xu mpih_sqr_n(prodp, up, size, tspace); \
35*2a598d0bSHerbert Xu } while (0);
36*2a598d0bSHerbert Xu
37*2a598d0bSHerbert Xu /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
38*2a598d0bSHerbert Xu * both with SIZE limbs, and store the result at PRODP. 2 * SIZE limbs are
39*2a598d0bSHerbert Xu * always stored. Return the most significant limb.
40*2a598d0bSHerbert Xu *
41*2a598d0bSHerbert Xu * Argument constraints:
42*2a598d0bSHerbert Xu * 1. PRODP != UP and PRODP != VP, i.e. the destination
43*2a598d0bSHerbert Xu * must be distinct from the multiplier and the multiplicand.
44*2a598d0bSHerbert Xu *
45*2a598d0bSHerbert Xu *
46*2a598d0bSHerbert Xu * Handle simple cases with traditional multiplication.
47*2a598d0bSHerbert Xu *
48*2a598d0bSHerbert Xu * This is the most critical code of multiplication. All multiplies rely
49*2a598d0bSHerbert Xu * on this, both small and huge. Small ones arrive here immediately. Huge
50*2a598d0bSHerbert Xu * ones arrive here as this is the base case for Karatsuba's recursive
51*2a598d0bSHerbert Xu * algorithm below.
52*2a598d0bSHerbert Xu */
53*2a598d0bSHerbert Xu
54*2a598d0bSHerbert Xu static mpi_limb_t
mul_n_basecase(mpi_ptr_t prodp,mpi_ptr_t up,mpi_ptr_t vp,mpi_size_t size)55*2a598d0bSHerbert Xu mul_n_basecase(mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
56*2a598d0bSHerbert Xu {
57*2a598d0bSHerbert Xu mpi_size_t i;
58*2a598d0bSHerbert Xu mpi_limb_t cy;
59*2a598d0bSHerbert Xu mpi_limb_t v_limb;
60*2a598d0bSHerbert Xu
61*2a598d0bSHerbert Xu /* Multiply by the first limb in V separately, as the result can be
62*2a598d0bSHerbert Xu * stored (not added) to PROD. We also avoid a loop for zeroing. */
63*2a598d0bSHerbert Xu v_limb = vp[0];
64*2a598d0bSHerbert Xu if (v_limb <= 1) {
65*2a598d0bSHerbert Xu if (v_limb == 1)
66*2a598d0bSHerbert Xu MPN_COPY(prodp, up, size);
67*2a598d0bSHerbert Xu else
68*2a598d0bSHerbert Xu MPN_ZERO(prodp, size);
69*2a598d0bSHerbert Xu cy = 0;
70*2a598d0bSHerbert Xu } else
71*2a598d0bSHerbert Xu cy = mpihelp_mul_1(prodp, up, size, v_limb);
72*2a598d0bSHerbert Xu
73*2a598d0bSHerbert Xu prodp[size] = cy;
74*2a598d0bSHerbert Xu prodp++;
75*2a598d0bSHerbert Xu
76*2a598d0bSHerbert Xu /* For each iteration in the outer loop, multiply one limb from
77*2a598d0bSHerbert Xu * U with one limb from V, and add it to PROD. */
78*2a598d0bSHerbert Xu for (i = 1; i < size; i++) {
79*2a598d0bSHerbert Xu v_limb = vp[i];
80*2a598d0bSHerbert Xu if (v_limb <= 1) {
81*2a598d0bSHerbert Xu cy = 0;
82*2a598d0bSHerbert Xu if (v_limb == 1)
83*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp, prodp, up, size);
84*2a598d0bSHerbert Xu } else
85*2a598d0bSHerbert Xu cy = mpihelp_addmul_1(prodp, up, size, v_limb);
86*2a598d0bSHerbert Xu
87*2a598d0bSHerbert Xu prodp[size] = cy;
88*2a598d0bSHerbert Xu prodp++;
89*2a598d0bSHerbert Xu }
90*2a598d0bSHerbert Xu
91*2a598d0bSHerbert Xu return cy;
92*2a598d0bSHerbert Xu }
93*2a598d0bSHerbert Xu
94*2a598d0bSHerbert Xu static void
mul_n(mpi_ptr_t prodp,mpi_ptr_t up,mpi_ptr_t vp,mpi_size_t size,mpi_ptr_t tspace)95*2a598d0bSHerbert Xu mul_n(mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp,
96*2a598d0bSHerbert Xu mpi_size_t size, mpi_ptr_t tspace)
97*2a598d0bSHerbert Xu {
98*2a598d0bSHerbert Xu if (size & 1) {
99*2a598d0bSHerbert Xu /* The size is odd, and the code below doesn't handle that.
100*2a598d0bSHerbert Xu * Multiply the least significant (size - 1) limbs with a recursive
101*2a598d0bSHerbert Xu * call, and handle the most significant limb of S1 and S2
102*2a598d0bSHerbert Xu * separately.
103*2a598d0bSHerbert Xu * A slightly faster way to do this would be to make the Karatsuba
104*2a598d0bSHerbert Xu * code below behave as if the size were even, and let it check for
105*2a598d0bSHerbert Xu * odd size in the end. I.e., in essence move this code to the end.
106*2a598d0bSHerbert Xu * Doing so would save us a recursive call, and potentially make the
107*2a598d0bSHerbert Xu * stack grow a lot less.
108*2a598d0bSHerbert Xu */
109*2a598d0bSHerbert Xu mpi_size_t esize = size - 1; /* even size */
110*2a598d0bSHerbert Xu mpi_limb_t cy_limb;
111*2a598d0bSHerbert Xu
112*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(prodp, up, vp, esize, tspace);
113*2a598d0bSHerbert Xu cy_limb = mpihelp_addmul_1(prodp + esize, up, esize, vp[esize]);
114*2a598d0bSHerbert Xu prodp[esize + esize] = cy_limb;
115*2a598d0bSHerbert Xu cy_limb = mpihelp_addmul_1(prodp + esize, vp, size, up[esize]);
116*2a598d0bSHerbert Xu prodp[esize + size] = cy_limb;
117*2a598d0bSHerbert Xu } else {
118*2a598d0bSHerbert Xu /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
119*2a598d0bSHerbert Xu *
120*2a598d0bSHerbert Xu * Split U in two pieces, U1 and U0, such that
121*2a598d0bSHerbert Xu * U = U0 + U1*(B**n),
122*2a598d0bSHerbert Xu * and V in V1 and V0, such that
123*2a598d0bSHerbert Xu * V = V0 + V1*(B**n).
124*2a598d0bSHerbert Xu *
125*2a598d0bSHerbert Xu * UV is then computed recursively using the identity
126*2a598d0bSHerbert Xu *
127*2a598d0bSHerbert Xu * 2n n n n
128*2a598d0bSHerbert Xu * UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
129*2a598d0bSHerbert Xu * 1 1 1 0 0 1 0 0
130*2a598d0bSHerbert Xu *
131*2a598d0bSHerbert Xu * Where B = 2**BITS_PER_MP_LIMB.
132*2a598d0bSHerbert Xu */
133*2a598d0bSHerbert Xu mpi_size_t hsize = size >> 1;
134*2a598d0bSHerbert Xu mpi_limb_t cy;
135*2a598d0bSHerbert Xu int negflg;
136*2a598d0bSHerbert Xu
137*2a598d0bSHerbert Xu /* Product H. ________________ ________________
138*2a598d0bSHerbert Xu * |_____U1 x V1____||____U0 x V0_____|
139*2a598d0bSHerbert Xu * Put result in upper part of PROD and pass low part of TSPACE
140*2a598d0bSHerbert Xu * as new TSPACE.
141*2a598d0bSHerbert Xu */
142*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(prodp + size, up + hsize, vp + hsize, hsize,
143*2a598d0bSHerbert Xu tspace);
144*2a598d0bSHerbert Xu
145*2a598d0bSHerbert Xu /* Product M. ________________
146*2a598d0bSHerbert Xu * |_(U1-U0)(V0-V1)_|
147*2a598d0bSHerbert Xu */
148*2a598d0bSHerbert Xu if (mpihelp_cmp(up + hsize, up, hsize) >= 0) {
149*2a598d0bSHerbert Xu mpihelp_sub_n(prodp, up + hsize, up, hsize);
150*2a598d0bSHerbert Xu negflg = 0;
151*2a598d0bSHerbert Xu } else {
152*2a598d0bSHerbert Xu mpihelp_sub_n(prodp, up, up + hsize, hsize);
153*2a598d0bSHerbert Xu negflg = 1;
154*2a598d0bSHerbert Xu }
155*2a598d0bSHerbert Xu if (mpihelp_cmp(vp + hsize, vp, hsize) >= 0) {
156*2a598d0bSHerbert Xu mpihelp_sub_n(prodp + hsize, vp + hsize, vp, hsize);
157*2a598d0bSHerbert Xu negflg ^= 1;
158*2a598d0bSHerbert Xu } else {
159*2a598d0bSHerbert Xu mpihelp_sub_n(prodp + hsize, vp, vp + hsize, hsize);
160*2a598d0bSHerbert Xu /* No change of NEGFLG. */
161*2a598d0bSHerbert Xu }
162*2a598d0bSHerbert Xu /* Read temporary operands from low part of PROD.
163*2a598d0bSHerbert Xu * Put result in low part of TSPACE using upper part of TSPACE
164*2a598d0bSHerbert Xu * as new TSPACE.
165*2a598d0bSHerbert Xu */
166*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(tspace, prodp, prodp + hsize, hsize,
167*2a598d0bSHerbert Xu tspace + size);
168*2a598d0bSHerbert Xu
169*2a598d0bSHerbert Xu /* Add/copy product H. */
170*2a598d0bSHerbert Xu MPN_COPY(prodp + hsize, prodp + size, hsize);
171*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp + size, prodp + size,
172*2a598d0bSHerbert Xu prodp + size + hsize, hsize);
173*2a598d0bSHerbert Xu
174*2a598d0bSHerbert Xu /* Add product M (if NEGFLG M is a negative number) */
175*2a598d0bSHerbert Xu if (negflg)
176*2a598d0bSHerbert Xu cy -=
177*2a598d0bSHerbert Xu mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace,
178*2a598d0bSHerbert Xu size);
179*2a598d0bSHerbert Xu else
180*2a598d0bSHerbert Xu cy +=
181*2a598d0bSHerbert Xu mpihelp_add_n(prodp + hsize, prodp + hsize, tspace,
182*2a598d0bSHerbert Xu size);
183*2a598d0bSHerbert Xu
184*2a598d0bSHerbert Xu /* Product L. ________________ ________________
185*2a598d0bSHerbert Xu * |________________||____U0 x V0_____|
186*2a598d0bSHerbert Xu * Read temporary operands from low part of PROD.
187*2a598d0bSHerbert Xu * Put result in low part of TSPACE using upper part of TSPACE
188*2a598d0bSHerbert Xu * as new TSPACE.
189*2a598d0bSHerbert Xu */
190*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(tspace, up, vp, hsize, tspace + size);
191*2a598d0bSHerbert Xu
192*2a598d0bSHerbert Xu /* Add/copy Product L (twice) */
193*2a598d0bSHerbert Xu
194*2a598d0bSHerbert Xu cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
195*2a598d0bSHerbert Xu if (cy)
196*2a598d0bSHerbert Xu mpihelp_add_1(prodp + hsize + size,
197*2a598d0bSHerbert Xu prodp + hsize + size, hsize, cy);
198*2a598d0bSHerbert Xu
199*2a598d0bSHerbert Xu MPN_COPY(prodp, tspace, hsize);
200*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize,
201*2a598d0bSHerbert Xu hsize);
202*2a598d0bSHerbert Xu if (cy)
203*2a598d0bSHerbert Xu mpihelp_add_1(prodp + size, prodp + size, size, 1);
204*2a598d0bSHerbert Xu }
205*2a598d0bSHerbert Xu }
206*2a598d0bSHerbert Xu
mpih_sqr_n_basecase(mpi_ptr_t prodp,mpi_ptr_t up,mpi_size_t size)207*2a598d0bSHerbert Xu void mpih_sqr_n_basecase(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size)
208*2a598d0bSHerbert Xu {
209*2a598d0bSHerbert Xu mpi_size_t i;
210*2a598d0bSHerbert Xu mpi_limb_t cy_limb;
211*2a598d0bSHerbert Xu mpi_limb_t v_limb;
212*2a598d0bSHerbert Xu
213*2a598d0bSHerbert Xu /* Multiply by the first limb in V separately, as the result can be
214*2a598d0bSHerbert Xu * stored (not added) to PROD. We also avoid a loop for zeroing. */
215*2a598d0bSHerbert Xu v_limb = up[0];
216*2a598d0bSHerbert Xu if (v_limb <= 1) {
217*2a598d0bSHerbert Xu if (v_limb == 1)
218*2a598d0bSHerbert Xu MPN_COPY(prodp, up, size);
219*2a598d0bSHerbert Xu else
220*2a598d0bSHerbert Xu MPN_ZERO(prodp, size);
221*2a598d0bSHerbert Xu cy_limb = 0;
222*2a598d0bSHerbert Xu } else
223*2a598d0bSHerbert Xu cy_limb = mpihelp_mul_1(prodp, up, size, v_limb);
224*2a598d0bSHerbert Xu
225*2a598d0bSHerbert Xu prodp[size] = cy_limb;
226*2a598d0bSHerbert Xu prodp++;
227*2a598d0bSHerbert Xu
228*2a598d0bSHerbert Xu /* For each iteration in the outer loop, multiply one limb from
229*2a598d0bSHerbert Xu * U with one limb from V, and add it to PROD. */
230*2a598d0bSHerbert Xu for (i = 1; i < size; i++) {
231*2a598d0bSHerbert Xu v_limb = up[i];
232*2a598d0bSHerbert Xu if (v_limb <= 1) {
233*2a598d0bSHerbert Xu cy_limb = 0;
234*2a598d0bSHerbert Xu if (v_limb == 1)
235*2a598d0bSHerbert Xu cy_limb = mpihelp_add_n(prodp, prodp, up, size);
236*2a598d0bSHerbert Xu } else
237*2a598d0bSHerbert Xu cy_limb = mpihelp_addmul_1(prodp, up, size, v_limb);
238*2a598d0bSHerbert Xu
239*2a598d0bSHerbert Xu prodp[size] = cy_limb;
240*2a598d0bSHerbert Xu prodp++;
241*2a598d0bSHerbert Xu }
242*2a598d0bSHerbert Xu }
243*2a598d0bSHerbert Xu
244*2a598d0bSHerbert Xu void
mpih_sqr_n(mpi_ptr_t prodp,mpi_ptr_t up,mpi_size_t size,mpi_ptr_t tspace)245*2a598d0bSHerbert Xu mpih_sqr_n(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size, mpi_ptr_t tspace)
246*2a598d0bSHerbert Xu {
247*2a598d0bSHerbert Xu if (size & 1) {
248*2a598d0bSHerbert Xu /* The size is odd, and the code below doesn't handle that.
249*2a598d0bSHerbert Xu * Multiply the least significant (size - 1) limbs with a recursive
250*2a598d0bSHerbert Xu * call, and handle the most significant limb of S1 and S2
251*2a598d0bSHerbert Xu * separately.
252*2a598d0bSHerbert Xu * A slightly faster way to do this would be to make the Karatsuba
253*2a598d0bSHerbert Xu * code below behave as if the size were even, and let it check for
254*2a598d0bSHerbert Xu * odd size in the end. I.e., in essence move this code to the end.
255*2a598d0bSHerbert Xu * Doing so would save us a recursive call, and potentially make the
256*2a598d0bSHerbert Xu * stack grow a lot less.
257*2a598d0bSHerbert Xu */
258*2a598d0bSHerbert Xu mpi_size_t esize = size - 1; /* even size */
259*2a598d0bSHerbert Xu mpi_limb_t cy_limb;
260*2a598d0bSHerbert Xu
261*2a598d0bSHerbert Xu MPN_SQR_N_RECURSE(prodp, up, esize, tspace);
262*2a598d0bSHerbert Xu cy_limb = mpihelp_addmul_1(prodp + esize, up, esize, up[esize]);
263*2a598d0bSHerbert Xu prodp[esize + esize] = cy_limb;
264*2a598d0bSHerbert Xu cy_limb = mpihelp_addmul_1(prodp + esize, up, size, up[esize]);
265*2a598d0bSHerbert Xu
266*2a598d0bSHerbert Xu prodp[esize + size] = cy_limb;
267*2a598d0bSHerbert Xu } else {
268*2a598d0bSHerbert Xu mpi_size_t hsize = size >> 1;
269*2a598d0bSHerbert Xu mpi_limb_t cy;
270*2a598d0bSHerbert Xu
271*2a598d0bSHerbert Xu /* Product H. ________________ ________________
272*2a598d0bSHerbert Xu * |_____U1 x U1____||____U0 x U0_____|
273*2a598d0bSHerbert Xu * Put result in upper part of PROD and pass low part of TSPACE
274*2a598d0bSHerbert Xu * as new TSPACE.
275*2a598d0bSHerbert Xu */
276*2a598d0bSHerbert Xu MPN_SQR_N_RECURSE(prodp + size, up + hsize, hsize, tspace);
277*2a598d0bSHerbert Xu
278*2a598d0bSHerbert Xu /* Product M. ________________
279*2a598d0bSHerbert Xu * |_(U1-U0)(U0-U1)_|
280*2a598d0bSHerbert Xu */
281*2a598d0bSHerbert Xu if (mpihelp_cmp(up + hsize, up, hsize) >= 0)
282*2a598d0bSHerbert Xu mpihelp_sub_n(prodp, up + hsize, up, hsize);
283*2a598d0bSHerbert Xu else
284*2a598d0bSHerbert Xu mpihelp_sub_n(prodp, up, up + hsize, hsize);
285*2a598d0bSHerbert Xu
286*2a598d0bSHerbert Xu /* Read temporary operands from low part of PROD.
287*2a598d0bSHerbert Xu * Put result in low part of TSPACE using upper part of TSPACE
288*2a598d0bSHerbert Xu * as new TSPACE. */
289*2a598d0bSHerbert Xu MPN_SQR_N_RECURSE(tspace, prodp, hsize, tspace + size);
290*2a598d0bSHerbert Xu
291*2a598d0bSHerbert Xu /* Add/copy product H */
292*2a598d0bSHerbert Xu MPN_COPY(prodp + hsize, prodp + size, hsize);
293*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp + size, prodp + size,
294*2a598d0bSHerbert Xu prodp + size + hsize, hsize);
295*2a598d0bSHerbert Xu
296*2a598d0bSHerbert Xu /* Add product M (if NEGFLG M is a negative number). */
297*2a598d0bSHerbert Xu cy -= mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace, size);
298*2a598d0bSHerbert Xu
299*2a598d0bSHerbert Xu /* Product L. ________________ ________________
300*2a598d0bSHerbert Xu * |________________||____U0 x U0_____|
301*2a598d0bSHerbert Xu * Read temporary operands from low part of PROD.
302*2a598d0bSHerbert Xu * Put result in low part of TSPACE using upper part of TSPACE
303*2a598d0bSHerbert Xu * as new TSPACE. */
304*2a598d0bSHerbert Xu MPN_SQR_N_RECURSE(tspace, up, hsize, tspace + size);
305*2a598d0bSHerbert Xu
306*2a598d0bSHerbert Xu /* Add/copy Product L (twice). */
307*2a598d0bSHerbert Xu cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
308*2a598d0bSHerbert Xu if (cy)
309*2a598d0bSHerbert Xu mpihelp_add_1(prodp + hsize + size,
310*2a598d0bSHerbert Xu prodp + hsize + size, hsize, cy);
311*2a598d0bSHerbert Xu
312*2a598d0bSHerbert Xu MPN_COPY(prodp, tspace, hsize);
313*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize,
314*2a598d0bSHerbert Xu hsize);
315*2a598d0bSHerbert Xu if (cy)
316*2a598d0bSHerbert Xu mpihelp_add_1(prodp + size, prodp + size, size, 1);
317*2a598d0bSHerbert Xu }
318*2a598d0bSHerbert Xu }
319*2a598d0bSHerbert Xu
320*2a598d0bSHerbert Xu
mpihelp_mul_n(mpi_ptr_t prodp,mpi_ptr_t up,mpi_ptr_t vp,mpi_size_t size)321*2a598d0bSHerbert Xu void mpihelp_mul_n(mpi_ptr_t prodp,
322*2a598d0bSHerbert Xu mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
323*2a598d0bSHerbert Xu {
324*2a598d0bSHerbert Xu if (up == vp) {
325*2a598d0bSHerbert Xu if (size < KARATSUBA_THRESHOLD)
326*2a598d0bSHerbert Xu mpih_sqr_n_basecase(prodp, up, size);
327*2a598d0bSHerbert Xu else {
328*2a598d0bSHerbert Xu mpi_ptr_t tspace;
329*2a598d0bSHerbert Xu tspace = mpi_alloc_limb_space(2 * size);
330*2a598d0bSHerbert Xu mpih_sqr_n(prodp, up, size, tspace);
331*2a598d0bSHerbert Xu mpi_free_limb_space(tspace);
332*2a598d0bSHerbert Xu }
333*2a598d0bSHerbert Xu } else {
334*2a598d0bSHerbert Xu if (size < KARATSUBA_THRESHOLD)
335*2a598d0bSHerbert Xu mul_n_basecase(prodp, up, vp, size);
336*2a598d0bSHerbert Xu else {
337*2a598d0bSHerbert Xu mpi_ptr_t tspace;
338*2a598d0bSHerbert Xu tspace = mpi_alloc_limb_space(2 * size);
339*2a598d0bSHerbert Xu mul_n(prodp, up, vp, size, tspace);
340*2a598d0bSHerbert Xu mpi_free_limb_space(tspace);
341*2a598d0bSHerbert Xu }
342*2a598d0bSHerbert Xu }
343*2a598d0bSHerbert Xu }
344*2a598d0bSHerbert Xu
345*2a598d0bSHerbert Xu int
mpihelp_mul_karatsuba_case(mpi_ptr_t prodp,mpi_ptr_t up,mpi_size_t usize,mpi_ptr_t vp,mpi_size_t vsize,struct karatsuba_ctx * ctx)346*2a598d0bSHerbert Xu mpihelp_mul_karatsuba_case(mpi_ptr_t prodp,
347*2a598d0bSHerbert Xu mpi_ptr_t up, mpi_size_t usize,
348*2a598d0bSHerbert Xu mpi_ptr_t vp, mpi_size_t vsize,
349*2a598d0bSHerbert Xu struct karatsuba_ctx *ctx)
350*2a598d0bSHerbert Xu {
351*2a598d0bSHerbert Xu mpi_limb_t cy;
352*2a598d0bSHerbert Xu
353*2a598d0bSHerbert Xu if (!ctx->tspace || ctx->tspace_size < vsize) {
354*2a598d0bSHerbert Xu if (ctx->tspace)
355*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tspace);
356*2a598d0bSHerbert Xu ctx->tspace = mpi_alloc_limb_space(2 * vsize);
357*2a598d0bSHerbert Xu if (!ctx->tspace)
358*2a598d0bSHerbert Xu return -ENOMEM;
359*2a598d0bSHerbert Xu ctx->tspace_size = vsize;
360*2a598d0bSHerbert Xu }
361*2a598d0bSHerbert Xu
362*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(prodp, up, vp, vsize, ctx->tspace);
363*2a598d0bSHerbert Xu
364*2a598d0bSHerbert Xu prodp += vsize;
365*2a598d0bSHerbert Xu up += vsize;
366*2a598d0bSHerbert Xu usize -= vsize;
367*2a598d0bSHerbert Xu if (usize >= vsize) {
368*2a598d0bSHerbert Xu if (!ctx->tp || ctx->tp_size < vsize) {
369*2a598d0bSHerbert Xu if (ctx->tp)
370*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tp);
371*2a598d0bSHerbert Xu ctx->tp = mpi_alloc_limb_space(2 * vsize);
372*2a598d0bSHerbert Xu if (!ctx->tp) {
373*2a598d0bSHerbert Xu if (ctx->tspace)
374*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tspace);
375*2a598d0bSHerbert Xu ctx->tspace = NULL;
376*2a598d0bSHerbert Xu return -ENOMEM;
377*2a598d0bSHerbert Xu }
378*2a598d0bSHerbert Xu ctx->tp_size = vsize;
379*2a598d0bSHerbert Xu }
380*2a598d0bSHerbert Xu
381*2a598d0bSHerbert Xu do {
382*2a598d0bSHerbert Xu MPN_MUL_N_RECURSE(ctx->tp, up, vp, vsize, ctx->tspace);
383*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp, prodp, ctx->tp, vsize);
384*2a598d0bSHerbert Xu mpihelp_add_1(prodp + vsize, ctx->tp + vsize, vsize,
385*2a598d0bSHerbert Xu cy);
386*2a598d0bSHerbert Xu prodp += vsize;
387*2a598d0bSHerbert Xu up += vsize;
388*2a598d0bSHerbert Xu usize -= vsize;
389*2a598d0bSHerbert Xu } while (usize >= vsize);
390*2a598d0bSHerbert Xu }
391*2a598d0bSHerbert Xu
392*2a598d0bSHerbert Xu if (usize) {
393*2a598d0bSHerbert Xu if (usize < KARATSUBA_THRESHOLD) {
394*2a598d0bSHerbert Xu mpi_limb_t tmp;
395*2a598d0bSHerbert Xu if (mpihelp_mul(ctx->tspace, vp, vsize, up, usize, &tmp)
396*2a598d0bSHerbert Xu < 0)
397*2a598d0bSHerbert Xu return -ENOMEM;
398*2a598d0bSHerbert Xu } else {
399*2a598d0bSHerbert Xu if (!ctx->next) {
400*2a598d0bSHerbert Xu ctx->next = kzalloc(sizeof *ctx, GFP_KERNEL);
401*2a598d0bSHerbert Xu if (!ctx->next)
402*2a598d0bSHerbert Xu return -ENOMEM;
403*2a598d0bSHerbert Xu }
404*2a598d0bSHerbert Xu if (mpihelp_mul_karatsuba_case(ctx->tspace,
405*2a598d0bSHerbert Xu vp, vsize,
406*2a598d0bSHerbert Xu up, usize,
407*2a598d0bSHerbert Xu ctx->next) < 0)
408*2a598d0bSHerbert Xu return -ENOMEM;
409*2a598d0bSHerbert Xu }
410*2a598d0bSHerbert Xu
411*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp, prodp, ctx->tspace, vsize);
412*2a598d0bSHerbert Xu mpihelp_add_1(prodp + vsize, ctx->tspace + vsize, usize, cy);
413*2a598d0bSHerbert Xu }
414*2a598d0bSHerbert Xu
415*2a598d0bSHerbert Xu return 0;
416*2a598d0bSHerbert Xu }
417*2a598d0bSHerbert Xu
mpihelp_release_karatsuba_ctx(struct karatsuba_ctx * ctx)418*2a598d0bSHerbert Xu void mpihelp_release_karatsuba_ctx(struct karatsuba_ctx *ctx)
419*2a598d0bSHerbert Xu {
420*2a598d0bSHerbert Xu struct karatsuba_ctx *ctx2;
421*2a598d0bSHerbert Xu
422*2a598d0bSHerbert Xu if (ctx->tp)
423*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tp);
424*2a598d0bSHerbert Xu if (ctx->tspace)
425*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tspace);
426*2a598d0bSHerbert Xu for (ctx = ctx->next; ctx; ctx = ctx2) {
427*2a598d0bSHerbert Xu ctx2 = ctx->next;
428*2a598d0bSHerbert Xu if (ctx->tp)
429*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tp);
430*2a598d0bSHerbert Xu if (ctx->tspace)
431*2a598d0bSHerbert Xu mpi_free_limb_space(ctx->tspace);
432*2a598d0bSHerbert Xu kfree(ctx);
433*2a598d0bSHerbert Xu }
434*2a598d0bSHerbert Xu }
435*2a598d0bSHerbert Xu
436*2a598d0bSHerbert Xu /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
437*2a598d0bSHerbert Xu * and v (pointed to by VP, with VSIZE limbs), and store the result at
438*2a598d0bSHerbert Xu * PRODP. USIZE + VSIZE limbs are always stored, but if the input
439*2a598d0bSHerbert Xu * operands are normalized. Return the most significant limb of the
440*2a598d0bSHerbert Xu * result.
441*2a598d0bSHerbert Xu *
442*2a598d0bSHerbert Xu * NOTE: The space pointed to by PRODP is overwritten before finished
443*2a598d0bSHerbert Xu * with U and V, so overlap is an error.
444*2a598d0bSHerbert Xu *
445*2a598d0bSHerbert Xu * Argument constraints:
446*2a598d0bSHerbert Xu * 1. USIZE >= VSIZE.
447*2a598d0bSHerbert Xu * 2. PRODP != UP and PRODP != VP, i.e. the destination
448*2a598d0bSHerbert Xu * must be distinct from the multiplier and the multiplicand.
449*2a598d0bSHerbert Xu */
450*2a598d0bSHerbert Xu
451*2a598d0bSHerbert Xu int
mpihelp_mul(mpi_ptr_t prodp,mpi_ptr_t up,mpi_size_t usize,mpi_ptr_t vp,mpi_size_t vsize,mpi_limb_t * _result)452*2a598d0bSHerbert Xu mpihelp_mul(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t usize,
453*2a598d0bSHerbert Xu mpi_ptr_t vp, mpi_size_t vsize, mpi_limb_t *_result)
454*2a598d0bSHerbert Xu {
455*2a598d0bSHerbert Xu mpi_ptr_t prod_endp = prodp + usize + vsize - 1;
456*2a598d0bSHerbert Xu mpi_limb_t cy;
457*2a598d0bSHerbert Xu struct karatsuba_ctx ctx;
458*2a598d0bSHerbert Xu
459*2a598d0bSHerbert Xu if (vsize < KARATSUBA_THRESHOLD) {
460*2a598d0bSHerbert Xu mpi_size_t i;
461*2a598d0bSHerbert Xu mpi_limb_t v_limb;
462*2a598d0bSHerbert Xu
463*2a598d0bSHerbert Xu if (!vsize) {
464*2a598d0bSHerbert Xu *_result = 0;
465*2a598d0bSHerbert Xu return 0;
466*2a598d0bSHerbert Xu }
467*2a598d0bSHerbert Xu
468*2a598d0bSHerbert Xu /* Multiply by the first limb in V separately, as the result can be
469*2a598d0bSHerbert Xu * stored (not added) to PROD. We also avoid a loop for zeroing. */
470*2a598d0bSHerbert Xu v_limb = vp[0];
471*2a598d0bSHerbert Xu if (v_limb <= 1) {
472*2a598d0bSHerbert Xu if (v_limb == 1)
473*2a598d0bSHerbert Xu MPN_COPY(prodp, up, usize);
474*2a598d0bSHerbert Xu else
475*2a598d0bSHerbert Xu MPN_ZERO(prodp, usize);
476*2a598d0bSHerbert Xu cy = 0;
477*2a598d0bSHerbert Xu } else
478*2a598d0bSHerbert Xu cy = mpihelp_mul_1(prodp, up, usize, v_limb);
479*2a598d0bSHerbert Xu
480*2a598d0bSHerbert Xu prodp[usize] = cy;
481*2a598d0bSHerbert Xu prodp++;
482*2a598d0bSHerbert Xu
483*2a598d0bSHerbert Xu /* For each iteration in the outer loop, multiply one limb from
484*2a598d0bSHerbert Xu * U with one limb from V, and add it to PROD. */
485*2a598d0bSHerbert Xu for (i = 1; i < vsize; i++) {
486*2a598d0bSHerbert Xu v_limb = vp[i];
487*2a598d0bSHerbert Xu if (v_limb <= 1) {
488*2a598d0bSHerbert Xu cy = 0;
489*2a598d0bSHerbert Xu if (v_limb == 1)
490*2a598d0bSHerbert Xu cy = mpihelp_add_n(prodp, prodp, up,
491*2a598d0bSHerbert Xu usize);
492*2a598d0bSHerbert Xu } else
493*2a598d0bSHerbert Xu cy = mpihelp_addmul_1(prodp, up, usize, v_limb);
494*2a598d0bSHerbert Xu
495*2a598d0bSHerbert Xu prodp[usize] = cy;
496*2a598d0bSHerbert Xu prodp++;
497*2a598d0bSHerbert Xu }
498*2a598d0bSHerbert Xu
499*2a598d0bSHerbert Xu *_result = cy;
500*2a598d0bSHerbert Xu return 0;
501*2a598d0bSHerbert Xu }
502*2a598d0bSHerbert Xu
503*2a598d0bSHerbert Xu memset(&ctx, 0, sizeof ctx);
504*2a598d0bSHerbert Xu if (mpihelp_mul_karatsuba_case(prodp, up, usize, vp, vsize, &ctx) < 0)
505*2a598d0bSHerbert Xu return -ENOMEM;
506*2a598d0bSHerbert Xu mpihelp_release_karatsuba_ctx(&ctx);
507*2a598d0bSHerbert Xu *_result = *prod_endp;
508*2a598d0bSHerbert Xu return 0;
509*2a598d0bSHerbert Xu }
510