485 lines · c
1// SPDX-License-Identifier: GPL-2.0-or-later2/* mpihelp-mul.c - MPI helper functions3 * 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 the10 * way the data is stored; this is to support the abstraction11 * of an optional secure memory allocation which may be used12 * 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 are39 * always stored. Return the most significant limb.40 *41 * Argument constraints:42 * 1. PRODP != UP and PRODP != VP, i.e. the destination43 * 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 rely49 * on this, both small and huge. Small ones arrive here immediately. Huge50 * ones arrive here as this is the base case for Karatsuba's recursive51 * algorithm below.52 */53 54static mpi_limb_t55mul_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 be62 * 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 else68 MPN_ZERO(prodp, size);69 cy = 0;70 } else71 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 from77 * 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 } else85 cy = mpihelp_addmul_1(prodp, up, size, v_limb);86 87 prodp[size] = cy;88 prodp++;89 }90 91 return cy;92}93 94static void95mul_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 recursive101 * call, and handle the most significant limb of S1 and S2102 * separately.103 * A slightly faster way to do this would be to make the Karatsuba104 * code below behave as if the size were even, and let it check for105 * 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 the107 * 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 that121 * U = U0 + U1*(B**n),122 * and V in V1 and V0, such that123 * V = V0 + V1*(B**n).124 *125 * UV is then computed recursively using the identity126 *127 * 2n n n n128 * UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V129 * 1 1 1 0 0 1 0 0130 *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 TSPACE140 * 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 TSPACE164 * 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 else180 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 TSPACE188 * 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 207void 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 be214 * 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 else220 MPN_ZERO(prodp, size);221 cy_limb = 0;222 } else223 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 from229 * 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 } else237 cy_limb = mpihelp_addmul_1(prodp, up, size, v_limb);238 239 prodp[size] = cy_limb;240 prodp++;241 }242}243 244void245mpih_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 recursive250 * call, and handle the most significant limb of S1 and S2251 * separately.252 * A slightly faster way to do this would be to make the Karatsuba253 * code below behave as if the size were even, and let it check for254 * 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 the256 * 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 TSPACE274 * 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 else284 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 TSPACE288 * 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 TSPACE303 * 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 320int321mpihelp_mul_karatsuba_case(mpi_ptr_t prodp,322 mpi_ptr_t up, mpi_size_t usize,323 mpi_ptr_t vp, mpi_size_t vsize,324 struct karatsuba_ctx *ctx)325{326 mpi_limb_t cy;327 328 if (!ctx->tspace || ctx->tspace_size < vsize) {329 if (ctx->tspace)330 mpi_free_limb_space(ctx->tspace);331 ctx->tspace = mpi_alloc_limb_space(2 * vsize);332 if (!ctx->tspace)333 return -ENOMEM;334 ctx->tspace_size = vsize;335 }336 337 MPN_MUL_N_RECURSE(prodp, up, vp, vsize, ctx->tspace);338 339 prodp += vsize;340 up += vsize;341 usize -= vsize;342 if (usize >= vsize) {343 if (!ctx->tp || ctx->tp_size < vsize) {344 if (ctx->tp)345 mpi_free_limb_space(ctx->tp);346 ctx->tp = mpi_alloc_limb_space(2 * vsize);347 if (!ctx->tp) {348 if (ctx->tspace)349 mpi_free_limb_space(ctx->tspace);350 ctx->tspace = NULL;351 return -ENOMEM;352 }353 ctx->tp_size = vsize;354 }355 356 do {357 MPN_MUL_N_RECURSE(ctx->tp, up, vp, vsize, ctx->tspace);358 cy = mpihelp_add_n(prodp, prodp, ctx->tp, vsize);359 mpihelp_add_1(prodp + vsize, ctx->tp + vsize, vsize,360 cy);361 prodp += vsize;362 up += vsize;363 usize -= vsize;364 } while (usize >= vsize);365 }366 367 if (usize) {368 if (usize < KARATSUBA_THRESHOLD) {369 mpi_limb_t tmp;370 if (mpihelp_mul(ctx->tspace, vp, vsize, up, usize, &tmp)371 < 0)372 return -ENOMEM;373 } else {374 if (!ctx->next) {375 ctx->next = kzalloc(sizeof *ctx, GFP_KERNEL);376 if (!ctx->next)377 return -ENOMEM;378 }379 if (mpihelp_mul_karatsuba_case(ctx->tspace,380 vp, vsize,381 up, usize,382 ctx->next) < 0)383 return -ENOMEM;384 }385 386 cy = mpihelp_add_n(prodp, prodp, ctx->tspace, vsize);387 mpihelp_add_1(prodp + vsize, ctx->tspace + vsize, usize, cy);388 }389 390 return 0;391}392 393void mpihelp_release_karatsuba_ctx(struct karatsuba_ctx *ctx)394{395 struct karatsuba_ctx *ctx2;396 397 if (ctx->tp)398 mpi_free_limb_space(ctx->tp);399 if (ctx->tspace)400 mpi_free_limb_space(ctx->tspace);401 for (ctx = ctx->next; ctx; ctx = ctx2) {402 ctx2 = ctx->next;403 if (ctx->tp)404 mpi_free_limb_space(ctx->tp);405 if (ctx->tspace)406 mpi_free_limb_space(ctx->tspace);407 kfree(ctx);408 }409}410 411/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)412 * and v (pointed to by VP, with VSIZE limbs), and store the result at413 * PRODP. USIZE + VSIZE limbs are always stored, but if the input414 * operands are normalized. Return the most significant limb of the415 * result.416 *417 * NOTE: The space pointed to by PRODP is overwritten before finished418 * with U and V, so overlap is an error.419 *420 * Argument constraints:421 * 1. USIZE >= VSIZE.422 * 2. PRODP != UP and PRODP != VP, i.e. the destination423 * must be distinct from the multiplier and the multiplicand.424 */425 426int427mpihelp_mul(mpi_ptr_t prodp, mpi_ptr_t up, mpi_size_t usize,428 mpi_ptr_t vp, mpi_size_t vsize, mpi_limb_t *_result)429{430 mpi_ptr_t prod_endp = prodp + usize + vsize - 1;431 mpi_limb_t cy;432 struct karatsuba_ctx ctx;433 434 if (vsize < KARATSUBA_THRESHOLD) {435 mpi_size_t i;436 mpi_limb_t v_limb;437 438 if (!vsize) {439 *_result = 0;440 return 0;441 }442 443 /* Multiply by the first limb in V separately, as the result can be444 * stored (not added) to PROD. We also avoid a loop for zeroing. */445 v_limb = vp[0];446 if (v_limb <= 1) {447 if (v_limb == 1)448 MPN_COPY(prodp, up, usize);449 else450 MPN_ZERO(prodp, usize);451 cy = 0;452 } else453 cy = mpihelp_mul_1(prodp, up, usize, v_limb);454 455 prodp[usize] = cy;456 prodp++;457 458 /* For each iteration in the outer loop, multiply one limb from459 * U with one limb from V, and add it to PROD. */460 for (i = 1; i < vsize; i++) {461 v_limb = vp[i];462 if (v_limb <= 1) {463 cy = 0;464 if (v_limb == 1)465 cy = mpihelp_add_n(prodp, prodp, up,466 usize);467 } else468 cy = mpihelp_addmul_1(prodp, up, usize, v_limb);469 470 prodp[usize] = cy;471 prodp++;472 }473 474 *_result = cy;475 return 0;476 }477 478 memset(&ctx, 0, sizeof ctx);479 if (mpihelp_mul_karatsuba_case(prodp, up, usize, vp, vsize, &ctx) < 0)480 return -ENOMEM;481 mpihelp_release_karatsuba_ctx(&ctx);482 *_result = *prod_endp;483 return 0;484}485