Efficient implementation of polynomial multiplication

[Décio Luiz Gazzoni Filho]

-----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 ...
> b0   zero-pad  b1  zero-pad  b2   zero-pad ...
>
> 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-----