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