mpi/mpih-mul.c

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