/Users/mourrain/Devel/mmx/algebramix/src/kronecker_integer.cpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : kronecker_integer.cpp
* DESCRIPTION: Breaking big integers into parts and recomposing them
* COPYRIGHT  : (C) 2007  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <numerix/integer.hpp>
namespace mmx {

/******************************************************************************
* Naive encoding and decoding of integer polynomials using Kronecker's method
******************************************************************************/

#if !defined(__GNU_MP__)

void
encode_kronecker (integer& dest, const integer* src, nat n, xnat bits) {
  if (n == 0) dest= 0;
  else if (n == 1) dest= src[0];
  else {
    nat h= n>>1;
    integer aux;
    encode_kronecker (aux , src+h, n-h, bits);
    encode_kronecker (dest, src  , h  , bits);
    dest += (aux << (h*bits));
  }
}

void
decode_kronecker (integer* dest, nat n, xnat bits, const integer& src) {
  if (n == 0);
  else if (n == 1) dest[0]= src;
  else if (src > 0) {
    nat h= n>>1;
    integer aux= src >> (h*bits);
    if (src[h*bits-1]) aux += 1;
    decode_kronecker (dest+h, n-h, bits, aux);
    aux= src - (aux << (h*bits));
    decode_kronecker (dest, h, bits, aux);
  }
  else {
    integer bis= -src;
    nat h= n>>1;
    integer aux= bis >> (h*bits);
    if (bis[h*bits-1]) aux += 1;
    decode_kronecker (dest+h, n-h, bits, -aux);
    aux= bis - (aux << (h*bits));
    decode_kronecker (dest, h, bits, -aux);
  }
}

#else
/******************************************************************************
* Specific implementation for GMP
******************************************************************************/

// from polynomial_gmp.cpp
void mpz_encode_kronecker
       (mpz_t dest, const mpz_t* src, unsigned long n, unsigned long bits);
void mpz_decode_kronecker
       (mpz_t* dest, unsigned long n, unsigned long bits, const mpz_t src_bis);

static inline mpz_t*
mmx_new_mpz (nat n) {
  return (mpz_t*) mmx_malloc (n * sizeof (mpz_t));
}

static inline void
mmx_delete_mpz (mpz_t* Ptr, nat n) {
  mmx_free ((void*) Ptr, n * sizeof (mpz_t));
}

void
encode_kronecker (integer& dest, const integer* src, nat n, xnat bits) {
  nat i;
  mpz_t* Src= mmx_new_mpz (n);
  for (i=0; i<n; i++) {
    const mpz_t& aux= *(src[i]);
    Src[i]->_mp_alloc= aux->_mp_alloc;
    Src[i]->_mp_size = aux->_mp_size;
    Src[i]->_mp_d    = aux->_mp_d;
  }
  mpz_encode_kronecker (*dest, Src, n, bits);
  mmx_delete_mpz (Src, n);
}

void
decode_kronecker (integer* dest, nat n, xnat bits, const integer& src) {
  nat i;
  mpz_t* Dest= mmx_new_mpz (n);
  for (i=0; i<n; i++) mpz_init (Dest[i]);
  mpz_decode_kronecker (Dest, n, bits, *src);
  for (i=0; i<n; i++) {
    mpz_t& aux= *(dest[i]);
    aux->_mp_alloc= Dest[i]->_mp_alloc;
    aux->_mp_size = Dest[i]->_mp_size;
    aux->_mp_d    = Dest[i]->_mp_d;
  }
  mmx_delete_mpz (Dest, n);
}
#endif // __GNU_MP__

/******************************************************************************
* Multiply two integer polynomials stored in segments using Kronecker's method
******************************************************************************/

xnat max_bit_size (const integer* src, nat n) {
  xnat m= 0;
  for (nat i=0; i<n; i++)
    m= max (m, bit_size (src[i]));
  return m;
}

void
mul_kronecker (integer* dest,
               const integer* src1, nat n1,
               const integer* src2, nat n2)
{
  /*
  mmout << "n1= " << n1 << ", " << "n2= " << n2 << "\n";
  for (nat i=0; i<n1; i++)
    mmout << "  s" << i << "\t" << src1[i] << "\n";
  for (nat i=0; i<n2; i++)
    mmout << "  t" << i << "\t" << src2[i] << "\n";
  */
  if (n1 == 0 && n2 == 0) return;
  for (nat i= 0; i < n1 + n2 - 1; i++) dest[i]= 0;
  while (n1 > 0 && src1[n1-1] == 0) n1--;
  while (n2 > 0 && src2[n2-1] == 0) n2--;
  if (n1 == 0 || n2 == 0) return;
  
  xnat bits1= max_bit_size (src1, n1);
  xnat bits2= max_bit_size (src2, n2);
  xnat bits= bits1 + bits2 + bit_size (integer (min (n1, n2))) + 1;
  integer aux1, aux2;
  //mmout << "Encoding\n";
  //nat start= mmx_time ();
  encode_kronecker (aux1, src1, n1, bits);
  encode_kronecker (aux2, src2, n2, bits);
  //mmout << "Done in " << mmx_time () - start << "ms\n";
  //mmout << "Multiplying\n";
  //start= mmx_time ();
  integer aux= aux1*aux2;
  //mmout << "Done in " << mmx_time () - start << "ms\n";
  //mmout << "Decoding\n";
  //start= mmx_time ();
  decode_kronecker (dest, n1+n2-1, bits, aux);
  //mmout << "Done in " << mmx_time () - start << "ms\n";

  /*
  integer* dest1= mmx_new<integer> (n1+n2-1);
  decode_kronecker_naive (dest1, n1+n2-1, bits, aux);
  decode_kronecker (dest, n1+n2-1, bits, aux);
  for (nat i=0; i<n1+n2-1; i++) {
    mmout << "naive\t" << i << "\t" << dest1[i] << "\n";
    mmout << "fast\t" << i << "\t" << dest[i] << "\n";
  }
  mmx_delete<integer> (dest1, n1+n2-1);
  */
}

void
square_kronecker (integer* dest, const integer* src, nat n) {
  if (n == 0) return;
  for (nat i= 0; i < 2 * n - 1; i++) dest[i]= 0;
  while (n > 0 && src[n-1] == 0) n--;
  if (n == 0) return;
  if (n == 1) { dest[0]= square (src[0]); return; }

  xnat bits1= max_bit_size (src, n);
  xnat bits = (bits1 << 1) + bit_size (integer (n)) + 1;

  integer aux1;
  encode_kronecker (aux1, src, n, bits);
  integer aux= square (aux1);
  decode_kronecker (dest, n+n-1, bits, aux);
}

} // namespace mmx