Reduced storage for toom32
Niels Möller
nisse at lysator.liu.se
Mon Nov 16 23:05:05 CET 2009
The following version of toom32_mul uses 2n + 2 limbs of scratch
space, rather than 4n + 2 (both figures not including any scratch
needed by the pointwise multiplies).
Almost a drop in replacement, but requires s + t >= n + 2, which is
slightly stricter than the old code.
The allocation is quite clever, if I may say so. It has the drawback
of some more function call/loop setup overhead since some adds or subs
on 2n limbs have to be done n limbs at time due to the allocation layout.
Regards,
/Niels
/* mpn_toom32_mul -- Multiply {ap,an} and {bp,bn} where an is nominally 1.5
times as large as bn. Or more accurately, bn < an < 3bn.
Contributed to the GNU project by Torbjorn Granlund.
Improvements by Marco Bodrato and Niels Möller.
The idea of applying toom to unbalanced multiplication is due to Marco
Bodrato and Alberto Zanoni.
THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
Copyright 2006, 2007, 2008, 2009 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 Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, 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 Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
#include "gmp.h"
#include "gmp-impl.h"
/* Evaluate in: -1, 0, +1, +inf
<-s-><--n--><--n-->
___ ______ ______
|a2_|___a1_|___a0_|
|_b1_|___b0_|
<-t--><--n-->
v0 = a0 * b0 # A(0)*B(0)
v1 = (a0+ a1+ a2)*(b0+ b1) # A(1)*B(1) ah <= 2 bh <= 1
vm1 = (a0- a1+ a2)*(b0- b1) # A(-1)*B(-1) |ah| <= 1 bh = 0
vinf= a2 * b1 # A(inf)*B(inf)
*/
#define TOOM32_MUL_N_REC(p, a, b, n, ws) \
do { \
mpn_mul_n (p, a, b, n); \
} while (0)
void
mpn_toom32_mul (mp_ptr pp,
mp_srcptr ap, mp_size_t an,
mp_srcptr bp, mp_size_t bn,
mp_ptr scratch)
{
mp_size_t n, s, t;
mp_limb_t cy;
int vm1_neg;
int hi;
#define a0 ap
#define a1 (ap + n)
#define a2 (ap + 2 * n)
#define b0 bp
#define b1 (bp + n)
/* Required, to guarantee that s + t >= n + 2.
Also implies that bn >= 9
FIXME: Too strict? E.g., an = 15, bn = 9, ==> n = 5, s = 5, t = 4
seem to be ok. */
ASSERT (bn + 4 <= an && an + 14 <= 3*bn);
n = 1 + (2 * an >= 3 * bn ? (an - 1) / (size_t) 3 : (bn - 1) >> 1);
s = an - 2 * n;
t = bn - n;
ASSERT (0 < s && s <= n);
ASSERT (0 < t && t <= n);
ASSERT (s + t >= n + 2);
/* Product area of size an + bn = 3*n + s + t >= 4*n + 2. */
#define ap1 (pp) /* n + 1 */
#define bp1 (pp + n + 1) /* n + 1 */
#define am1 (pp + 2*n + 2) /* n, most significant bit stored elsewhere */
#define bm1 (pp + 3*n + 2) /* n */
#define v1 (scratch) /* 2n + 1 (+1) */
#define vm1 (pp) /* 2n + 1 */
#define scratch_out (scratch + 2*n + 2)
/* Scratch need: 2*n + 2 + scratch for the recursive multiplications. */
/* Compute ap1 = a0 + a1 + a3, am1 = a0 - a1 + a3 */
ap1[n] = mpn_add (ap1, a0, n, a2, s);
#if HAVE_NATIVE_mpn_addsub_n
if (ap1[n] == 0 && mpn_cmp (ap1, a1, n) < 0)
{
ap1[n] = mpn_addsub_n (ap1, am1, a1, ap1, n) >> 1;
hi = 0;
vm1_neg = 1;
}
else
{
cy = mpn_addsub_n (ap1, am1, ap1, a1, n);
hi = ap1[n] - (cy & 1);
ap1[n] += (cy >> 1);
vm1_neg = 0;
}
#else
if (ap1[n] == 0 && mpn_cmp (ap1, a1, n) < 0)
{
ASSERT_NOCARRY (mpn_sub_n (am1, a1, ap1, n));
hi = 0;
vm1_neg = 1;
}
else
{
hi = ap1[n] - mpn_sub_n (am1, ap1, a1, n);
vm1_neg = 0;
}
ap1[n] += mpn_add_n (ap1, ap1, a1, n);
#endif
/* Compute bp1 = b0 + b1 and bm1 = b0 - b1. */
if (t == n)
{
#if HAVE_NATIVE_mpn_addsub_n
if (mpn_cmp (b0, b1, n) < 0)
{
cy = mpn_addsub_n (bp1, bm1, b1, b0, n);
vm1_neg ^= 1;
}
else
{
cy = mpn_addsub_n (bp1, bm1, b0, b1, n);
}
bp1[n] = cy >> 1;
#else
bp1[n] = mpn_add_n (bp1, b0, b1, n);
if (mpn_cmp (b0, b1, n) < 0)
{
ASSERT_NOCARRY (mpn_sub_n (bm1, b1, b0, n));
vm1_neg ^= 1;
}
else
{
ASSERT_NOCARRY (mpn_sub_n (bm1, b0, b1, n));
}
#endif
}
else
{
/* FIXME: Should still use addsub for the main part. */
bp1[n] = mpn_add (bp1, b0, n, b1, t);
if (mpn_zero_p (b0 + t, n - t) && mpn_cmp (b0, b1, t) < 0)
{
ASSERT_NOCARRY (mpn_sub_n (bm1, b1, b0, t));
MPN_ZERO (bm1 + t, n - t);
vm1_neg ^= 1;
}
else
{
ASSERT_NOCARRY (mpn_sub (bm1, b0, n, b1, t));
}
}
TOOM32_MUL_N_REC (v1, ap1, bp1, n+1, scratch_out);
TOOM32_MUL_N_REC (vm1, am1, bm1, n, scratch_out);
if (hi)
hi = mpn_add_n (vm1+n, vm1+n, bm1, n);
vm1[2*n] = hi;
/* v1 <-- (v1 + vm1) / 2 = x0 + x2 */
if (vm1_neg)
{
#if HAVE_NATIVE_mpn_rsh1sub_n
mpn_rsh1sub_n (v1, v1, vm1, 2*n+1);
#else
mpn_sub_n (v1, v1, vm1, 2*n+1);
ASSERT_NOCARRY (mpn_rshift (v1, v1, 2*n+1, 1));
#endif
}
else
{
#if HAVE_NATIVE_mpn_rsh1add_n
mpn_rsh1add_n (v1, v1, vm1, 2*n+1);
#else
mpn_add_n (v1, v1, vm1, 2*n+1);
ASSERT_NOCARRY (mpn_rshift (v1, v1, 2*n+1, 1));
#endif
}
/* We get x1 + x3 = (x0 + x2) - (x0 - x1 + x2 - x3), and hence
y = x1 + x3 + (x0 + x2) * B
= (x0 + x2) * B + (x0 + x2) - vm1.
y is 3*n + 1 limbs, y = y0 + y1 B + y2 B^2. We store them as
follows: y0 at scratch, y1 at pp + 2*n, and y2 at scratch + n
(already in place, except for carry propagation).
We thus add
B^3 B^2 B 1
| | | |
+-----+----+
+ | x0 + x2 |
+----+-----+----+
+ | x0 + x2 |
+----------+
- | vm1 |
--+----++----+----+-
| y2 | y1 | y0 |
+-----+----+----+
Since we store y0 at the same location as the low half of x0 + x2, we
need to do the middle sum first. */
hi = vm1[2*n];
cy = mpn_add_n (pp + 2*n, v1, v1 + n, n);
MPN_INCR_U (v1 + n, n + 1, cy + v1[2*n]);
if (vm1_neg)
{
cy = mpn_add_n (v1, v1, vm1, n);
hi += mpn_add_nc (pp + 2*n, pp + 2*n, vm1 + n, n, cy);
MPN_INCR_U (v1 + n, n+1, hi);
}
else
{
cy = mpn_sub_n (v1, v1, vm1, n);
hi += mpn_sub_nc (pp + 2*n, pp + 2*n, vm1 + n, n, cy);
MPN_DECR_U (v1 + n, n+1, hi);
}
TOOM32_MUL_N_REC (pp, a0, b0, n, scratch_out);
/* vinf, s+t limbs. Use mpn_mul for now, to handle unbalanced operands */
if (s > t) mpn_mul (pp+3*n, a2, s, b1, t);
else mpn_mul (pp+3*n, b1, t, a2, s);
/* Remaining interpolation.
y * B + x0 + x3 B^3 - x0 B^2 - x3 B
= (x1 + x3) B + (x0 + x2) B^2 + x0 + x3 B^3 - x0 B^2 - x3 B
= y0 B + y1 B^2 + y3 B^3 + Lx0 + H x0 B
+ L x3 B^3 + H x3 B^4 - Lx0 B^2 - H x0 B^3 - L x3 B - H x3 B^2
= L x0 + (y0 + H x0 - L x3) B + (y1 - L x0 - H x3) B^2
+ (y2 - (H x0 - L x3)) B^3 + H x3 B^4
B^4 B^3 B^2 B 1
| | | | | |
+-------+ +---------+---------+
| Hx3 | | Hx0-Lx3 | Lx0 |
+------+----------+---------+---------+---------+
| y2 | y1 | y0 |
++---------+---------+---------+
-| Hx0-Lx3 | - Lx0 |
+---------+---------+
| - Hx3 |
+--------+
We must take into account the carry from Hx0 - Lx3.
*/
cy = mpn_sub_n (pp + n, pp + n, pp+3*n, n);
hi = scratch[2*n] + cy;
cy = mpn_sub_nc (pp + 2*n, pp + 2*n, pp, n, cy);
hi -= mpn_sub_nc (pp + 3*n, scratch + n, pp + n, n, cy);
hi += mpn_add (pp + n, pp + n, 3*n, scratch, n);
hi -= mpn_sub (pp + 2*n, pp + 2*n, 2*n, pp + 4*n, s+t-n);
if (hi < 0)
MPN_DECR_U (pp + 4*n, s+t-n, -hi);
else
MPN_INCR_U (pp + 4*n, s+t-n, hi);
}
More information about the gmp-devel
mailing list