stdlib/mul_n.c

   1 /* mpn_mul_n -- Multiply two natural numbers of length n.
   2
   3 Copyright (C) 1991, 1992, 1993, 1994, 1996 Free Software Foundation, Inc.
   4
   5 This file is part of the GNU MP Library.
   6
   7 The GNU MP Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The GNU MP Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include <gmp.h>
  23 #include "gmp-impl.h"
  24
  25 /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
  26    both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
  27    always stored.  Return the most significant limb.
  28
  29    Argument constraints:
  30    1. PRODP != UP and PRODP != VP, i.e. the destination
  31       must be distinct from the multiplier and the multiplicand.  */
  32
  33 /* If KARATSUBA_THRESHOLD is not already defined, define it to a
  34    value which is good on most machines.  */
  35 #ifndef KARATSUBA_THRESHOLD
  36 #define KARATSUBA_THRESHOLD 32
  37 #endif
  38
  39 /* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
  40 #if KARATSUBA_THRESHOLD < 2
  41 #undef KARATSUBA_THRESHOLD
  42 #define KARATSUBA_THRESHOLD 2
  43 #endif
  44
  45 /* Handle simple cases with traditional multiplication.
  46
  47    This is the most critical code of multiplication.  All multiplies rely
  48    on this, both small and huge.  Small ones arrive here immediately.  Huge
  49    ones arrive here as this is the base case for Karatsuba's recursive
  50    algorithm below.  */
  51
  52 void
  53 #if __STDC__
  54 impn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
  55 #else
  56 impn_mul_n_basecase (prodp, up, vp, size)
  57      mp_ptr prodp;
  58      mp_srcptr up;
  59      mp_srcptr vp;
  60      mp_size_t size;
  61 #endif
  62 {
  63   mp_size_t i;
  64   mp_limb_t cy_limb;
  65   mp_limb_t v_limb;
  66
  67   /* Multiply by the first limb in V separately, as the result can be
  68      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
  69   v_limb = vp[0];
  70   if (v_limb <= 1)
  71     {
  72       if (v_limb == 1)
  73         MPN_COPY (prodp, up, size);
  74       else
  75         MPN_ZERO (prodp, size);
  76       cy_limb = 0;
  77     }
  78   else
  79     cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
  80
  81   prodp[size] = cy_limb;
  82   prodp++;
  83
  84   /* For each iteration in the outer loop, multiply one limb from
  85      U with one limb from V, and add it to PROD.  */
  86   for (i = 1; i < size; i++)
  87     {
  88       v_limb = vp[i];
  89       if (v_limb <= 1)
  90         {
  91           cy_limb = 0;
  92           if (v_limb == 1)
  93             cy_limb = mpn_add_n (prodp, prodp, up, size);
  94         }
  95       else
  96         cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
  97
  98       prodp[size] = cy_limb;
  99       prodp++;
 100     }
 101 }
 102
 103 void
 104 #if __STDC__
 105 impn_mul_n (mp_ptr prodp,
 106              mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
 107 #else
 108 impn_mul_n (prodp, up, vp, size, tspace)
 109      mp_ptr prodp;
 110      mp_srcptr up;
 111      mp_srcptr vp;
 112      mp_size_t size;
 113      mp_ptr tspace;
 114 #endif
 115 {
 116   if ((size & 1) != 0)
 117     {
 118       /* The size is odd, the code code below doesn't handle that.
 119          Multiply the least significant (size - 1) limbs with a recursive
 120          call, and handle the most significant limb of S1 and S2
 121          separately.  */
 122       /* A slightly faster way to do this would be to make the Karatsuba
 123          code below behave as if the size were even, and let it check for
 124          odd size in the end.  I.e., in essence move this code to the end.
 125          Doing so would save us a recursive call, and potentially make the
 126          stack grow a lot less.  */
 127
 128       mp_size_t esize = size - 1;       /* even size */
 129       mp_limb_t cy_limb;
 130
 131       MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
 132       cy_limb = mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
 133       prodp[esize + esize] = cy_limb;
 134       cy_limb = mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
 135
 136       prodp[esize + size] = cy_limb;
 137     }
 138   else
 139     {
 140       /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
 141
 142          Split U in two pieces, U1 and U0, such that
 143          U = U0 + U1*(B**n),
 144          and V in V1 and V0, such that
 145          V = V0 + V1*(B**n).
 146
 147          UV is then computed recursively using the identity
 148
 149                 2n   n          n                     n
 150          UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
 151                         1 1        1  0   0  1              0 0
 152
 153          Where B = 2**BITS_PER_MP_LIMB.  */
 154
 155       mp_size_t hsize = size >> 1;
 156       mp_limb_t cy;
 157       int negflg;
 158
 159       /*** Product H.    ________________  ________________
 160                         |_____U1 x V1____||____U0 x V0_____|  */
 161       /* Put result in upper part of PROD and pass low part of TSPACE
 162          as new TSPACE.  */
 163       MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
 164
 165       /*** Product M.    ________________
 166                         |_(U1-U0)(V0-V1)_|  */
 167       if (mpn_cmp (up + hsize, up, hsize) >= 0)
 168         {
 169           mpn_sub_n (prodp, up + hsize, up, hsize);
 170           negflg = 0;
 171         }
 172       else
 173         {
 174           mpn_sub_n (prodp, up, up + hsize, hsize);
 175           negflg = 1;
 176         }
 177       if (mpn_cmp (vp + hsize, vp, hsize) >= 0)
 178         {
 179           mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
 180           negflg ^= 1;
 181         }
 182       else
 183         {
 184           mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
 185           /* No change of NEGFLG.  */
 186         }
 187       /* Read temporary operands from low part of PROD.
 188          Put result in low part of TSPACE using upper part of TSPACE
 189          as new TSPACE.  */
 190       MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
 191
 192       /*** Add/copy product H.  */
 193       MPN_COPY (prodp + hsize, prodp + size, hsize);
 194       cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
 195
 196       /*** Add product M (if NEGFLG M is a negative number).  */
 197       if (negflg)
 198         cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
 199       else
 200         cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
 201
 202       /*** Product L.    ________________  ________________
 203                         |________________||____U0 x V0_____|  */
 204       /* Read temporary operands from low part of PROD.
 205          Put result in low part of TSPACE using upper part of TSPACE
 206          as new TSPACE.  */
 207       MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
 208
 209       /*** Add/copy Product L (twice).  */
 210
 211       cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
 212       if (cy)
 213         mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
 214
 215       MPN_COPY (prodp, tspace, hsize);
 216       cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
 217       if (cy)
 218         mpn_add_1 (prodp + size, prodp + size, size, 1);
 219     }
 220 }
 221
 222 void
 223 #if __STDC__
 224 impn_sqr_n_basecase (mp_ptr prodp, mp_srcptr up, mp_size_t size)
 225 #else
 226 impn_sqr_n_basecase (prodp, up, size)
 227      mp_ptr prodp;
 228      mp_srcptr up;
 229      mp_size_t size;
 230 #endif
 231 {
 232   mp_size_t i;
 233   mp_limb_t cy_limb;
 234   mp_limb_t v_limb;
 235
 236   /* Multiply by the first limb in V separately, as the result can be
 237      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
 238   v_limb = up[0];
 239   if (v_limb <= 1)
 240     {
 241       if (v_limb == 1)
 242         MPN_COPY (prodp, up, size);
 243       else
 244         MPN_ZERO (prodp, size);
 245       cy_limb = 0;
 246     }
 247   else
 248     cy_limb = mpn_mul_1 (prodp, up, size, v_limb);
 249
 250   prodp[size] = cy_limb;
 251   prodp++;
 252
 253   /* For each iteration in the outer loop, multiply one limb from
 254      U with one limb from V, and add it to PROD.  */
 255   for (i = 1; i < size; i++)
 256     {
 257       v_limb = up[i];
 258       if (v_limb <= 1)
 259         {
 260           cy_limb = 0;
 261           if (v_limb == 1)
 262             cy_limb = mpn_add_n (prodp, prodp, up, size);
 263         }
 264       else
 265         cy_limb = mpn_addmul_1 (prodp, up, size, v_limb);
 266
 267       prodp[size] = cy_limb;
 268       prodp++;
 269     }
 270 }
 271
 272 void
 273 #if __STDC__
 274 impn_sqr_n (mp_ptr prodp,
 275              mp_srcptr up, mp_size_t size, mp_ptr tspace)
 276 #else
 277 impn_sqr_n (prodp, up, size, tspace)
 278      mp_ptr prodp;
 279      mp_srcptr up;
 280      mp_size_t size;
 281      mp_ptr tspace;
 282 #endif
 283 {
 284   if ((size & 1) != 0)
 285     {
 286       /* The size is odd, the code code below doesn't handle that.
 287          Multiply the least significant (size - 1) limbs with a recursive
 288          call, and handle the most significant limb of S1 and S2
 289          separately.  */
 290       /* A slightly faster way to do this would be to make the Karatsuba
 291          code below behave as if the size were even, and let it check for
 292          odd size in the end.  I.e., in essence move this code to the end.
 293          Doing so would save us a recursive call, and potentially make the
 294          stack grow a lot less.  */
 295
 296       mp_size_t esize = size - 1;       /* even size */
 297       mp_limb_t cy_limb;
 298
 299       MPN_SQR_N_RECURSE (prodp, up, esize, tspace);
 300       cy_limb = mpn_addmul_1 (prodp + esize, up, esize, up[esize]);
 301       prodp[esize + esize] = cy_limb;
 302       cy_limb = mpn_addmul_1 (prodp + esize, up, size, up[esize]);
 303
 304       prodp[esize + size] = cy_limb;
 305     }
 306   else
 307     {
 308       mp_size_t hsize = size >> 1;
 309       mp_limb_t cy;
 310
 311       /*** Product H.    ________________  ________________
 312                         |_____U1 x U1____||____U0 x U0_____|  */
 313       /* Put result in upper part of PROD and pass low part of TSPACE
 314          as new TSPACE.  */
 315       MPN_SQR_N_RECURSE (prodp + size, up + hsize, hsize, tspace);
 316
 317       /*** Product M.    ________________
 318                         |_(U1-U0)(U0-U1)_|  */
 319       if (mpn_cmp (up + hsize, up, hsize) >= 0)
 320         {
 321           mpn_sub_n (prodp, up + hsize, up, hsize);
 322         }
 323       else
 324         {
 325           mpn_sub_n (prodp, up, up + hsize, hsize);
 326         }
 327
 328       /* Read temporary operands from low part of PROD.
 329          Put result in low part of TSPACE using upper part of TSPACE
 330          as new TSPACE.  */
 331       MPN_SQR_N_RECURSE (tspace, prodp, hsize, tspace + size);
 332
 333       /*** Add/copy product H.  */
 334       MPN_COPY (prodp + hsize, prodp + size, hsize);
 335       cy = mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
 336
 337       /*** Add product M (if NEGFLG M is a negative number).  */
 338       cy -= mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
 339
 340       /*** Product L.    ________________  ________________
 341                         |________________||____U0 x U0_____|  */
 342       /* Read temporary operands from low part of PROD.
 343          Put result in low part of TSPACE using upper part of TSPACE
 344          as new TSPACE.  */
 345       MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
 346
 347       /*** Add/copy Product L (twice).  */
 348
 349       cy += mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
 350       if (cy)
 351         mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
 352
 353       MPN_COPY (prodp, tspace, hsize);
 354       cy = mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
 355       if (cy)
 356         mpn_add_1 (prodp + size, prodp + size, size, 1);
 357     }
 358 }
 359
 360 /* This should be made into an inline function in gmp.h.  */
 361 void
 362 #if __STDC__
 363 mpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
 364 #else
 365 mpn_mul_n (prodp, up, vp, size)
 366      mp_ptr prodp;
 367      mp_srcptr up;
 368      mp_srcptr vp;
 369      mp_size_t size;
 370 #endif
 371 {
 372   TMP_DECL (marker);
 373   TMP_MARK (marker);
 374   if (up == vp)
 375     {
 376       if (size < KARATSUBA_THRESHOLD)
 377         {
 378           impn_sqr_n_basecase (prodp, up, size);
 379         }
 380       else
 381         {
 382           mp_ptr tspace;
 383           tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
 384           impn_sqr_n (prodp, up, size, tspace);
 385         }
 386     }
 387   else
 388     {
 389       if (size < KARATSUBA_THRESHOLD)
 390         {
 391           impn_mul_n_basecase (prodp, up, vp, size);
 392         }
 393       else
 394         {
 395           mp_ptr tspace;
 396           tspace = (mp_ptr) TMP_ALLOC (2 * size * BYTES_PER_MP_LIMB);
 397           impn_mul_n (prodp, up, vp, size, tspace);
 398         }
 399     }
 400   TMP_FREE (marker);
 401 }