Back to index

glibc  2.9
mul.c
Go to the documentation of this file.
00001 /* mpn_mul -- Multiply two natural numbers.
00002 
00003 Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
00004 
00005 This file is part of the GNU MP Library.
00006 
00007 The GNU MP Library is free software; you can redistribute it and/or modify
00008 it under the terms of the GNU Lesser General Public License as published by
00009 the Free Software Foundation; either version 2.1 of the License, or (at your
00010 option) any later version.
00011 
00012 The GNU MP Library is distributed in the hope that it will be useful, but
00013 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
00014 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00015 License for more details.
00016 
00017 You should have received a copy of the GNU Lesser General Public License
00018 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
00019 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00020 MA 02111-1307, USA. */
00021 
00022 #include <gmp.h>
00023 #include "gmp-impl.h"
00024 
00025 /* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
00026    and v (pointed to by VP, with VSIZE limbs), and store the result at
00027    PRODP.  USIZE + VSIZE limbs are always stored, but if the input
00028    operands are normalized.  Return the most significant limb of the
00029    result.
00030 
00031    NOTE: The space pointed to by PRODP is overwritten before finished
00032    with U and V, so overlap is an error.
00033 
00034    Argument constraints:
00035    1. USIZE >= VSIZE.
00036    2. PRODP != UP and PRODP != VP, i.e. the destination
00037       must be distinct from the multiplier and the multiplicand.  */
00038 
00039 /* If KARATSUBA_THRESHOLD is not already defined, define it to a
00040    value which is good on most machines.  */
00041 #ifndef KARATSUBA_THRESHOLD
00042 #define KARATSUBA_THRESHOLD 32
00043 #endif
00044 
00045 mp_limb_t
00046 #if __STDC__
00047 mpn_mul (mp_ptr prodp,
00048         mp_srcptr up, mp_size_t usize,
00049         mp_srcptr vp, mp_size_t vsize)
00050 #else
00051 mpn_mul (prodp, up, usize, vp, vsize)
00052      mp_ptr prodp;
00053      mp_srcptr up;
00054      mp_size_t usize;
00055      mp_srcptr vp;
00056      mp_size_t vsize;
00057 #endif
00058 {
00059   mp_ptr prod_endp = prodp + usize + vsize - 1;
00060   mp_limb_t cy;
00061   mp_ptr tspace;
00062   TMP_DECL (marker);
00063 
00064   if (vsize < KARATSUBA_THRESHOLD)
00065     {
00066       /* Handle simple cases with traditional multiplication.
00067 
00068         This is the most critical code of the entire function.  All
00069         multiplies rely on this, both small and huge.  Small ones arrive
00070         here immediately.  Huge ones arrive here as this is the base case
00071         for Karatsuba's recursive algorithm below.  */
00072       mp_size_t i;
00073       mp_limb_t cy_limb;
00074       mp_limb_t v_limb;
00075 
00076       if (vsize == 0)
00077        return 0;
00078 
00079       /* Multiply by the first limb in V separately, as the result can be
00080         stored (not added) to PROD.  We also avoid a loop for zeroing.  */
00081       v_limb = vp[0];
00082       if (v_limb <= 1)
00083        {
00084          if (v_limb == 1)
00085            MPN_COPY (prodp, up, usize);
00086          else
00087            MPN_ZERO (prodp, usize);
00088          cy_limb = 0;
00089        }
00090       else
00091        cy_limb = mpn_mul_1 (prodp, up, usize, v_limb);
00092 
00093       prodp[usize] = cy_limb;
00094       prodp++;
00095 
00096       /* For each iteration in the outer loop, multiply one limb from
00097         U with one limb from V, and add it to PROD.  */
00098       for (i = 1; i < vsize; i++)
00099        {
00100          v_limb = vp[i];
00101          if (v_limb <= 1)
00102            {
00103              cy_limb = 0;
00104              if (v_limb == 1)
00105               cy_limb = mpn_add_n (prodp, prodp, up, usize);
00106            }
00107          else
00108            cy_limb = mpn_addmul_1 (prodp, up, usize, v_limb);
00109 
00110          prodp[usize] = cy_limb;
00111          prodp++;
00112        }
00113       return cy_limb;
00114     }
00115 
00116   TMP_MARK (marker);
00117 
00118   tspace = (mp_ptr) TMP_ALLOC (2 * vsize * BYTES_PER_MP_LIMB);
00119   MPN_MUL_N_RECURSE (prodp, up, vp, vsize, tspace);
00120 
00121   prodp += vsize;
00122   up += vsize;
00123   usize -= vsize;
00124   if (usize >= vsize)
00125     {
00126       mp_ptr tp = (mp_ptr) TMP_ALLOC (2 * vsize * BYTES_PER_MP_LIMB);
00127       do
00128        {
00129          MPN_MUL_N_RECURSE (tp, up, vp, vsize, tspace);
00130          cy = mpn_add_n (prodp, prodp, tp, vsize);
00131          mpn_add_1 (prodp + vsize, tp + vsize, vsize, cy);
00132          prodp += vsize;
00133          up += vsize;
00134          usize -= vsize;
00135        }
00136       while (usize >= vsize);
00137     }
00138 
00139   /* True: usize < vsize.  */
00140 
00141   /* Make life simple: Recurse.  */
00142 
00143   if (usize != 0)
00144     {
00145       mpn_mul (tspace, vp, vsize, up, usize);
00146       cy = mpn_add_n (prodp, prodp, tspace, vsize);
00147       mpn_add_1 (prodp + vsize, tspace + vsize, usize, cy);
00148     }
00149 
00150   TMP_FREE (marker);
00151   return *prod_endp;
00152 }