/******************************************************************************
* MODULE     : p_adic_test.cpp
* DESCRIPTION: Test p_adic numbers
* 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 <algebramix/p_adic.hpp>

using namespace mmx;

#define ZZp modular<modulus<integer> >

#define F1 modular<modulus<signed char, modulus_int_preinverse<3> >,    \
                   modular_fixed<int,2> >

#define F3 modular<modulus<signed char, modulus_int_preinverse<3> >,    \
                   modular_fixed<int,7> >

#define F29 modular<modulus<int, modulus_int_preinverse<29> >,  \
                    modular_fixed<int,268435459> >

#define F30 modular<modulus<int, modulus_int_preinverse<30> >,  \
                    modular_fixed<int,536870923> >

#define VN series_carry_p_adic<series_carry_naive>
#define VL series_carry_lift<series_carry_naive>
#define VR series_carry_relaxed<series_carry_naive> 
#define VM series_carry_modular_int_naive<series_carry_naive>
#define VO4 series_carry_monoblock<VN,4,VN,4>
#define VO32 series_carry_monoblock<VN,4,VN,4>
#define VB4 series_carry_blocks<VN,8,VN,4>
#define VB32 series_carry_blocks<VN,32,VN,32>

#define P_adic series<Modular, V>

#define TMPL template<typename Modular, typename V>

TMPL void
simple_test () {
  P_adic::set_output_order (10);

  P_adic x (101), y, z;
  mmout << "x= " << x << "\n";
  mmout << "-x= " << -x << "\n";
  mmout << "x + (-x)= " << x + (-x) << "\n";
  mmout << "x * x= " << x * x << "\n";
  mmout << "5 * x= " << Modular(5) * x << "\n";
  mmout << "(5 * x) / 5= " << (Modular(5) * x) / Modular(5) << "\n";

  P_adic p (vec (Modular (0), Modular (1)));
  P_adic a = 1 - p;
  mmout << "a     = " << a << "\n";
  P_adic b = 1 / a;
  mmout << "1 / a   = " << b << "\n";
  mmout << "a * (1 / a) = " << a * b << "\n";

  mmout << "a / b = " << a / b << "\n";
  mmout << "(a / b) * b = " << (a / b) * b << "\n";
  mmout << "(a / b) * b - a = " << (a / b) * b - a << "\n\n";

  mmout << "1 / 3 = " << P_adic (1) / P_adic (3) << "\n";
  mmout << "3 * (1 / 3) = "
        << Modular (3) * (P_adic (1) / P_adic (3)) << "\n\n";
  
  mmout << "1 / 5 = " << P_adic (1) / P_adic (5) << "\n";
  mmout << "5 * (1 / 5) = " << 5 * (P_adic (1) / P_adic (5)) << "\n\n";
}

template<typename Op,typename Modular,typename V> void
test_binary (nat n) {
  typedef series<Modular, series_carry_naive> P_adic_naive;
  mmout << "Testing " << Op::name () << " with "
        << cpp_demangle (typeid (P_adic).name ()) << "\n";
  for (nat i = 1; i < n; i++) {
    mmout << "in size " << i << "\r" << flush_now;
    P_adic::set_output_order (i);
    P_adic_naive::set_output_order (i);
    P_adic::set_cancel_order (i);
    vector<Modular> v (Modular (), i);
    for (nat j= 0; j < i; j++) v[j]= Modular (rand ());
    P_adic x (v); P_adic_naive xx (v);
    for (nat j= 0; j < i; j++) v[j]= Modular (rand ());
    v[0]= 1;
    
    P_adic y (v); P_adic_naive yy (y);
    P_adic z= Op::op (x, y);
    (void) z[i-1];
    P_adic_naive zz= Op::op (xx, yy);
    if (as<P_adic> (zz) != z) {
      mmerr << "\nError with \n"
            << "x= " << x << "\n"
            << "y= " << y << "\n"
            << "x * y=     " << z << "\n"
            << "instead of " << zz << "\n";
      exit (-1);
    }
  }
  mmout << "\n";
}

TMPL void
test_ring (nat n) {
  test_binary<sub_op,Modular,V> (n);
  test_binary<add_op,Modular,V> (n);
  test_binary<mul_op,Modular,V> (n);
  test_binary<div_op,Modular,V> (n);
}

TMPL void
test_separable_root (int r) {
  mmout << "Testing separable root with "
        << cpp_demangle (typeid (P_adic).name ()) << ", r= " << r << "\n";
  integer p= * Modular::get_modulus ();
  P_adic::set_cancel_order (50);
  vector<Modular> v (Modular (), 50);
  v[0]= 3; for (nat j= 1; j < 50; j++) v[j]= Modular (rand ());
  P_adic x (v); x = binpow (x, r);
  P_adic y = separable_root (x, r);
  P_adic z = x - binpow (y, r);
  if (z != P_adic (Modular (0))) {
    mmerr << "x= " << x << "\n";
    mmerr << "y= " << y << "\n";
    mmerr << "x - y^" << p << " = " << z << "\n";
    exit (-1);
  }
}

TMPL void
test_pth_root () {
  integer p= * Modular::get_modulus ();
  mmout << "Testing p-th root with "
        << cpp_demangle (typeid (P_adic).name ()) << "\n";
  P_adic::set_cancel_order (50);
  vector<Modular> v (Modular (), 50);
  v[0]= (p+1) >> 1; for (nat j= 1; j < 50; j++) v[j]= Modular (j);
  P_adic x (v); x = binpow (x, p);
  P_adic y = pth_root (x);
  P_adic z = x - binpow (y, p);
  if (z != P_adic (Modular (0))) {
    mmerr << "x= " << x << "\n";
    mmerr << "y= " << y << "\n";
    mmerr << "x - y^" << p << " = " << z << "\n";
    exit (-1);
  }
}

int
main () {
  integer p (536870923);
  integer q= binpow (p, 32);
  ZZp::set_modulus (q);

  try {
    //simple_test<F3,VN> ();

    test_ring<F3,VB32> (100);
    test_ring<F29,VB32> (100);
    test_ring<F30,VB32> (100);
    //test_ring<ZZp,VB32> (100);

    test_separable_root<F3,VB32> (2);
    test_separable_root<F3,VB32> (3);
    test_separable_root<F3,VB32> (5);
    test_separable_root<F3,VB32> (13);
    test_separable_root<F3,VB32> (29);

    test_pth_root<F1,VB32> ();
    test_pth_root<F3,VB32> ();
    test_pth_root<F29,VB32> ();
    test_pth_root<F30,VB32> ();
  }
  catch (exception& e) {
    mmerr << e << "\n";
    return -1;
  }
}