lib/gmp/mpn/mul.c

   1 /* mpn_mul -- Multiply two natural numbers.
   2
   3 Copyright (C) 1991, 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, with USIZE limbs)
  27    and v (pointed to by VP, with VSIZE limbs), and store the result at
  28    PRODP.  USIZE + VSIZE limbs are always stored, but if the input
  29    operands are normalized.  Return the most significant limb of the
  30    result.
  31
  32    NOTE: The space pointed to by PRODP is overwritten before finished
  33    with U and V, so overlap is an error.
  34
  35    Argument constraints:
  36    1. USIZE >= VSIZE.
  37    2. PRODP != UP and PRODP != VP, i.e. the destination
  38       must be distinct from the multiplier and the multiplicand.  */
  39
  40 /* If KARATSUBA_THRESHOLD is not already defined, define it to a
  41    value which is good on most machines.  */
  42 #ifndef KARATSUBA_THRESHOLD
  43 #define KARATSUBA_THRESHOLD 32
  44 #endif
  45
  46 mp_limb_t
  47 #if __STDC__
  48 mpn_mul (mp_ptr prodp,
  49          mp_srcptr up, mp_size_t usize,
  50          mp_srcptr vp, mp_size_t vsize)
  51 #else
  52 mpn_mul (prodp, up, usize, vp, vsize)
  53      mp_ptr prodp;
  54      mp_srcptr up;
  55      mp_size_t usize;
  56      mp_srcptr vp;
  57      mp_size_t vsize;
  58 #endif
  59 {
  60   mp_ptr prod_endp = prodp + usize + vsize - 1;
  61   mp_limb_t cy;
  62   mp_ptr tspace;
  63   TMP_DECL (marker);
  64
  65   if (vsize < KARATSUBA_THRESHOLD)
  66     {
  67       /* Handle simple cases with traditional multiplication.
  68
  69          This is the most critical code of the entire function.  All
  70          multiplies rely on this, both small and huge.  Small ones arrive
  71          here immediately.  Huge ones arrive here as this is the base case
  72          for Karatsuba's recursive algorithm below.  */
  73       mp_size_t i;
  74       mp_limb_t cy_limb;
  75       mp_limb_t v_limb;
  76
  77       if (vsize == 0)
  78         return 0;
  79
  80       /* Multiply by the first limb in V separately, as the result can be
  81          stored (not added) to PROD.  We also avoid a loop for zeroing.  */
  82       v_limb = vp[0];
  83       if (v_limb <= 1)
  84         {
  85           if (v_limb == 1)
  86             MPN_COPY (prodp, up, usize);
  87           else
  88             MPN_ZERO (prodp, usize);
  89           cy_limb = 0;
  90         }
  91       else
  92         cy_limb = mpn_mul_1 (prodp, up, usize, v_limb);
  93
  94       prodp[usize] = cy_limb;
  95       prodp++;
  96
  97       /* For each iteration in the outer loop, multiply one limb from
  98          U with one limb from V, and add it to PROD.  */
  99       for (i = 1; i < vsize; i++)
 100         {
 101           v_limb = vp[i];
 102           if (v_limb <= 1)
 103             {
 104               cy_limb = 0;
 105               if (v_limb == 1)
 106                 cy_limb = mpn_add_n (prodp, prodp, up, usize);
 107             }
 108           else
 109             cy_limb = mpn_addmul_1 (prodp, up, usize, v_limb);
 110
 111           prodp[usize] = cy_limb;
 112           prodp++;
 113         }
 114       return cy_limb;
 115     }
 116
 117   TMP_MARK (marker);
 118
 119   tspace = (mp_ptr) TMP_ALLOC (2 * vsize * BYTES_PER_MP_LIMB);
 120   MPN_MUL_N_RECURSE (prodp, up, vp, vsize, tspace);
 121
 122   prodp += vsize;
 123   up += vsize;
 124   usize -= vsize;
 125   if (usize >= vsize)
 126     {
 127       mp_ptr tp = (mp_ptr) TMP_ALLOC (2 * vsize * BYTES_PER_MP_LIMB);
 128       do
 129         {
 130           MPN_MUL_N_RECURSE (tp, up, vp, vsize, tspace);
 131           cy = mpn_add_n (prodp, prodp, tp, vsize);
 132           mpn_add_1 (prodp + vsize, tp + vsize, vsize, cy);
 133           prodp += vsize;
 134           up += vsize;
 135           usize -= vsize;
 136         }
 137       while (usize >= vsize);
 138     }
 139
 140   /* True: usize < vsize.  */
 141
 142   /* Make life simple: Recurse.  */
 143
 144   if (usize != 0)
 145     {
 146       mpn_mul (tspace, vp, vsize, up, usize);
 147       cy = mpn_add_n (prodp, prodp, tspace, vsize);
 148       mpn_add_1 (prodp + vsize, tspace + vsize, usize, cy);
 149     }
 150
 151   TMP_FREE (marker);
 152   return *prod_endp;
 153 }