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