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 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 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 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 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 int 321 mpihelp_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 393 void 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 at 413 * PRODP. USIZE + VSIZE limbs are always stored, but if the input 414 * operands are normalized. Return the most significant limb of the 415 * result. 416 * 417 * NOTE: The space pointed to by PRODP is overwritten before finished 418 * 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 destination 423 * must be distinct from the multiplier and the multiplicand. 424 */ 425 426 int 427 mpihelp_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 be 444 * 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 else 450 MPN_ZERO(prodp, usize); 451 cy = 0; 452 } else 453 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 from 459 * 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 } else 468 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