/Users/mourrain/Devel/mmx/algebramix/include/algebramix/crt_int.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : crt_int.hpp
* DESCRIPTION: CRT subroutines for hardware integers
* COPYRIGHT  : (C) 2008  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.
******************************************************************************/

#ifndef __MMX_CRT_INT_HPP
#define __MMX_CRT_INT_HPP
#include <numerix/integer.hpp>
#include <numerix/modular_int.hpp>
#include <algebramix/crt_naive.hpp>

namespace mmx {

/******************************************************************************
* Variant for hardware integers
******************************************************************************/

DEFINE_VARIANT(crt_int,crt_signed<crt_naive>)
DEFINE_VARIANT(crt_unsigned_int,crt_naive)

#define DECLARE_HELPER(C)                       \
  STMPL struct crt_naive_variant_helper<C> {    \
    typedef crt_int CV; };
DECLARE_HELPER (signed char)
DECLARE_HELPER (short int)
DECLARE_HELPER (int)
DECLARE_HELPER (long int)
DECLARE_HELPER (long long int)
#undef DECLARE_HELPER

#define DECLARE_HELPER(C)                       \
  STMPL struct crt_naive_variant_helper<C> {    \
    typedef crt_unsigned_int CV; };
DECLARE_HELPER (unsigned char)
DECLARE_HELPER (unsigned short int)
DECLARE_HELPER (unsigned int)
DECLARE_HELPER (unsigned long int)
DECLARE_HELPER (unsigned long long int)
#undef DECLARE_HELPER

/******************************************************************************
* Special transform taking care of casual overflows
******************************************************************************/

template<typename V>
struct implementation<crt_transform,V,crt_int> :
    public implementation<crt_project,V> {
  typedef implementation<crt_project,V> Crt;

  template<typename C, typename M, typename I> static inline void
  direct (I* c, const C& a, const M* p, nat n) {
    for (nat i= 0; i < n; i++)
      c[i]= Crt::mod (a, p[i]); }

  template<typename C, typename M, typename I, typename K> static inline void
  inverse (C& a, const I* c, const M* p, const K& P,
           const C* q, const I* m, nat n) {
    a= 0; I t;
    for (nat i= 0; i < n; i++) {
      mul_mod (t, m[i], c[i], p[i]);
      add_mod (a, t * q[i], P); } }
};

template<typename V>
struct implementation<crt_transform,V,crt_unsigned_int> :
    public implementation<crt_transform,V,crt_int> {};

/******************************************************************************
* Sequence of prime numbers
******************************************************************************/

#define TMPL template<typename M>

struct prime_sequence_int {
private:
  TMPL static inline bool
  is_prime (const vector<M>& v, nat n) {
    for (nat i= 0; true; i++) {
      nat p= (nat) *v[i];
      if (n % p == 0) return false;
      if (p*p >= n) break;
    }
    return true; }

public:
  TMPL static bool
  extend (vector<M>& v, nat i) {
    if (N(v) == 0) v= vec (M(2), M(3));
    nat p= *(v[N(v) - 1]) + 2;
    while (N(v) < i) {
      if (is_prime (v, p)) v << M(p);
      p += 2;
    }
    return true; }
};

#define DECLARE_HELPER(I)                               \
template<typename V>                                    \
 struct coprime_moduli_helper<modulus<I, V> > {         \
  typedef prime_sequence_int sequence;                  \
};
DECLARE_HELPER (signed char)
DECLARE_HELPER (unsigned char)
DECLARE_HELPER (short int)
DECLARE_HELPER (unsigned short int)
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

/******************************************************************************
* Covering
******************************************************************************/

template<typename C, typename M, typename W>
struct moduli_unsigned_integer_helper {
  // C is an unsigned int type
  template<typename V> static bool
  covering (vector<M, V>& v, nat s) {
    static coprime_moduli_sequence<M,W> seq;
    v= vector<M, V> ();
    C p (1), q; nat i= 0;
    while (bit_size (p) < s + 1) {
      if (seq [i] == 0) return false;
      v << seq [i];
      q = p * C (* seq [i]);
      if (quo (q, p) != C (* seq [i])) return false; // means overflow!
      p = q;
      i++;
    }
    return true; }
};

template<typename C, typename M, typename W>
struct moduli_signed_integer_helper {
  // C is a signed int type
  template<typename V> static bool
  covering (vector<M, V>& v, nat s) {
    static coprime_moduli_sequence<M,W> seq;
    v= vector<M, V> ();
    C p (1), q; nat i= 0;
    while (bit_size (p) < s + 2) { // s + 2 is needed for sign decoding!
      if (seq [i] == 0) return false;
      v << seq [i];
      q = p * C (* (seq [i]));
      if (quo (q, p) != C (* seq [i])) return false; // means overflow!
      p = q;
      i++;
    }
    return true; }
};

#define DECLARE_HELPER(C)                                               \
  template<typename M, typename W>                                      \
  struct moduli_helper<C,M,W> : moduli_signed_integer_helper<C,M,W> {};
DECLARE_HELPER (signed char)
DECLARE_HELPER (short int)
DECLARE_HELPER (int)
DECLARE_HELPER (long int)
DECLARE_HELPER (long long int)
#undef DECLARE_HELPER

#define DECLARE_HELPER(C)                                               \
  template<typename M, typename W>                                      \
  struct moduli_helper<C,M,W> : moduli_unsigned_integer_helper<C,M,W> {};
DECLARE_HELPER (unsigned char)
DECLARE_HELPER (unsigned short int)
DECLARE_HELPER (unsigned int)
DECLARE_HELPER (unsigned long int)
DECLARE_HELPER (unsigned long long int)
#undef DECLARE_HELPER

/******************************************************************************
* Naive primality test
******************************************************************************/

bool is_prime (nat n);
bool is_probable_prime (unsigned long int n);

/******************************************************************************
* Sequence of prime numbers of bit size s
******************************************************************************/

template <nat s>
struct fft_prime_sequence_int {
  // prime numbers of bit-size at most s,
  // starting with primes suitable to fft.
private:
  static nat
  bit_reverse (nat i, nat n) {
    if (n <= 1) return i;
    nat n2= n >> 1;
    nat n1= n - n2;
    return (bit_reverse (i & ((1 << n2) - 1), n2) << n1) +
      bit_reverse (i >> n2, n1);
  }
public:
  template<typename C, typename V> static bool
  extend (vector<modulus<C,V> >& v, nat n) {
    // construct prime moduli of s bits exactly,
    // starting with those which have a maximal number of
    // primitive 2^k-th roots of unity
    typedef modulus<C,V> M;
    const nat b= V::template maximum_size_helper<C>::value;
    ASSERT (s <= b, "bitsize overflow");
    ASSERT (b <= 8 * sizeof (nat), "bitsize overflow");
    if (s < 2 && n > 0) return false;
    const C m= 1 + MMX_SAFE_LEFT_SHIFT_INT (C, 1, (s-1));
    const C l= MMX_SAFE_LEFT_SHIFT_INT (C, 1, s);
    for (nat i= N(v); i < n; i++) {
      C j= i == 0 ? m : (bit_reverse (* v[i-1], s) + 2), k;
      for (; j <= l - 1 && j >= m; j += 2)
        if (is_prime (k= bit_reverse (j, s))) break;
      if (j > l - 1 || j < m) return false;
      v << M(k); }
    return true; }
};

/******************************************************************************
* Sequence of probable prime numbers of bit size s
******************************************************************************/

template <nat s>
struct probable_prime_sequence_int {
  // probable prime numbers of bit-size at least s

  template<typename C, typename V> static bool
  extend (vector<modulus<C,V> >& v, nat n) {
    // construct probable prime moduli of s bits exactly.
    typedef modulus<C,V> M;
    const nat b= V::template maximum_size_helper<C>::value;
    ASSERT (s <= b, "bitsize overflow");
    if (s < 2 && n > 0) return false;
    integer a (N(v) == 0 ? (integer (1) << (s-1)) : integer (* v[N(v)-1]));
    integer p= probable_next_prime (a);
    if (bit_size (p) > b) return false;
    v << M(as <C> (p));
    return true; }
};

#undef TMPL
} // namespace mmx
#endif // __MMX_CRT_INT_HPP