rts/gmp/mpn/generic/mul.c

   1 /* mpn_mul -- Multiply two natural numbers.
   2
   3    THE HELPER FUNCTIONS IN THIS FILE (meaning everything except mpn_mul)
   4    ARE INTERNAL FUNCTIONS WITH MUTABLE INTERFACES.  IT IS ONLY SAFE TO REACH
   5    THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST GUARANTEED
   6    THAT THEY'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
   7
   8
   9 Copyright (C) 1991, 1993, 1994, 1996, 1997, 1999, 2000 Free Software
  10 Foundation, Inc.
  11
  12 This file is part of the GNU MP Library.
  13
  14 The GNU MP Library is free software; you can redistribute it and/or modify
  15 it under the terms of the GNU Lesser General Public License as published by
  16 the Free Software Foundation; either version 2.1 of the License, or (at your
  17 option) any later version.
  18
  19 The GNU MP Library is distributed in the hope that it will be useful, but
  20 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  21 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  22 License for more details.
  23
  24 You should have received a copy of the GNU Lesser General Public License
  25 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  26 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  27 MA 02111-1307, USA. */
  28
  29 #include "gmp.h"
  30 #include "gmp-impl.h"
  31
  32 /* Multiply the natural numbers u (pointed to by UP, with UN limbs) and v
  33    (pointed to by VP, with VN limbs), and store the result at PRODP.  The
  34    result is UN + VN limbs.  Return the most significant limb of the result.
  35
  36    NOTE: The space pointed to by PRODP is overwritten before finished with U
  37    and V, so overlap is an error.
  38
  39    Argument constraints:
  40    1. UN >= VN.
  41    2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from
  42       the multiplier and the multiplicand.  */
  43
  44 void
  45 #if __STDC__
  46 mpn_sqr_n (mp_ptr prodp,
  47          mp_srcptr up, mp_size_t un)
  48 #else
  49 mpn_sqr_n (prodp, up, un)
  50      mp_ptr prodp;
  51      mp_srcptr up;
  52      mp_size_t un;
  53 #endif
  54 {
  55   if (un < KARATSUBA_SQR_THRESHOLD)
  56     { /* plain schoolbook multiplication */
  57       if (un == 0)
  58         return;
  59       mpn_sqr_basecase (prodp, up, un);
  60     }
  61   else if (un < TOOM3_SQR_THRESHOLD)
  62     { /* karatsuba multiplication */
  63       mp_ptr tspace;
  64       TMP_DECL (marker);
  65       TMP_MARK (marker);
  66       tspace = (mp_ptr) TMP_ALLOC (2 * (un + BITS_PER_MP_LIMB) * BYTES_PER_MP_LIMB);
  67       mpn_kara_sqr_n (prodp, up, un, tspace);
  68       TMP_FREE (marker);
  69     }
  70 #if WANT_FFT || TUNE_PROGRAM_BUILD
  71   else if (un < FFT_SQR_THRESHOLD)
  72 #else
  73   else
  74 #endif
  75     { /* toom3 multiplication */
  76       mp_ptr tspace;
  77       TMP_DECL (marker);
  78       TMP_MARK (marker);
  79       tspace = (mp_ptr) TMP_ALLOC (2 * (un + BITS_PER_MP_LIMB) * BYTES_PER_MP_LIMB);
  80       mpn_toom3_sqr_n (prodp, up, un, tspace);
  81       TMP_FREE (marker);
  82     }
  83 #if WANT_FFT || TUNE_PROGRAM_BUILD
  84   else
  85     {
  86       /* schoenhage multiplication */
  87       mpn_mul_fft_full (prodp, up, un, up, un);
  88     }
  89 #endif
  90 }
  91
  92 mp_limb_t
  93 #if __STDC__
  94 mpn_mul (mp_ptr prodp,
  95          mp_srcptr up, mp_size_t un,
  96          mp_srcptr vp, mp_size_t vn)
  97 #else
  98 mpn_mul (prodp, up, un, vp, vn)
  99      mp_ptr prodp;
 100      mp_srcptr up;
 101      mp_size_t un;
 102      mp_srcptr vp;
 103      mp_size_t vn;
 104 #endif
 105 {
 106   mp_size_t l;
 107   mp_limb_t c;
 108
 109   if (up == vp && un == vn)
 110     {
 111       mpn_sqr_n (prodp, up, un);
 112       return prodp[2 * un - 1];
 113     }
 114
 115   if (vn < KARATSUBA_MUL_THRESHOLD)
 116     { /* long multiplication */
 117       mpn_mul_basecase (prodp, up, un, vp, vn);
 118       return prodp[un + vn - 1];
 119     }
 120
 121   mpn_mul_n (prodp, up, vp, vn);
 122   if (un != vn)
 123     { mp_limb_t t;
 124       mp_ptr ws;
 125       TMP_DECL (marker);
 126       TMP_MARK (marker);
 127
 128       prodp += vn;
 129       l = vn;
 130       up += vn;
 131       un -= vn;
 132
 133       if (un < vn)
 134         {
 135           /* Swap u's and v's. */
 136           MPN_SRCPTR_SWAP (up,un, vp,vn);
 137         }
 138
 139       ws = (mp_ptr) TMP_ALLOC (((vn >= KARATSUBA_MUL_THRESHOLD ? vn : un) + vn)
 140                                * BYTES_PER_MP_LIMB);
 141
 142       t = 0;
 143       while (vn >= KARATSUBA_MUL_THRESHOLD)
 144         {
 145           mpn_mul_n (ws, up, vp, vn);
 146           if (l <= 2*vn)
 147             {
 148               t += mpn_add_n (prodp, prodp, ws, l);
 149               if (l != 2*vn)
 150                 {
 151                   t = mpn_add_1 (prodp + l, ws + l, 2*vn - l, t);
 152                   l = 2*vn;
 153                 }
 154             }
 155           else
 156             {
 157               c = mpn_add_n (prodp, prodp, ws, 2*vn);
 158               t += mpn_add_1 (prodp + 2*vn, prodp + 2*vn, l - 2*vn, c);
 159             }
 160           prodp += vn;
 161           l -= vn;
 162           up += vn;
 163           un -= vn;
 164           if (un < vn)
 165             {
 166               /* Swap u's and v's. */
 167               MPN_SRCPTR_SWAP (up,un, vp,vn);
 168             }
 169         }
 170
 171       if (vn)
 172         {
 173           mpn_mul_basecase (ws, up, un, vp, vn);
 174           if (l <= un + vn)
 175             {
 176               t += mpn_add_n (prodp, prodp, ws, l);
 177               if (l != un + vn)
 178                 t = mpn_add_1 (prodp + l, ws + l, un + vn - l, t);
 179             }
 180           else
 181             {
 182               c = mpn_add_n (prodp, prodp, ws, un + vn);
 183               t += mpn_add_1 (prodp + un + vn, prodp + un + vn, l - un - vn, c);
 184             }
 185         }
 186
 187     TMP_FREE (marker);
 188   }
 189   return prodp[un + vn - 1];
 190 }