Efficient implementation of polynomial multiplication

Décio Luiz Gazzoni Filho decio at decpp.net
Fri Mar 5 18:13:57 CET 2004

```-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Replying to my own message...

I wrote something that seems to be working. Could anyone give me some advice
on how hopelessly non-portable this is? Here `coeffs' can be thought of an
array of mpz_classes with `deg' + 1 elements, and `arg' is an array of
mpz_classes with `arg.deg' + 1 elements. `limbs' is the number of required
limbs and `b' is the number of required bits, rounded to the nearest multiply
of `mp_bits_per_limb'.

- -------------- BEGIN CODE ---------------

mpz_t X, Y;
mpz_init2(X,(deg+1)*b);
mpz_init2(Y,(arg.deg+1*b);

// This code builds the first polynomial, whose coefficients are stored in
// coeffs[].

unsigned long j, sz;
for (unsigned long i = 0; i < deg; ++i)
{
sz = mpz_size(coeffs[i].get_mpz_t());
for (j = 0; j < sz; ++j)
X->_mp_d[i*limbs+j] = mpz_getlimbn(coeffs[i].get_mpz_t(),j);
for (; j < limbs; ++j)
X->_mp_d[i*limbs+j] = 0;
}
sz = mpz_size(coeffs[deg].get_mpz_t());
for (j = 0; j < sz; ++j)
X->_mp_d[deg*limbs+j] = mpz_getlimbn(coeffs[deg].get_mpz_t(),j);

X->_mp_size = deg*limbs+j;

// This code build the second polynomial, whose coefficients are stored in
// arg[].

for (unsigned long i = 0; i < arg.deg; ++i)
{
sz = mpz_size(arg[i].get_mpz_t());
for (j = 0; j < sz; ++j)
Y->_mp_d[i*limbs+j] = mpz_getlimbn(arg[i].get_mpz_t(),j);
for (; j < limbs; ++j)
Y->_mp_d[i*limbs+j] = 0;
}
sz = mpz_size(arg[arg.deg].get_mpz_t());
for (j = 0; j < sz; ++j)
Y->_mp_d[arg.deg*limbs+j] = mpz_getlimbn(arg[arg.deg].get_mpz_t(),j);

Y->_mp_size = arg.deg*limbs+j;

// Multiply.
mpz_t Z;
mpz_init(Z);
mpz_mul(Z,X,Y);
mpz_clear(X);
mpz_clear(Y);

// Now remove the coefficients from Z into coeffs, which stores the result.

// After the changedegree() call, coeffs can store deg + arg.deg + 1
// elements. deg is also updated to be deg + arg.deg.
changedegree(deg + arg.deg);
mp_limb_t* old_mp_d = Z->_mp_d;
unsigned long old_mp_size = Z->_mp_size;
for (unsigned long i = 0; i <= deg; ++i)
{
mpz_fdiv_r_2exp(coeffs[i].get_mpz_t(),Z,b);
Z->_mp_d += limbs;
Z->_mp_size -= limbs;
coeffs[i] %= modulo;
}

Z->_mp_d = old_mp_d;
Z->_mp_size = old_mp_size;

mpz_clear(Z);

- ----------------- END CODE ---------------

Thanks

Décio

On Friday 05 March 2004 13:12, Décio Luiz Gazzoni Filho wrote:
> Hello,
>
> I'm working on a project using GMP, and I'm stuck trying to make a fast
> routine for multiplying polynomials (with coefficients from Z/pZ, i.e. the
> size of the coefficients is bounded.) The technique I'm used is called
> `binary segmentation' and consists of evaluating the polynomials f(x), g(x)
> at x = 2^b, with b > p^2*max(deg(f),deg(g)) so that there's no risk of
> carry from one polynomial coefficient to the other. So here's how the
> resulting integers f(2^n) = a0 + a1*2^n + a2*2^(2n) + ... and g(2^n) = b0 =
> b1*2^n + b2*2^(2n) + ... look like (sorry for the ASCII art):
>
>
> 0              n             2n
> a0   00...00   a1  00...00   a2   00...00 ...
>
> So my question, what is the most efficient way of building this integer
> (and later removing coefficients from it at the same 0,n,2n,...
> boundaries), while still being portable? Of course I could use lshifts and
> adds to insert the coefficients, and mpz_fdiv_r_2exp and rshifts to remove
> them, but this is hopelessly slow: my routine spends the bulk of its time
> doing the shifts.
>
> I'm thinking either something that manipulates the mpz_t->_mp_d array
> directly, or something using mpz_import/mpz_export... anyway, I hope
> someone on the list can help me with this.
>
> I was going to put my current code here, but as it stands it's not
> generating meaningful results.
>
> Thanks for the help
>
> Décio
> _______________________________________________
> gmp-discuss mailing list
> gmp-discuss at swox.com
> https://gmplib.org/mailman/listinfo/gmp-discuss
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.3 (GNU/Linux)

iD8DBQFASO2bFXvAfvngkOIRAps4AJ99RMNKjs2cY1P+SjQvsZJ64UyNmQCgiN4m
lf8+4sPUTyoCFUD3nIdRebI=
=f6b5
-----END PGP SIGNATURE-----
```