crt_test.cpp

/******************************************************************************
* MODULE     : crt_test.cpp
* DESCRIPTION: Test CRT
* COPYRIGHT  : (C) 2009  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/rational.hpp>
#include <algebramix/polynomial.hpp>
#include <algebramix/crt_naive.hpp>
#include <algebramix/crt_int.hpp>
#include <algebramix/crt_integer.hpp>
#include <algebramix/crt_polynomial.hpp>
#include <algebramix/crt_blocks.hpp>

using namespace mmx;

template<typename C, typename V> void
sample (polynomial<C,V>& a, int n) {
  a= C(0);
  for (int i=0; i<n; i++) a += polynomial<C,V> (C(random ()), i);
}

void
sample (integer& a, nat n) {
  if (n == 0) { a= 0; return; }
  integer b (1);
  while (bit_size (b) < n) { a= b; b = b * random (); }
  if (random () & 1) a= -a;
}

void
sample (int& a, nat n) {
  unsigned int mask= ~ ((~ 0u) << n);
  a= ((unsigned int) random ()) & mask;
}

template<typename Crter> void
test_crt (nat m= 500) {
  typedef typename Crter::base C;
  typedef typename Crter::modulus_base I;  
  typedef typename Crter::modulus_base_variant V;  
  typedef modulus<I,V> M;

  mmout << "Testing " << typeid (Crter).name () << "\n";
  for (nat s=0; s < m;  s++) {
    mmout << "\rup to size " << s << "   " << flush_now;
    C a, b; sample (a, s);
    vector<M> p; bool t= moduli_helper<C, M>::covering (p, s);
    nat n= N(p);
    if (!t) { mmout << " --> overflow in comoduli covering!"; break; }
    Crter crter (p,false);
    vector<I> c (I (), n);
    crter.direct_transform (seg (c), a);
    crter.inverse_transform (b, seg (c));
    if (a != b) {
      mmout << "\nError detected with:\n";
      mmout << "p= " << p << "\n";
      mmout << "a= " << a << "\n";
      mmout << "b= " << b << "\n";
      exit (-1);
    }
    mmout << "\r";
  }
  mmout << "\n";
}

int
main () {
  try {
    mmout << "Default covering with primes: ";
    vector<modulus<int> > v;
    moduli_helper<int>::covering (v, 10);
    mmout << v << "\n";

    mmout << "Sequence of primes with 23 bits: ";
    coprime_moduli_sequence<modulus<int>, fft_prime_sequence_int<23> > seq;
    mmout << range (seq,0,100) << "\n";
    mmout << range (seq,99,100) << "\n";
    
    test_crt<crt_naive_transformer<int> > ();
    test_crt<crt_naive_transformer<integer> > ();
    test_crt<crt_dicho_transformer<integer> > ();
#if defined (__GNU_MP__)
    test_crt<crt_gmp_transformer<int> > ();
#endif

    test_crt<crt_blocks_transformer<
      crt_naive_transformer<integer>,
      crt_dicho_transformer<integer>,4 > > ();

  }
  catch (exception& e) {
    mmout << e << "\n";
    return -1;
  }
  return 0;
}