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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_integer.cpp
* DESCRIPTION: Subroutines for efficient Chinese remaindering with GMP
* 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.
******************************************************************************/

#include <basix/port.hpp>
#include <algebramix/crt_integer.hpp>
namespace mmx {
#if defined (__GNU_MP__)
typedef mp_limb_t C;

/******************************************************************************
* Extra routine for debugging
******************************************************************************/

/*
static void
mpn_print (const port& out, const C* src, nat n) {
  mpz_t i;
  i->_mp_alloc= n;
  i->_mp_size = n;
  i->_mp_d    = const_cast<C*> (src);
  char* s= mpz_get_str (NULL, 10, i);
  out << s;
  mmx_free (s, strlen (s) + 1);  
}
*/

/******************************************************************************
* Building the auxiliary information for encoding and decoding
******************************************************************************/

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

static inline void
mpn_clear (C* dest, nat n) {
  for (; n != 0; n--) {
    *dest= 0;
    dest++;
  }
}

static inline void
mpn_copy (C* dest, const C* src, nat n) {
  for (; n != 0; n--) {
    *dest= *src;
    dest++; src++;
  }
}

static inline void
mpn_copy (C* dest, const C* src, nat n1, nat n2) {
  mpn_copy (dest, src, n2);
  mpn_clear (dest + n2, n1 - n2);
}

static void
mpn_cofactor (C* d, const C* s1, const C* s2, nat n1, nat n2) {
  // Compute the first cofactor [d, n2] of [s1, n1] and [s2, n2].
  // We have 0 <= d < s2 and d s1 + e s2 = 1 for the second cofactor e.
  mpz_t src1;
  mpz_t src2;
  mpz_t g;
  mpz_t c1;
  mpz_t c2;
  src1->_mp_alloc= n1;
  src1->_mp_size = mpn_size (s1, n1);
  src1->_mp_d    = const_cast<C*> (s1);
  src2->_mp_alloc= n2;
  src2->_mp_size = mpn_size (s2, n2);
  src2->_mp_d    = const_cast<C*> (s2);
  mpz_init (g);
  mpz_init (c1);
  mpz_init (c2);
  mpz_gcdext (g, c1, c2, src1, src2);
  if (mpz_sgn (c1) < 0) mpz_add (c1, c1, src2);
  mpn_copy (d, c1->_mp_d, n2, c1->_mp_size);
  mpz_clear (g);
  mpz_clear (c1);
  mpz_clear (c2);
}

static void
build_converters (C* enc, const C* src, nat n, nat inc) {
  // At the end, [enc, n] contains the product of all moduli.
  // If n = n1 + n2 >= 2, then 'enc + inc' and 'enc + inc + n1'
  // recursively contains the converter information for
  // the first n1 and the second n2 moduli.
  // Furthermore, [enc + inc/2, n2] contains the first cofactor
  // for the corresponding products
  //mmout << "Build " << n << ", " << inc << "\n";
  if (n == 1) *enc= (C) *src;
  else {
    nat n2= n >> 1;
    nat n1= n - n2;
    C* dec= enc + (inc >> 1);
    build_converters (enc + inc, src, n1, inc);
    build_converters (enc + inc + n1, src + n1, n2, inc);
    mpn_mul (enc, enc + inc, n1, enc + inc + n1, n2);
    mpn_cofactor (dec, enc + inc, enc + inc + n1, n1, n2);
  }
}

vector<C>
mpz_setup_crt (const vector<C>& mods) {
  nat n= N(mods);
  if (n == 0) return vector<C> ();
  nat inc= n << 1;
  nat steps= 1;
  for (nat i= n-1; i != 0; steps++, i >>= 1);
  vector<C> cv= fill<C> (0, steps * inc);
  C* enc= seg (cv);
  build_converters (enc, seg (mods), n, inc);
  return cv;
}

integer
mpz_moduli_product (const vector<C>& mods, const vector<mp_limb_t>& cv) {
  nat n= N(mods);
  if (n == 0) return 1;
  integer i= raw_integer (n);
  mpn_copy ((*i)->_mp_d, seg (cv), n);
  (*i)->_mp_size= mpn_size (seg (cv), n);
  return i;
}
  
/******************************************************************************
* Encoding and decoding integers using Chinese remaindering
******************************************************************************/

inline const C* seg (const integer& i) { return (*i)->_mp_d; }
inline C* seg (integer& i) { return (*i)->_mp_d; }

static void
mpn_mod (C* dest, C* tmp,
         const C* s1, const C* s2, nat n1, nat n2) {
  // Compute dest := s1 mod s2. The argument tmp contains auxiliary storage
  nat eff_n1= n1, eff_n2= n2;
  while (eff_n1 > 0 && s1[eff_n1-1] == 0) eff_n1--;
  while (eff_n2 > 0 && s2[eff_n2-1] == 0) eff_n2--;
  if (eff_n1 < eff_n2) mpn_copy (dest, s1, n2, eff_n1);
  else {
    mpn_tdiv_qr (tmp, dest, 0, s1, eff_n1, s2, eff_n2);
    if (n2 > eff_n2) mpn_clear (dest + eff_n2, n2 - eff_n2);
  }
}

static void
encode (C* dest, C* tmp, const C* enc, nat n, nat inc) {
  if (n == 1) return;
  nat n2= n >> 1;
  nat n1= n - n2;
  mpn_copy (tmp, dest, n);
  mpn_mod (dest, tmp + n, tmp, enc, n, n1);
  mpn_mod (dest + n1, tmp + n, tmp, enc + n1, n, n2);
  encode (dest, tmp, enc + inc, n1, inc);
  encode (dest + n1, tmp, enc + n1 + inc, n2, inc);
}

void
mpz_encode_crt (C* dest, const integer& src,
                const vector<C>& mods, const vector<C>& cv)
{
  nat n= N(mods);
  if (n == 0) return;
  const C* enc= seg (cv);
  C* tmp = mmx_new<C> (n << 1);
  mpn_copy (dest, seg (src), n, limb_size (src));
  encode (dest, tmp, enc + (n << 1), n, n << 1);
  if (sign (src) < 0)
    for (nat i=0; i<n; i++)
      if (dest[i] != 0)
        dest[i]= mods[i] - dest[i];
  mmx_delete<C> (tmp, n << 1);

  /*
  // optional check for correctness
  integer check= chinese_decode (dest, mods, cv);
  if (check != src)
    mmout << "Check failed " << src << " -> " << check << " " << mods << "\n";
  */
}

static bool
mpn_is_zero (const C* src, nat n) {
  for (nat i= 0; i<n; i++)
    if (src[i] != 0) return false;
  return true;
}

static void
mpn_reconstruct (C* dest, C* temp,
                 const C* x1, const C* x2,
                 const C* m1, const C* m2, const C* c1,
                 nat n1, nat n2)
{
  mpn_copy (temp, x2, n1, n2);
  C borrow= mpn_sub_n (dest, temp, x1, n1);
  if (borrow) mpn_sub_n (dest, x1, temp, n1);
  mpn_mul (temp, dest, n1, c1, n2);
  mpn_mod (dest, temp + n1 + n2, temp, m2, n1 + n2, n2);
  if (borrow && !mpn_is_zero (dest, n2)) {
    mpn_copy (temp, dest, n2);
    mpn_sub_n (dest, m2, temp, n2);
  }
  mpn_mul (temp, m1, n1, dest, n2);
  mpn_add (dest, temp, n1 + n2, x1, n1);
}

static void
decode (C* dest, C* tmp, const C* src,
        const C* enc, nat n, nat inc)
{
  if (n == 1) *dest= (C) *src;
  else {
    //typedef implementation<vector_linear,vector_naive> NVec;
    //mmout << "Decode "; NVec::print (mmout, src, n); mmout << "\n";
    nat n2= n >> 1;
    nat n1= n - n2;
    const C* dec= enc + (inc >> 1);
    decode (dest, tmp, src, enc + inc, n1, inc);
    decode (dest + n1, tmp, src + n1, enc + n1 + inc, n2, inc);
    mpn_copy (tmp, dest, n);
    mpn_reconstruct (dest, tmp + n, tmp, tmp + n1,
                     enc + inc, enc + inc + n1, dec, n1, n2);
    //mmout << "Decoded "; mpn_print (mmout, dest, n); mmout << "\n";
  }
}

integer
mpz_decode_crt (const C* src,
                const vector<C>& mods, const vector<C>& cv) {
  nat n= N(mods);
  if (n == 0) return 0;
  integer i= raw_integer (n);
  const C* enc= seg (cv);
  C* tmp = mmx_new<C> (n << 2);
  decode (seg (i), tmp, src, enc, n, n << 1);
  mpn_sub_n (tmp, enc, seg (i), n);
  if (mpn_cmp (seg (i), tmp, n) <= 0)
    (*i)->_mp_size= mpn_size (seg (i), n);
  else {
    mpn_copy (seg (i), tmp, n);
    (*i)->_mp_size= -mpn_size (seg (i), n);
  }
  mmx_delete<C> (tmp, n << 2);

  /*
  // optional check for correctness
  C* check= mmx_new<C> (n);
  chinese_encode (check, i, mods, cv);
  for (nat k=0; k<n; k++)
    if (check[k] != src[k])
      mmout << "Check failed " << src[k] << " -> " << check[k] << "\n";
  mmx_delete<C> (check, n);
  */

  return i;
}

#endif // __GNU_MP__
} // namespace mmx