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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : kronecker_gmp.cpp
* DESCRIPTION: Breaking big integers into parts and recomposing them
* COPYRIGHT  : (C) 2007  Joris van der Hoeven and Gregoire Lecerf
*******************************************************************************
* 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 {
#if defined(__GNU_MP__)

#define ASSERT_RANGE(cond) assert(cond)
//#define ASSERT_RANGE(cond) VERIFY(cond,"limb out of range")

static inline void
mpn_set (mp_limb_t* dest, const mp_limb_t* src, mp_size_t size) {
  mp_size_t i;
  for (i=0; i<size; i++) dest[i]= src[i];
}

template<typename I> static inline void
mpn_get2 (I& dest, const mp_limb_t* src, unsigned long n,
          unsigned long beg, unsigned long end) {
  // write bits of src from index beg to end-1 into dest
  typedef typename unsigned_of_helper<I>::type U;
  dest= 0; 
  if (end <= beg) return;
  const xnat shift_l= bit_size (BITS_PER_LIMB) - 1;
  const mp_limb_t ceil_l= (- (mp_limb_t) 1) << shift_l;
  const mp_limb_t mod_l= (- (mp_limb_t) 1) - ceil_l;
  end= min (beg + 8 * sizeof (I), min (end, BITS_PER_LIMB * n));
  unsigned long pos= beg;
  while (pos < end) {
    unsigned long i= pos >> shift_l, o= pos & mod_l;
    dest |= ((I) (src[i] >> o)) << (pos - beg);
    pos += BITS_PER_LIMB - o;
  }
  dest &= (((U) -1) >> (8 * sizeof (I) - (end - beg)));
}

static inline void
mpn_get2 (mp_limb_t* dest, const mp_limb_t* src, unsigned long n,
          unsigned long beg, unsigned long end) {
  // write the bits from beg to end-1 of src into dest
  unsigned long i= 0;
  while (beg < end) {
    mpn_get2 (dest[i], src, n, beg, min (end, beg + BITS_PER_LIMB));
    i++; beg += BITS_PER_LIMB;
  }
}

/******************************************************************************
* Encoding a polynomial of integers using Kronecker's method
******************************************************************************/

static inline void
mpn_lshift_bis (mp_limb_t* dest, mp_size_t dest_size,
                const mp_limb_t* src , mp_size_t src_size,
                unsigned int count)
{
  assert (dest_size >= src_size && dest_size > 0);
  if (dest_size == src_size) {
    if (count == 0) mpn_set (dest, src, src_size);
    else (void) mpn_lshift (dest, src, src_size, count);
  }
  else {
    mp_size_t i;
    if (src_size == 0) dest[src_size]= 0;
    else {
      if (count == 0) {
        mpn_set (dest, src, src_size);
        dest[src_size]= 0;
      }
      else dest[src_size]= mpn_lshift (dest, src, src_size, count);
    }
    for (i=src_size+1; i<dest_size; i++) dest[i]= 0;
  }
}

void
mpz_encode_kronecker (mpz_t dest,
                      const mpz_t* src, unsigned long n, unsigned long bits)
{
  // src is an array of length n of integers which fit into bits-1 bits.
  // We set dest := src[0]+ src[1] 2^{bits} + ... + src[n-1] 2^{(n-1) bits}.

  unsigned long i;
  mp_size_t dest_size;
  mp_size_t cleaned;
  mpz_t aux;
  int dest_sgn;

  while (n>0 && mpz_sgn (src[n-1]) == 0) n--;
  mpz_set_si (dest, 0);
  if (n == 0) return;
  dest_sgn= mpz_sgn (src[n-1]);

  mpz_realloc2 (dest, n * bits);
  dest_size  = dest->_mp_alloc;
  cleaned= dest_size;
  //mmout << "dest_size= " << dest_size << "\n";

  mpz_init2 (aux, bits + BITS_PER_LIMB);
  for (i=n-1; i != ((unsigned long) -1); i--) {
    mp_size_t    j;
    mp_size_t    bit     = i * bits;
    unsigned int shift   = bit & (BITS_PER_LIMB-1);
    mp_size_t    offset  = bit / BITS_PER_LIMB;
    mp_size_t    extra   = cleaned - offset;
    mp_size_t    aux_size= aux->_mp_alloc;
    mp_size_t    src_ssz = src[i]->_mp_size;
    mp_size_t    src_size= (src_ssz>0? src_ssz: -src_ssz);
    int          src_sgn = (src_ssz>=0? 1: -1);
    if (aux_size + offset > dest_size) aux_size--;
    /*
    mmout << "Step " << i << "\n";
    mmout << "  bit     = " << bit << "\n";
    mmout << "  shift   = " << shift << "\n";
    mmout << "  offset  = " << offset << "\n";
    mmout << "  extra   = " << extra << "\n";
    mmout << "  aux_size= " << aux_size << "\n";
    mmout << "  src_size= " << src_size << "\n";
    mmout << "  dest_sgn= " << dest_sgn << "\n";
    mmout << "  src_sgn = " << src_sgn << "\n";
    */

    // ASSERT_RANGE (aux->_mp_alloc >= aux_size); // Useless, GL.
    mpn_lshift_bis (aux->_mp_d, aux_size, src[i]->_mp_d, src_size, shift);
    if (src_sgn == dest_sgn) {
      for (j=0; j<extra; j++)
        dest->_mp_d[j+offset]= aux->_mp_d[j];
      if (aux_size > extra) {
        ASSERT_RANGE (dest->_mp_alloc >= cleaned + aux_size - extra);
        mpn_add_n (dest->_mp_d + cleaned, dest->_mp_d + cleaned,
                   aux->_mp_d + extra, aux_size - extra);
      }
    }
    else {
      ASSERT_RANGE (dest->_mp_alloc >= offset + extra);
      for (j=0; j<extra; j++)
        dest->_mp_d[j+offset]= 0;
      ASSERT_RANGE (dest->_mp_alloc >= dest_size);
      mpn_sub (dest->_mp_d + offset,
               dest->_mp_d + offset, dest_size - offset,
               aux->_mp_d, aux_size);
    }
    cleaned= offset;
  }
  mpz_clear (aux);

  while (dest->_mp_d[dest_size-1] == 0) dest_size--;
  dest->_mp_size= (dest_sgn>0? dest_size: -dest_size);
}

template<typename I> static inline void
mpz_encode_kronecker_uint (mpz_t dest,
                           const I* src, unsigned long n, unsigned long bits)
{
  // src is an array of length n of unsigned hardware integers which fit
  // into e bits. We set dest := src[0]+ src[1] 2^{bits} + ... +
  // src[n-1] 2^{(n-1) bits}.
  while (n > 0 && src[n-1] == 0) n--;
  mpz_set_si (dest, 0);
  if (n == 0) return;
  static const xnat k= 8 * sizeof (I);
  const xnat shift_k= bit_size (k) - 1;
  const unsigned long ceil_k= (- (unsigned long) 1) << shift_k;
  const unsigned long mod_k= (- (unsigned long) 1) - ceil_k;
  mpz_realloc2 (dest, (n * bits + k - 1) & ceil_k);
  for (unsigned long i= 0; i < (unsigned long) dest->_mp_alloc; i++)
    dest->_mp_d[i]= 0;
  unsigned long pos= 0, j, o;
  I* aux= (I*) (void*) dest->_mp_d;
  for (unsigned long i= 0; i < n; i++, pos += bits) {
    j= pos >> shift_k, o= pos & mod_k; 
                aux[j    ] |= src[i] << o;
    if (o != 0) aux[j + 1] |= src[i] >> (k - o);
  }
  dest->_mp_size= dest->_mp_alloc;
  while (dest->_mp_size > 0 && dest->_mp_d[dest->_mp_size - 1] == 0)
    dest->_mp_size--;
}

#define IMPLEMENTATION_HELPER(I)                                        \
  void mpz_encode_kronecker (mpz_t dest,                                \
                             const I* src, unsigned long n,             \
                             unsigned long bits) {                      \
    mpz_encode_kronecker_uint (dest, src, n, bits); }

IMPLEMENTATION_HELPER(unsigned char)
IMPLEMENTATION_HELPER(unsigned short)
IMPLEMENTATION_HELPER(unsigned int)
IMPLEMENTATION_HELPER(unsigned long int)
IMPLEMENTATION_HELPER(unsigned long long int)
#undef IMPLEMENTATION_HELPER

/******************************************************************************
* Decoding a polynomial of integers using Kronecker's method
******************************************************************************/

static inline void
mpn_rshift_bis (mp_limb_t* dest, mp_size_t dest_size,
                const mp_limb_t* src , mp_size_t src_size,
                unsigned int count)
{
  assert (dest_size > 0);
  if (src_size >= dest_size) {
    if (count == 0) mpn_set (dest, src, dest_size);
    else (void) mpn_rshift (dest, src, dest_size, count);
  }
  else {
    mp_size_t i;
    if (src_size != 0) {
      if (count == 0) mpn_set (dest, src, src_size);
      else mpn_rshift (dest, src, src_size, count);
    }
    for (i=src_size; i<dest_size; i++) dest[i]= 0;
  }
}

void
mpz_decode_kronecker (mpz_t* dest, unsigned long n, unsigned long bits,
                      const mpz_t src_bis)
{
  // let src be an integer of the form a_0 + ... + a_{n-1} 2^{(n-1) bits},
  // where the a_i are positive or negative integers of at most bits-1 bits.
  // On exit, we have dest[i] := a_i for each i.

  unsigned long i, src_bit_size;
  mp_size_t src_size, aux_size;
  int src_sgn;
  mpz_t src, aux, tmp;
  mp_limb_t mask= (((mp_limb_t) 1) << (bits & (BITS_PER_LIMB-1))) - 1;
  if (mask == 0) mask= (mp_limb_t) -1;
  mpz_init2 (src, (mpz_size (src_bis) + 1) * BITS_PER_LIMB);
  mpz_set (src, src_bis);
  ASSERT_RANGE (((long int) src->_mp_alloc) > ((long int) mpz_size (src)));
  src->_mp_d[mpz_size (src)]= 0;

  /* handle trailing zeros */
  if (n == 0) return;
  src_bit_size= mpz_sizeinbase (src, 2);
  if (src_bit_size < (n-1) * bits) {
    unsigned long new_n= (src_bit_size / bits) + 1;
    for (i= new_n; i<n; i++) mpz_set_si (dest[i], 0);
    n= new_n;
  }
  ASSERT_RANGE (src_bit_size < n*bits);
  ASSERT_RANGE (src_bit_size >= (n-1) * bits);

  src_size= mpz_size (src) + 1;
  src_sgn= mpz_sgn (src);
  mpz_init2 (aux, bits + BITS_PER_LIMB);
  aux_size= aux->_mp_alloc;
  mpz_init2 (tmp, bits + BITS_PER_LIMB);
  mpz_set_si (tmp, 1);
  mpz_mul_2exp (tmp, tmp, bits);

  /*
  mmout << "mask\t= " << mask << "\n";
  mmout << "n\t= " << n << "\n";
  mmout << "src_sz\t= " << src_size << "\n";
  mmout << "src_sgn\t= " << src_sgn << "\n";
  mmout << "aux_sz\t= " << aux_size << "\n";
  */

  for (i=0; i<n; i++) {
    int          dest_sgn = 1;
    mp_size_t    bit      = i * bits;
    unsigned int shift    = bit & (BITS_PER_LIMB-1);
    mp_size_t    offset   = bit / BITS_PER_LIMB;

    /*
    mmout << "Stage " << i << "\n";
    mmout << "  bit\t= " << bit << "\n";
    mmout << "  total\t= " << src_bit_size << "\n";
    mmout << "  shift\t= " << shift << "\n";
    mmout << "  off\t= " << offset << "\n";
    */

    /* put the raw next coefficient in aux */
    ASSERT_RANGE (aux->_mp_alloc >= aux_size);
    mpn_rshift_bis (aux->_mp_d, aux_size,
                    src->_mp_d+offset, src_size-offset, shift);
    ASSERT_RANGE (aux->_mp_alloc >= aux_size);
    aux->_mp_d[aux_size-2] &= mask;
    aux->_mp_size= aux_size - 1;
    while (aux->_mp_size > 0 && aux->_mp_d[aux->_mp_size-1] == 0)
      aux->_mp_size--;

    /* check whether the next coefficient has a negative sign */
    if (i != n-1) {
      unsigned int k= (bits-1) & (BITS_PER_LIMB-1);
      mp_size_t    l= (bits-1) / BITS_PER_LIMB;
      dest_sgn= ((aux->_mp_d[l] & ((mp_limb_t) 1<<k)) == 0? 1: -1);
      if (dest_sgn < 0) {
        /* in-place update of src in case of a negative sign */
        mp_size_t    next_bit   = (i+1) * bits;
        unsigned int next_shift = next_bit & (BITS_PER_LIMB-1);
        mp_size_t    next_offset= next_bit / BITS_PER_LIMB;
        mp_limb_t    plus       = ((mp_limb_t) 1) << next_shift;
        ASSERT_RANGE (src->_mp_alloc >= src_size);
        if (mpn_add (src->_mp_d + next_offset,
                     src->_mp_d + next_offset, src_size - next_offset,
                     &plus, 1) != 0)
          ASSERT (false, "unexpected carry (mpz_decode_kronecker)");
      }
      //mmout << "  sgn\t= " << dest_sgn << "\n";
    }

    /* put the sign-corrected next coefficient in dest[i] */
    if (dest_sgn >= 0) mpz_set (dest[i], aux);
    else mpz_sub (dest[i], tmp, aux);
    if (dest_sgn != src_sgn) dest[i]->_mp_size= -dest[i]->_mp_size;
  }

  mpz_clear (src);
  mpz_clear (aux);
  mpz_clear (tmp);
}

template<typename I> static inline void
mpz_decode_kronecker_uint (I* dest, unsigned long n, unsigned long bits,
                           const mpz_t src,
                           unsigned int e= 8 * sizeof (I)) {
  ASSERT (2*e <= bits, "overflow");
  const xnat k= 8 * sizeof (I);
  const xnat shift_k= bit_size (k) - 1;
  const unsigned long ceil_k= (- (unsigned long) 1) << shift_k;
  const unsigned long mod_k= (- (unsigned long) 1) - ceil_k;
  unsigned long i= 0, j, o, sz= BITS_PER_LIMB * src->_mp_size, pos= 0;
  const I* aux= (I*) (void*) src->_mp_d;
  for (; pos + (k << 1) < sz && i < n; i++, pos += bits) {
    j= pos >> shift_k, o= pos & mod_k;
    if (o == 0) 
      dest[i]= aux[j] >> o;
    else
      dest[i]= (aux[j] >> o) | (aux[j+1] << (k - o));
  }
  for (; pos < sz && i < n; i++, pos += bits)
    mpn_get2 (dest[i], src->_mp_d, src->_mp_size, pos, pos+bits);
  for (; i < n; i++) dest[i]= 0;
}
#define IMPLEMENTATION_HELPER(I)                                        \
  void mpz_decode_kronecker (I* dest, unsigned long n,                  \
                             unsigned long bits, const mpz_t src) {     \
    mpz_decode_kronecker_uint (dest, n, bits, src); }

IMPLEMENTATION_HELPER(unsigned char)
IMPLEMENTATION_HELPER(unsigned short)
IMPLEMENTATION_HELPER(unsigned int)
IMPLEMENTATION_HELPER(unsigned long int)
IMPLEMENTATION_HELPER(unsigned long long int)
#undef IMPLEMENTATION_HELPER

/******************************************************************************
* Decoding a polynomial of integers using Kronecker's method
******************************************************************************/

template<typename L,typename I> static inline void
_mpz_decode_kronecker_mod_uint (I* dest,
                                unsigned long n, unsigned long bits,
                                const mpz_t src, const I& p) {
  // Kronecker decoding modulo p in case when L is a hardware integer
  // type with at least 'bits' bits.
  const xnat k= 8 * sizeof (L);
  ASSERT (k >= bits, "overflow");
  const xnat shift_k= bit_size (k) - 1;
  const L mask= ((L) -1) >> (k - bits);
  const unsigned long ceil_k= (- (unsigned long) 1) << shift_k;
  const unsigned long mod_k= (- (unsigned long) 1) - ceil_k;
  unsigned long i= 0, j, o, sz= BITS_PER_LIMB * src->_mp_size, pos= 0;
  const L* aux= (L*) (void*) src->_mp_d;
  for (; pos + (k << 1) < sz && i < n; i++, pos += bits) {
    j= pos >> shift_k, o= pos & mod_k;
    if (o == 0) 
      dest[i]= ((aux[j] >> o) & mask) % p;
    else
      dest[i]= (((aux[j] >> o) | (aux[j+1] << (k - o))) & mask) % p;
  }
  for (; pos < sz && i < n; i++, pos += bits) {
    L tmp; mpn_get2 (tmp, src->_mp_d, src->_mp_size, pos, pos+bits);
    dest[i]= tmp % p;
  }
  for (; i < n; i++) dest[i]= 0;
}

template<typename I> static inline void
_mpz_decode_kronecker_mod_integer (I* dest, unsigned long n,
                                   unsigned long bits, const mpz_t src,
                                   const I& p)
{
  // let src be an integer of the form a_0 + ... + a_{n-1} 2^{(n-1) bits},
  // where the a_i are positive or negative integers of at most bits-1 bits.
  // On exit, we have dest[i] := a_i mod p for each i.
  ASSERT (bit_size (p) <= BITS_PER_LIMB, "overflow");
  const mp_size_t e= (bits + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
  unsigned long sz= src->_mp_size * BITS_PER_LIMB, i= 0;
  mp_limb_t* aux= mmx_new<mp_limb_t> (2*e);
  mp_limb_t* q= aux + e;
  for (unsigned long pos= 0; pos < sz && i < n; i++, pos += bits) {
    mpn_get2 (aux, src->_mp_d, src->_mp_size, pos, pos + bits);
    dest[i]= mpn_divrem_1 (q, 0, aux, e, (mp_limb_t) p);
  }
  for (; i < n; i++) dest[i]= 0;
  mmx_delete<mp_limb_t> (aux, 2*e);
}

template<typename I> static inline void
_mpz_decode_kronecker_mod_integer (I* dest, unsigned long n,
                                   unsigned long bits, const mpz_t src,
                                   const mpz_t p)
{
  // let src be an integer of the form a_0 + ... + a_{n-1} 2^{(n-1) bits},
  // where the a_i are positive or negative integers of at most bits-1 bits.
  // On exit, we have dest[i] := a_i mod p for each i.
  const mp_size_t e= (bits + BITS_PER_LIMB - 1) / BITS_PER_LIMB;
  unsigned long sz= src->_mp_size * BITS_PER_LIMB, i= 0;
  mp_limb_t* aux= mmx_new<mp_limb_t> (2*e + p->_mp_size);
  mp_limb_t* q= aux + e, * r= q + e;
  for (unsigned long pos= 0; pos < sz && i < n; i++, pos += bits) {
    mpn_get2 (aux, src->_mp_d, src->_mp_size, pos, pos + bits);
    mpn_tdiv_qr (q, r, 0, aux, e, p->_mp_d, p->_mp_size);
    dest[i]= r[0];
    for (xnat j= 1; j < (xnat) p->_mp_size; j++)
      dest[i] |= ((I) r[j]) << (j * BITS_PER_LIMB);
  }
  for (; i < n; i++) dest[i]= 0;
  mmx_delete<mp_limb_t> (aux, 2*e + p->_mp_size);
}

template<typename I> static inline void
mpz_decode_kronecker_mod_int (I* dest, unsigned long n, unsigned long bits,
                              const mpz_t src, const I& p) {
  if (p == 0) {
    typedef typename unsigned_of_helper<I>::type U;
    mpz_decode_kronecker ((U*) (void*) dest, n, bits, src);
    return;
  }
  if (bits <= 8 * sizeof (unsigned char)) {
    typedef unsigned char L;
    _mpz_decode_kronecker_mod_uint<L> (dest, n, bits, src, p);
  }
  else if (bits <= 8 * sizeof (unsigned short int)) {
    typedef unsigned short int L;
    _mpz_decode_kronecker_mod_uint<L> (dest, n, bits, src, p);
  }
  else if (bits <= 8 * sizeof (unsigned int)) {
    typedef unsigned int L;
    _mpz_decode_kronecker_mod_uint<L> (dest, n, bits, src, p);
  }
  else if (bits <= 8 * sizeof (unsigned long int)) {
    typedef unsigned long int L;
    _mpz_decode_kronecker_mod_uint<L> (dest, n, bits, src, p);
  }
  else if (bits <= 8 * sizeof (unsigned long long int)) {
    typedef unsigned long long int L;
    _mpz_decode_kronecker_mod_uint<L> (dest, n, bits, src, p);
  }
  else
    if (bit_size (p) <= BITS_PER_LIMB)
      _mpz_decode_kronecker_mod_integer (dest, n, bits, src, p);
    else {
      integer _p (p);
      _mpz_decode_kronecker_mod_integer (dest, n, bits, src, * _p); } }

#define DECLARE_HELPER(I)                                               \
  void mpz_decode_kronecker_mod (I* dest, unsigned long n,              \
                                 unsigned long bits, const mpz_t src,   \
                                 const I& p) {                          \
    mpz_decode_kronecker_mod_int (dest, n, bits, src, p); }
DECLARE_HELPER(signed char)
DECLARE_HELPER(unsigned char)
DECLARE_HELPER(short int)
DECLARE_HELPER(unsigned short)
DECLARE_HELPER(int)
DECLARE_HELPER(unsigned int)
DECLARE_HELPER(long int)
DECLARE_HELPER(unsigned long int)
DECLARE_HELPER(long long int)
DECLARE_HELPER(unsigned long long int)
#undef DECLARE_HELPER

#endif // __GNU_MP__
} // namespace mmx