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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : base_double.cpp
* DESCRIPTION: Cutting integers into pieces for FFT-multiplication
* COPYRIGHT  : (C) 2008  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.
******************************************************************************/

#ifdef BASIX_HAVE_STDINT_H
#include <stdint.h>
typedef uint64_t lnat;
#else
typedef unsigned long long int lnat;
#endif
#include <numerix/integer.hpp>
#include <algebramix/fft_roots.hpp>

namespace mmx {
#if defined(__GNU_MP__)

/******************************************************************************
* Encoding integers as vectors of doubles of a fixed bitsize
******************************************************************************/

static inline nat
mpn_size (const mp_limb_t* src, nat n) {
  while (n > 0 && src[n-1] == 0) n--;
  return n;
}

void
mpz_base_encode2 (double* dest, nat n,
                   const mpz_t& src, unsigned long cap) {
  // dest: array of n nats
  // src: integer to be cut in n pieces of cap bits.
  VERIFY (cap <= BITS_PER_LIMB, "invalid capacity");
  static nat lbpl= log_2 (BITS_PER_LIMB);
  mp_limb_t* sd= src->_mp_d;
  nat sn= abs (src->_mp_size);
  nat i;
  unsigned long pos;
  unsigned long end= sn << lbpl;
  mp_limb_t mask= (((mp_limb_t) 1) << cap) - 1;
  double sgn= mpz_sgn (src);
  //mmout << "sign= " << sgn << "\n";
  //for (nat i=0; i<sn; i++)
  //mmout << "source " << i << "= " << sd[i] << "\n";
  for (i=0, pos=0; pos < end; i++, pos += cap) {
    nat index  = pos >> lbpl;
    nat offset = pos & (BITS_PER_LIMB - 1);
    mp_limb_t x= sd[index] >> offset;
    if (offset + cap > BITS_PER_LIMB && pos + cap < end && offset != 0)
      x += sd[index+1] << (BITS_PER_LIMB - offset);
    dest[i]= sgn * ((double) (x & mask));
  }
  for (; i<n; i++) dest[i]= 0;
}

void
mpz_base_decode2 (mpz_t& dest, const double* src, nat n, unsigned long cap) {
  VERIFY (cap <= BITS_PER_LIMB, "invalid capacity");
  nat i;
  double sgn= 0;
  for (i= n-1; i!=((nat) (-1)); i--)
    if (abs (src[i]) > 0.5) break;
  if (i == ((nat) (-1))) {
    mpz_set_si (dest, 0);
    return;
  }
  sgn= sign (src[i]);
  static nat lbpl= log_2 (BITS_PER_LIMB);
  static mp_limb_t hibit= ((mp_limb_t) 1) << (BITS_PER_LIMB - 1);
  mp_limb_t* dd= dest->_mp_d;
  nat dn= dest->_mp_alloc;
  for (nat i=0; i<dn; i++) dd[i]= 0;
  unsigned long pos;
  nat extra= (BITS_PER_LIMB == 32? 3: 2);
  VERIFY (extra == 3 || BITS_PER_LIMB >= 64, "unsupported architecture");
  ASSERT ((dn + 1 - extra) * BITS_PER_LIMB >= n * cap, "oveflow");
  //ASSERT ((dn - extra) * BITS_PER_LIMB >= n * cap, "oveflow");
  for (i=0, pos=0; i<n; i++, pos += cap) {
    nat index = pos >> lbpl;
    nat offset= pos & (BITS_PER_LIMB - 1);
    double f= floor (sgn * src[i] + 0.5);
    lnat   x= (lnat) abs (f);
    mp_limb_t pp[extra];
    pp[0]= (mp_limb_t) x;
    if (BITS_PER_LIMB == 32) {
      pp[1]= (mp_limb_t) ((x >> 16) >> 16);
      pp[2]= 0;
    }
    else pp[1]= 0;
    if (offset != 0) mpn_lshift (pp, pp, extra, offset);
    //mmout << "Stage " << i << "\n";
    //mmout << "  B= " << BITS_PER_LIMB << "\n";
    //mmout << "  I= " << index << "\n";
    //mmout << "  O= " << offset << "\n";
    //mmout << "  f= " << f << "\n";
    //mmout << "  x= " << x << "\n";
    //mmout << "  p= " << pp[0] << ", " << pp[1] << ", " << pp[2] << "\n";
    if (f > 0) mpn_add_n (dd + index, dd + index, pp, extra);
    else mpn_sub_n (dd + index, dd + index, pp, extra);
    //mmout << "  r= " << dd[index] << ", " << dd[index+1] << ", " << dd[index+2] << "\n";
    nat new_index= (pos + cap) >> lbpl;
    if (new_index > index && ((dd[index+extra-1] & hibit) != 0))
      dd[index+extra]= (mp_limb_t) (-1);
  }
  //for (nat i=0; i<dn; i++)
  //mmout << "result " << i << "= " << dd[i] << "\n";
  dest->_mp_size = mpn_size (dd, dn);
  if (sgn < 0) mpz_neg (dest, dest);
}

#endif
} // namespace mmx