--- /dev/null
+/* mpn_mul -- Multiply two natural numbers.
+
+Copyright (C) 1991, 1992 Free Software Foundation, Inc.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2, or (at your option)
+any later version.
+
+The GNU MP Library is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the GNU MP Library; see the file COPYING. If not, write to
+the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+#ifdef GMP_DEBUG /* partain: was DEBUG */
+#define MPN_MUL_VERIFY(res_ptr,res_size,op1_ptr,op1_size,op2_ptr,op2_size) \
+ mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size)
+
+#include <stdio.h>
+static void
+mpn_mul_verify (res_ptr, res_size, op1_ptr, op1_size, op2_ptr, op2_size)
+ mp_ptr res_ptr, op1_ptr, op2_ptr;
+ mp_size res_size, op1_size, op2_size;
+{
+ mp_ptr tmp_ptr;
+ mp_size tmp_size;
+ tmp_ptr = alloca ((op1_size + op2_size) * BYTES_PER_MP_LIMB);
+ if (op1_size >= op2_size)
+ tmp_size = mpn_mul_classic (tmp_ptr,
+ op1_ptr, op1_size, op2_ptr, op2_size);
+ else
+ tmp_size = mpn_mul_classic (tmp_ptr,
+ op2_ptr, op2_size, op1_ptr, op1_size);
+ if (tmp_size != res_size
+ || mpn_cmp (tmp_ptr, res_ptr, tmp_size) != 0)
+ {
+ fprintf (stderr, "GNU MP internal error: Wrong result in mpn_mul.\n");
+ fprintf (stderr, "op1{%d} = ", op1_size); mpn_dump (op1_ptr, op1_size);
+ fprintf (stderr, "op2{%d} = ", op2_size); mpn_dump (op2_ptr, op2_size);
+ abort ();
+ }
+}
+#else
+#define MPN_MUL_VERIFY(a,b,c,d,e,f)
+#endif
+
+/* Multiply the natural numbers u (pointed to by UP, with USIZE limbs)
+ and v (pointed to by VP, with VSIZE limbs), and store the result at
+ PRODP. USIZE + VSIZE limbs are always stored, but if the input
+ operands are normalized, the return value will reflect the true
+ result size (which is either USIZE + VSIZE, or USIZE + VSIZE -1).
+
+ NOTE: The space pointed to by PRODP is overwritten before finished
+ with U and V, so overlap is an error.
+
+ Argument constraints:
+ 1. USIZE >= VSIZE.
+ 2. PRODP != UP and PRODP != VP, i.e. the destination
+ must be distinct from the multiplier and the multiplicand. */
+
+/* If KARATSUBA_THRESHOLD is not already defined, define it to a
+ value which is good on most machines. */
+#ifndef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 8
+#endif
+
+/* The code can't handle KARATSUBA_THRESHOLD smaller than 4. */
+#if KARATSUBA_THRESHOLD < 4
+#undef KARATSUBA_THRESHOLD
+#define KARATSUBA_THRESHOLD 4
+#endif
+
+mp_size
+#ifdef __STDC__
+mpn_mul (mp_ptr prodp,
+ mp_srcptr up, mp_size usize,
+ mp_srcptr vp, mp_size vsize)
+#else
+mpn_mul (prodp, up, usize, vp, vsize)
+ mp_ptr prodp;
+ mp_srcptr up;
+ mp_size usize;
+ mp_srcptr vp;
+ mp_size vsize;
+#endif
+{
+ mp_size n;
+ mp_size prod_size;
+ mp_limb cy;
+
+ if (vsize < KARATSUBA_THRESHOLD)
+ {
+ /* Handle simple cases with traditional multiplication.
+
+ This is the most critical code of the entire function. All
+ multiplies rely on this, both small and huge. Small ones arrive
+ here immediately. Huge ones arrive here as this is the base case
+ for the recursive algorithm below. */
+ mp_size i, j;
+ mp_limb prod_low, prod_high;
+ mp_limb cy_limb;
+ mp_limb v_limb;
+
+ if (vsize == 0)
+ return 0;
+
+ /* Offset UP and PRODP so that the inner loop can be faster. */
+ up += usize;
+ prodp += usize;
+
+ /* Multiply by the first limb in V separately, as the result can
+ be stored (not added) to PROD. We also avoid a loop for zeroing. */
+ v_limb = vp[0];
+ if (v_limb <= 1)
+ {
+ if (v_limb == 1)
+ MPN_COPY (prodp - usize, up - usize, usize);
+ else
+ MPN_ZERO (prodp - usize, usize);
+ cy_limb = 0;
+ }
+ else
+ {
+ cy_limb = 0;
+ j = -usize;
+ do
+ {
+ umul_ppmm (prod_high, prod_low, up[j], v_limb);
+ add_ssaaaa (cy_limb, prodp[j], prod_high, prod_low, 0, cy_limb);
+ j++;
+ }
+ while (j < 0);
+ }
+
+ prodp[0] = cy_limb;
+ prodp++;
+
+ /* For each iteration in the outer loop, multiply one limb from
+ U with one limb from V, and add it to PROD. */
+ for (i = 1; i < vsize; i++)
+ {
+ v_limb = vp[i];
+ if (v_limb <= 1)
+ {
+ cy_limb = 0;
+ if (v_limb == 1)
+ cy_limb = mpn_add (prodp - usize,
+ prodp - usize, usize, up - usize, usize);
+ }
+ else
+ {
+ cy_limb = 0;
+ j = -usize;
+
+ do
+ {
+ umul_ppmm (prod_high, prod_low, up[j], v_limb);
+ add_ssaaaa (cy_limb, prod_low,
+ prod_high, prod_low, 0, cy_limb);
+ add_ssaaaa (cy_limb, prodp[j],
+ cy_limb, prod_low, 0, prodp[j]);
+ j++;
+ }
+ while (j < 0);
+ }
+
+ prodp[0] = cy_limb;
+ prodp++;
+ }
+
+ return usize + vsize - (cy_limb == 0);
+ }
+
+ n = (usize + 1) / 2;
+
+ /* Is USIZE larger than 1.5 times VSIZE? Avoid Karatsuba's algorithm. */
+ if (2 * usize > 3 * vsize)
+ {
+ /* If U has at least twice as many limbs as V. Split U in two
+ pieces, U1 and U0, such that U = U0 + U1*(2**BITS_PER_MP_LIMB)**N,
+ and recursively multiply the two pieces separately with V. */
+
+ mp_size u0_size;
+ mp_ptr tmp;
+ mp_size tmp_size;
+
+ /* V1 (the high part of V) is zero. */
+
+ /* Calculate the length of U0. It is normally equal to n, but
+ of course not for sure. */
+ for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--)
+ ;
+
+ /* Perform (U0 * V). */
+ if (u0_size >= vsize)
+ prod_size = mpn_mul (prodp, up, u0_size, vp, vsize);
+ else
+ prod_size = mpn_mul (prodp, vp, vsize, up, u0_size);
+ MPN_MUL_VERIFY (prodp, prod_size, up, u0_size, vp, vsize);
+
+ /* We have to zero-extend the lower partial product to n limbs,
+ since the mpn_add some lines below expect the first n limbs
+ to be well defined. (This is normally a no-op. It may
+ do something when U1 has many leading 0 limbs.) */
+ while (prod_size < n)
+ prodp[prod_size++] = 0;
+
+ tmp = (mp_ptr) alloca ((usize + vsize - n) * BYTES_PER_MP_LIMB);
+
+ /* Perform (U1 * V). Make sure the first source argument to mpn_mul
+ is not less than the second source argument. */
+ if (vsize <= usize - n)
+ tmp_size = mpn_mul (tmp, up + n, usize - n, vp, vsize);
+ else
+ tmp_size = mpn_mul (tmp, vp, vsize, up + n, usize - n);
+ MPN_MUL_VERIFY (tmp, tmp_size, up + n, usize - n, vp, vsize);
+
+ /* In this addition hides a potentially large copying of TMP. */
+ if (prod_size - n >= tmp_size)
+ cy = mpn_add (prodp + n, prodp + n, prod_size - n, tmp, tmp_size);
+ else
+ cy = mpn_add (prodp + n, tmp, tmp_size, prodp + n, prod_size - n);
+ if (cy)
+ abort (); /* prodp[prod_size] = cy; */
+
+ alloca (0);
+ return tmp_size + n;
+ }
+ else
+ {
+ /* Karatsuba's divide-and-conquer algorithm.
+
+ Split U in two pieces, U1 and U0, such that
+ U = U0 + U1*(B**n),
+ and V in V1 and V0, such that
+ V = V0 + V1*(B**n).
+
+ UV is then computed recursively using the identity
+
+ 2n n n n
+ UV = (B + B )U V + B (U -U )(V -V ) + (B + 1)U V
+ 1 1 1 0 0 1 0 0
+
+ Where B = 2**BITS_PER_MP_LIMB.
+ */
+
+ /* It's possible to decrease the temporary allocation by using the
+ prodp area for temporary storage of the middle term, and doing
+ that recursive multiplication first. (Do this later.) */
+
+ mp_size u0_size;
+ mp_size v0_size;
+ mp_size u0v0_size;
+ mp_size u1v1_size;
+ mp_ptr temp;
+ mp_size temp_size;
+ mp_size utem_size;
+ mp_size vtem_size;
+ mp_ptr ptem;
+ mp_size ptem_size;
+ int negflg;
+ mp_ptr pp;
+
+ pp = (mp_ptr) alloca (4 * n * BYTES_PER_MP_LIMB);
+
+ /* Calculate the lengths of U0 and V0. They are normally equal
+ to n, but of course not for sure. */
+ for (u0_size = n; u0_size > 0 && up[u0_size - 1] == 0; u0_size--)
+ ;
+ for (v0_size = n; v0_size > 0 && vp[v0_size - 1] == 0; v0_size--)
+ ;
+
+ /*** 1. PROD]2n..0] := U0 x V0
+ (Recursive call to mpn_mul may NOT overwrite input operands.)
+ ________________ ________________
+ |________________||____U0 x V0_____| */
+
+ if (u0_size >= v0_size)
+ u0v0_size = mpn_mul (pp, up, u0_size, vp, v0_size);
+ else
+ u0v0_size = mpn_mul (pp, vp, v0_size, up, u0_size);
+ MPN_MUL_VERIFY (pp, u0v0_size, up, u0_size, vp, v0_size);
+
+ /* Zero-extend to 2n limbs. */
+ while (u0v0_size < 2 * n)
+ pp[u0v0_size++] = 0;
+
+
+ /*** 2. PROD]4n..2n] := U1 x V1
+ (Recursive call to mpn_mul may NOT overwrite input operands.)
+ ________________ ________________
+ |_____U1 x V1____||____U0 x V0_____| */
+
+ u1v1_size = mpn_mul (pp + 2*n,
+ up + n, usize - n,
+ vp + n, vsize - n);
+ MPN_MUL_VERIFY (pp + 2*n, u1v1_size,
+ up + n, usize - n, vp + n, vsize - n);
+ prod_size = 2 * n + u1v1_size;
+
+
+ /*** 3. PTEM]2n..0] := (U1-U0) x (V0-V1)
+ (Recursive call to mpn_mul may overwrite input operands.)
+ ________________
+ |_(U1-U0)(V0-V1)_| */
+
+ temp = (mp_ptr) alloca ((2 * n + 1) * BYTES_PER_MP_LIMB);
+ if (usize - n > u0_size
+ || (usize - n == u0_size
+ && mpn_cmp (up + n, up, u0_size) >= 0))
+ {
+ utem_size = usize - n
+ + mpn_sub (temp, up + n, usize - n, up, u0_size);
+ negflg = 0;
+ }
+ else
+ {
+ utem_size = u0_size
+ + mpn_sub (temp, up, u0_size, up + n, usize - n);
+ negflg = 1;
+ }
+ if (vsize - n > v0_size
+ || (vsize - n == v0_size
+ && mpn_cmp (vp + n, vp, v0_size) >= 0))
+ {
+ vtem_size = vsize - n
+ + mpn_sub (temp + n, vp + n, vsize - n, vp, v0_size);
+ negflg ^= 1;
+ }
+ else
+ {
+ vtem_size = v0_size
+ + mpn_sub (temp + n, vp, v0_size, vp + n, vsize - n);
+ /* No change of NEGFLG. */
+ }
+ ptem = (mp_ptr) alloca (2 * n * BYTES_PER_MP_LIMB);
+ if (utem_size >= vtem_size)
+ ptem_size = mpn_mul (ptem, temp, utem_size, temp + n, vtem_size);
+ else
+ ptem_size = mpn_mul (ptem, temp + n, vtem_size, temp, utem_size);
+ MPN_MUL_VERIFY (ptem, ptem_size, temp, utem_size, temp + n, vtem_size);
+
+ /*** 4. TEMP]2n..0] := PROD]2n..0] + PROD]4n..2n]
+ ________________
+ |_____U1 x V1____|
+ ________________
+ |_____U0_x_V0____| */
+
+ cy = mpn_add (temp, pp, 2*n, pp + 2*n, u1v1_size);
+ if (cy != 0)
+ {
+ temp[2*n] = cy;
+ temp_size = 2*n + 1;
+ }
+ else
+ {
+ /* Normalize temp. pp[2*n-1] might have been zero in the
+ mpn_add call above, and thus temp might be unnormalized. */
+ for (temp_size = 2*n; temp_size > 0 && temp[temp_size - 1] == 0;
+ temp_size--)
+ ;
+ }
+
+ if (prod_size - n >= temp_size)
+ cy = mpn_add (pp + n, pp + n, prod_size - n, temp, temp_size);
+ else
+ {
+ /* This is a weird special case that should not happen (often)! */
+ cy = mpn_add (pp + n, temp, temp_size, pp + n, prod_size - n);
+ prod_size = temp_size + n;
+ }
+ if (cy != 0)
+ {
+ pp[prod_size] = cy;
+ prod_size++;
+ }
+#ifdef GMP_DEBUG /* partain: was DEBUG */
+ if (prod_size > 4 * n)
+ abort();
+#endif
+ if (negflg)
+ prod_size = prod_size
+ + mpn_sub (pp + n, pp + n, prod_size - n, ptem, ptem_size);
+ else
+ {
+ if (prod_size - n < ptem_size)
+ abort();
+ cy = mpn_add (pp + n, pp + n, prod_size - n, ptem, ptem_size);
+ if (cy != 0)
+ {
+ pp[prod_size] = cy;
+ prod_size++;
+#ifdef GMP_DEBUG /* partain: was DEBUG */
+ if (prod_size > 4 * n)
+ abort();
+#endif
+ }
+ }
+
+ MPN_COPY (prodp, pp, prod_size);
+ alloca (0);
+ return prod_size;
+ }
+}
projects / ghc-hetmet.git / blobdiff