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