[BACK]Return to mul.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / mpn / generic

Annotation of OpenXM_contrib/gmp/mpn/generic/mul.c, Revision 1.1.1.2

1.1       maekawa     1: /* mpn_mul -- Multiply two natural numbers.
                      2:
1.1.1.2 ! maekawa     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.
1.1       maekawa    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
1.1.1.2 ! maekawa    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
1.1       maekawa    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
1.1.1.2 ! maekawa    21: or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
1.1       maekawa    22: License for more details.
                     23:
1.1.1.2 ! maekawa    24: You should have received a copy of the GNU Lesser General Public License
1.1       maekawa    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:
1.1.1.2 ! maekawa    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.
1.1       maekawa    35:
1.1.1.2 ! maekawa    36:    NOTE: The space pointed to by PRODP is overwritten before finished with U
        !            37:    and V, so overlap is an error.
1.1       maekawa    38:
                     39:    Argument constraints:
1.1.1.2 ! maekawa    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:     }
1.1       maekawa    89: #endif
1.1.1.2 ! maekawa    90: }
1.1       maekawa    91:
                     92: mp_limb_t
                     93: #if __STDC__
                     94: mpn_mul (mp_ptr prodp,
1.1.1.2 ! maekawa    95:         mp_srcptr up, mp_size_t un,
        !            96:         mp_srcptr vp, mp_size_t vn)
1.1       maekawa    97: #else
1.1.1.2 ! maekawa    98: mpn_mul (prodp, up, un, vp, vn)
1.1       maekawa    99:      mp_ptr prodp;
                    100:      mp_srcptr up;
1.1.1.2 ! maekawa   101:      mp_size_t un;
1.1       maekawa   102:      mp_srcptr vp;
1.1.1.2 ! maekawa   103:      mp_size_t vn;
1.1       maekawa   104: #endif
                    105: {
1.1.1.2 ! maekawa   106:   mp_size_t l;
        !           107:   mp_limb_t c;
1.1       maekawa   108:
1.1.1.2 ! maekawa   109:   if (up == vp && un == vn)
1.1       maekawa   110:     {
1.1.1.2 ! maekawa   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;
1.1       maekawa   132:
1.1.1.2 ! maekawa   133:       if (un < vn)
1.1       maekawa   134:        {
1.1.1.2 ! maekawa   135:          /* Swap u's and v's. */
        !           136:           MPN_SRCPTR_SWAP (up,un, vp,vn);
1.1       maekawa   137:        }
                    138:
1.1.1.2 ! maekawa   139:       ws = (mp_ptr) TMP_ALLOC (((vn >= KARATSUBA_MUL_THRESHOLD ? vn : un) + vn)
        !           140:                               * BYTES_PER_MP_LIMB);
1.1       maekawa   141:
1.1.1.2 ! maekawa   142:       t = 0;
        !           143:       while (vn >= KARATSUBA_MUL_THRESHOLD)
1.1       maekawa   144:        {
1.1.1.2 ! maekawa   145:          mpn_mul_n (ws, up, vp, vn);
        !           146:          if (l <= 2*vn)
1.1       maekawa   147:            {
1.1.1.2 ! maekawa   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:                }
1.1       maekawa   154:            }
                    155:          else
1.1.1.2 ! maekawa   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:            }
1.1       maekawa   169:        }
                    170:
1.1.1.2 ! maekawa   171:       if (vn)
1.1       maekawa   172:        {
1.1.1.2 ! maekawa   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:            }
1.1       maekawa   185:        }
                    186:
1.1.1.2 ! maekawa   187:     TMP_FREE (marker);
        !           188:   }
        !           189:   return prodp[un + vn - 1];
1.1       maekawa   190: }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>