/* _mpn_mul -- Multiply two natural numbers.
Copyright (C) 1991 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 <stdlib.h>
#include "gmp.h"
#include "gmp-impl.h"
#include "longlong.h"
/* disabled for now, _mpn_mul_classic has vanished somehow... */
#if defined(DEBUG) && 0
#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); /* debug only */
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 (res_ptr, res_size, tmp_ptr, tmp_size) != 0)
{
fprintf (stderr, "GNU MP internal error: Wrong result in _mpn_mul.\n");
fprintf (stderr, "op1{%d} = ", op1_size); _mpn_print (op1_ptr, op1_size);
fprintf (stderr, "op2{%d} = ", op2_size); _mpn_print (op2_ptr, op2_size);
abort();
}
}
_mpn_print (mp_ptr p, mp_size size)
{
mp_size i;
for (i = size - 1; i >= 0; i--)
printf ("%08X", p[i]);
puts ("");
}
#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. */
#define KARATSUBA_THRESHOLD 100000000 /* Ridiculous large, but we will never do alloca ! */
/* 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 _mpn_mul (mp_ptr prodp,
mp_srcptr up, mp_size usize,
mp_srcptr vp, mp_size vsize)
{
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. */
Int i, j;
mp_limb prod_low, prod_high;
mp_limb cy_limb;
mp_limb v_limb, c;
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];
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[j] = 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];
cy_limb = 0;
j = -usize;
/* Inner loops. Simulate the carry flag by jumping between
these loops. The first is used when there was no carry
in the previois iteration; the second when there was carry. */
do
{
umul_ppmm (prod_high, prod_low, up[j], v_limb);
add_ssaaaa (cy_limb, prod_low, prod_high, prod_low, 0, cy_limb);
c = prodp[j];
prod_low += c;
prodp[j] = prod_low;
if (prod_low < c)
goto cy_loop;
ncy_loop:
j++;
}
while (j < 0);
prodp[j] = cy_limb;
prodp++;
continue;
do
{
umul_ppmm (prod_high, prod_low, up[j], v_limb);
add_ssaaaa (cy_limb, prod_low, prod_high, prod_low, 0, cy_limb);
c = prodp[j];
prod_low += c + 1;
prodp[j] = prod_low;
if (prod_low > c)
goto ncy_loop;
cy_loop:
j++;
}
while (j < 0);
cy_limb += 1;
prodp[j] = 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 ((size_t)((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; */
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 ((size_t)(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 ((size_t)((2 * n + 1) * BYTES_PER_MP_LIMB));
if (_mpn_cmp (up + n, usize - 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 (_mpn_cmp (vp + n, vsize - 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 ((size_t)(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 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 DEBUG
if (prod_size > 4 * n)
abort();
#endif
}
}
MPN_COPY (prodp, pp, prod_size);
return prod_size;
}
}
|
