1 /* mpihelp-mul.c - MPI helper functions
2 * Copyright (C) 1994, 1996, 1998, 1999,
3 * 2000 Free Software Foundation, Inc.
5 * This file is part of GnuPG.
7 * GnuPG is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
12 * GnuPG is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
21 * Note: This code is heavily based on the GNU MP Library.
22 * Actually it's the same code with only minor changes in the
23 * way the data is stored; this is to support the abstraction
24 * of an optional secure memory allocation which may be used
25 * to avoid revealing of sensitive data due to paging etc.
26 * The GNU MP Library itself is published under the LGPL;
27 * however I decided to publish this code under the plain GPL.
30 #include <linux/string.h>
31 #include "mpi-internal.h"
36 #define MPN_MUL_N_RECURSE(prodp, up, vp, size, tspace) \
38 if( (size) < KARATSUBA_THRESHOLD ) \
39 mul_n_basecase (prodp, up, vp, size); \
41 mul_n (prodp, up, vp, size, tspace); \
44 #define MPN_SQR_N_RECURSE(prodp, up, size, tspace) \
46 if ((size) < KARATSUBA_THRESHOLD) \
47 mpih_sqr_n_basecase (prodp, up, size); \
49 mpih_sqr_n (prodp, up, size, tspace); \
55 /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
56 * both with SIZE limbs, and store the result at PRODP. 2 * SIZE limbs are
57 * always stored. Return the most significant limb.
59 * Argument constraints:
60 * 1. PRODP != UP and PRODP != VP, i.e. the destination
61 * must be distinct from the multiplier and the multiplicand.
64 * Handle simple cases with traditional multiplication.
66 * This is the most critical code of multiplication. All multiplies rely
67 * on this, both small and huge. Small ones arrive here immediately. Huge
68 * ones arrive here as this is the base case for Karatsuba's recursive
73 mul_n_basecase( mpi_ptr_t prodp, mpi_ptr_t up,
74 mpi_ptr_t vp, mpi_size_t size)
80 /* Multiply by the first limb in V separately, as the result can be
81 * stored (not added) to PROD. We also avoid a loop for zeroing. */
85 MPN_COPY( prodp, up, size );
87 MPN_ZERO( prodp, size );
91 cy = mpihelp_mul_1( prodp, up, size, v_limb );
96 /* For each iteration in the outer loop, multiply one limb from
97 * U with one limb from V, and add it to PROD. */
98 for( i = 1; i < size; i++ ) {
103 cy = mpihelp_add_n(prodp, prodp, up, size);
106 cy = mpihelp_addmul_1(prodp, up, size, v_limb);
117 mul_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp,
118 mpi_size_t size, mpi_ptr_t tspace )
121 /* The size is odd, and the code below doesn't handle that.
122 * Multiply the least significant (size - 1) limbs with a recursive
123 * call, and handle the most significant limb of S1 and S2
125 * A slightly faster way to do this would be to make the Karatsuba
126 * code below behave as if the size were even, and let it check for
127 * odd size in the end. I.e., in essence move this code to the end.
128 * Doing so would save us a recursive call, and potentially make the
129 * stack grow a lot less.
131 mpi_size_t esize = size - 1; /* even size */
134 MPN_MUL_N_RECURSE( prodp, up, vp, esize, tspace );
135 cy_limb = mpihelp_addmul_1( prodp + esize, up, esize, vp[esize] );
136 prodp[esize + esize] = cy_limb;
137 cy_limb = mpihelp_addmul_1( prodp + esize, vp, size, up[esize] );
138 prodp[esize + size] = cy_limb;
141 /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
143 * Split U in two pieces, U1 and U0, such that
144 * U = U0 + U1*(B**n),
145 * and V in V1 and V0, such that
146 * V = V0 + V1*(B**n).
148 * UV is then computed recursively using the identity
151 * UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
154 * Where B = 2**BITS_PER_MP_LIMB.
156 mpi_size_t hsize = size >> 1;
160 /* Product H. ________________ ________________
161 * |_____U1 x V1____||____U0 x V0_____|
162 * Put result in upper part of PROD and pass low part of TSPACE
165 MPN_MUL_N_RECURSE(prodp + size, up + hsize, vp + hsize, hsize, tspace);
167 /* Product M. ________________
170 if( mpihelp_cmp(up + hsize, up, hsize) >= 0 ) {
171 mpihelp_sub_n(prodp, up + hsize, up, hsize);
175 mpihelp_sub_n(prodp, up, up + hsize, hsize);
178 if( mpihelp_cmp(vp + hsize, vp, hsize) >= 0 ) {
179 mpihelp_sub_n(prodp + hsize, vp + hsize, vp, hsize);
183 mpihelp_sub_n(prodp + hsize, vp, vp + hsize, hsize);
184 /* No change of NEGFLG. */
186 /* Read temporary operands from low part of PROD.
187 * Put result in low part of TSPACE using upper part of TSPACE
190 MPN_MUL_N_RECURSE(tspace, prodp, prodp + hsize, hsize, tspace + size);
192 /* Add/copy product H. */
193 MPN_COPY (prodp + hsize, prodp + size, hsize);
194 cy = mpihelp_add_n( prodp + size, prodp + size,
195 prodp + size + hsize, hsize);
197 /* Add product M (if NEGFLG M is a negative number) */
199 cy -= mpihelp_sub_n(prodp + hsize, prodp + hsize, tspace, size);
201 cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
203 /* Product L. ________________ ________________
204 * |________________||____U0 x V0_____|
205 * Read temporary operands from low part of PROD.
206 * Put result in low part of TSPACE using upper part of TSPACE
209 MPN_MUL_N_RECURSE(tspace, up, vp, hsize, tspace + size);
211 /* Add/copy Product L (twice) */
213 cy += mpihelp_add_n(prodp + hsize, prodp + hsize, tspace, size);
215 mpihelp_add_1(prodp + hsize + size, prodp + hsize + size, hsize, cy);
217 MPN_COPY(prodp, tspace, hsize);
218 cy = mpihelp_add_n(prodp + hsize, prodp + hsize, tspace + hsize, hsize);
220 mpihelp_add_1(prodp + size, prodp + size, size, 1);
226 mpih_sqr_n_basecase( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size )
232 /* Multiply by the first limb in V separately, as the result can be
233 * stored (not added) to PROD. We also avoid a loop for zeroing. */
237 MPN_COPY( prodp, up, size );
239 MPN_ZERO(prodp, size);
243 cy_limb = mpihelp_mul_1( prodp, up, size, v_limb );
245 prodp[size] = cy_limb;
248 /* For each iteration in the outer loop, multiply one limb from
249 * U with one limb from V, and add it to PROD. */
250 for( i=1; i < size; i++) {
255 cy_limb = mpihelp_add_n(prodp, prodp, up, size);
258 cy_limb = mpihelp_addmul_1(prodp, up, size, v_limb);
260 prodp[size] = cy_limb;
267 mpih_sqr_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t size, mpi_ptr_t tspace)
270 /* The size is odd, and the code below doesn't handle that.
271 * Multiply the least significant (size - 1) limbs with a recursive
272 * call, and handle the most significant limb of S1 and S2
274 * A slightly faster way to do this would be to make the Karatsuba
275 * code below behave as if the size were even, and let it check for
276 * odd size in the end. I.e., in essence move this code to the end.
277 * Doing so would save us a recursive call, and potentially make the
278 * stack grow a lot less.
280 mpi_size_t esize = size - 1; /* even size */
283 MPN_SQR_N_RECURSE( prodp, up, esize, tspace );
284 cy_limb = mpihelp_addmul_1( prodp + esize, up, esize, up[esize] );
285 prodp[esize + esize] = cy_limb;
286 cy_limb = mpihelp_addmul_1( prodp + esize, up, size, up[esize] );
288 prodp[esize + size] = cy_limb;
291 mpi_size_t hsize = size >> 1;
294 /* Product H. ________________ ________________
295 * |_____U1 x U1____||____U0 x U0_____|
296 * Put result in upper part of PROD and pass low part of TSPACE
299 MPN_SQR_N_RECURSE(prodp + size, up + hsize, hsize, tspace);
301 /* Product M. ________________
304 if( mpihelp_cmp( up + hsize, up, hsize) >= 0 )
305 mpihelp_sub_n( prodp, up + hsize, up, hsize);
307 mpihelp_sub_n (prodp, up, up + hsize, hsize);
309 /* Read temporary operands from low part of PROD.
310 * Put result in low part of TSPACE using upper part of TSPACE
312 MPN_SQR_N_RECURSE(tspace, prodp, hsize, tspace + size);
314 /* Add/copy product H */
315 MPN_COPY(prodp + hsize, prodp + size, hsize);
316 cy = mpihelp_add_n(prodp + size, prodp + size,
317 prodp + size + hsize, hsize);
319 /* Add product M (if NEGFLG M is a negative number). */
320 cy -= mpihelp_sub_n (prodp + hsize, prodp + hsize, tspace, size);
322 /* Product L. ________________ ________________
323 * |________________||____U0 x U0_____|
324 * Read temporary operands from low part of PROD.
325 * Put result in low part of TSPACE using upper part of TSPACE
327 MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
329 /* Add/copy Product L (twice). */
330 cy += mpihelp_add_n (prodp + hsize, prodp + hsize, tspace, size);
332 mpihelp_add_1(prodp + hsize + size, prodp + hsize + size,
335 MPN_COPY(prodp, tspace, hsize);
336 cy = mpihelp_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
338 mpihelp_add_1 (prodp + size, prodp + size, size, 1);
343 /* This should be made into an inline function in gmp.h. */
345 mpihelp_mul_n( mpi_ptr_t prodp, mpi_ptr_t up, mpi_ptr_t vp, mpi_size_t size)
348 if( size < KARATSUBA_THRESHOLD )
349 mpih_sqr_n_basecase( prodp, up, size );
352 tspace = mpi_alloc_limb_space( 2 * size );
355 mpih_sqr_n( prodp, up, size, tspace );
356 mpi_free_limb_space( tspace );
360 if( size < KARATSUBA_THRESHOLD )
361 mul_n_basecase( prodp, up, vp, size );
364 tspace = mpi_alloc_limb_space( 2 * size );
367 mul_n (prodp, up, vp, size, tspace);
368 mpi_free_limb_space( tspace );
378 mpihelp_mul_karatsuba_case( mpi_ptr_t prodp,
379 mpi_ptr_t up, mpi_size_t usize,
380 mpi_ptr_t vp, mpi_size_t vsize,
381 struct karatsuba_ctx *ctx )
385 if( !ctx->tspace || ctx->tspace_size < vsize ) {
387 mpi_free_limb_space( ctx->tspace );
388 ctx->tspace = mpi_alloc_limb_space( 2 * vsize);
391 ctx->tspace_size = vsize;
394 MPN_MUL_N_RECURSE( prodp, up, vp, vsize, ctx->tspace );
399 if( usize >= vsize ) {
400 if( !ctx->tp || ctx->tp_size < vsize ) {
402 mpi_free_limb_space( ctx->tp );
403 ctx->tp = mpi_alloc_limb_space( 2 * vsize );
406 mpi_free_limb_space( ctx->tspace );
410 ctx->tp_size = vsize;
414 MPN_MUL_N_RECURSE( ctx->tp, up, vp, vsize, ctx->tspace );
415 cy = mpihelp_add_n( prodp, prodp, ctx->tp, vsize );
416 mpihelp_add_1( prodp + vsize, ctx->tp + vsize, vsize, cy );
420 } while( usize >= vsize );
424 if( usize < KARATSUBA_THRESHOLD ) {
426 if (mpihelp_mul( ctx->tspace, vp, vsize, up, usize, &tmp) < 0)
431 ctx->next = kmalloc( sizeof *ctx, GFP_KERNEL );
434 memset(ctx->next, 0, sizeof(ctx));
436 if (mpihelp_mul_karatsuba_case( ctx->tspace,
443 cy = mpihelp_add_n( prodp, prodp, ctx->tspace, vsize);
444 mpihelp_add_1( prodp + vsize, ctx->tspace + vsize, usize, cy );
452 mpihelp_release_karatsuba_ctx( struct karatsuba_ctx *ctx )
454 struct karatsuba_ctx *ctx2;
457 mpi_free_limb_space( ctx->tp );
459 mpi_free_limb_space( ctx->tspace );
460 for( ctx=ctx->next; ctx; ctx = ctx2 ) {
463 mpi_free_limb_space( ctx->tp );
465 mpi_free_limb_space( ctx->tspace );
470 /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
471 * and v (pointed to by VP, with VSIZE limbs), and store the result at
472 * PRODP. USIZE + VSIZE limbs are always stored, but if the input
473 * operands are normalized. Return the most significant limb of the
476 * NOTE: The space pointed to by PRODP is overwritten before finished
477 * with U and V, so overlap is an error.
479 * Argument constraints:
481 * 2. PRODP != UP and PRODP != VP, i.e. the destination
482 * must be distinct from the multiplier and the multiplicand.
486 mpihelp_mul( mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t usize,
487 mpi_ptr_t vp, mpi_size_t vsize,
490 mpi_ptr_t prod_endp = prodp + usize + vsize - 1;
492 struct karatsuba_ctx ctx;
494 if( vsize < KARATSUBA_THRESHOLD ) {
503 /* Multiply by the first limb in V separately, as the result can be
504 * stored (not added) to PROD. We also avoid a loop for zeroing. */
508 MPN_COPY( prodp, up, usize );
510 MPN_ZERO( prodp, usize );
514 cy = mpihelp_mul_1( prodp, up, usize, v_limb );
519 /* For each iteration in the outer loop, multiply one limb from
520 * U with one limb from V, and add it to PROD. */
521 for( i = 1; i < vsize; i++ ) {
526 cy = mpihelp_add_n(prodp, prodp, up, usize);
529 cy = mpihelp_addmul_1(prodp, up, usize, v_limb);
539 memset( &ctx, 0, sizeof ctx );
540 if (mpihelp_mul_karatsuba_case( prodp, up, usize, vp, vsize, &ctx ) < 0)
542 mpihelp_release_karatsuba_ctx( &ctx );
543 *_result = *prod_endp;