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>