polynomial_int_test.cpp

/******************************************************************************
* MODULE     : polynomial_int_test.cpp
* DESCRIPTION: Test polynomials
* COPYRIGHT  : (C) 2004  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.
******************************************************************************/

#include <numerix/rational.hpp>
#include <numerix/modular_integer.hpp>
#include <numerix/modular_int.hpp>
#include <algebramix/polynomial.hpp>
#include <algebramix/modular_polynomial.hpp>
#include <algebramix/polynomial_tft.hpp>
#include <algebramix/polynomial_ring_dicho.hpp>
#include <algebramix/polynomial_dicho.hpp>
#include <algebramix/polynomial_int.hpp>
#include <algebramix/polynomial_integer.hpp>
#include <algebramix/polynomial_quotient.hpp>
#include <algebramix/polynomial_polynomial.hpp>
#include <algebramix/polynomial_modular_int.hpp>
#include <algebramix/polynomial_modular_integer.hpp>

using namespace mmx;

#define TMPL template<typename C, typename V>
#define Polynomial polynomial<C,V>

#define TMPLW template<typename C, typename V, typename W>
#define Polynomial_polynomial polynomial<Polynomial,W>

/******************************************************************************
* Polynomial families
******************************************************************************/

struct sample {
  TMPL static void
  op (Polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial (1, deg);
    for (int i=0; i<deg; i++)
      P += Polynomial (C(random () - C(random ()) + 1), i); }

  TMPLW static void
  op (Polynomial_polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial_polynomial (Polynomial(1), deg);
    for (int i=0; i<deg; i++) {
      Polynomial c; op (c, deg);
      P += Polynomial_polynomial (c, i); } }
};

struct sample0 {
  TMPL static void
  op (Polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial (1, deg);
    for (int i=0; i<deg; i++)
      P += Polynomial (C(i + 1), i); }

  TMPLW static void
  op (Polynomial_polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial_polynomial (Polynomial(1), deg);
    for (int i=0; i<deg; i++) {
      Polynomial c; op (c, deg);
      P += Polynomial_polynomial (c, i); } }
};

template<typename S>
struct sample_vector {
  TMPL static vector<Polynomial>
  op (int deg) {
    vector<Polynomial> v;
    while (deg > 0) {
      nat d= 1 + abs (rand()) % deg;
      Polynomial (p); S::op (p, d);
      v << p;
      deg -= d;
    }
    return v; }
};

// Special families for subresultants

struct sample_subresultant0 {
  TMPL static void
  op (Polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial (2, deg);
    for (int i=0; i<deg; i++)
      P += Polynomial (C(i + 1), i); }

  TMPLW static void
  op (Polynomial_polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial_polynomial (Polynomial(1), deg);
    for (int i=0; i<deg; i++) {
      Polynomial c; op (c, deg);
      P += Polynomial_polynomial (c, i); } }
};

struct sample_subresultant1 {
  TMPL static void
  op (Polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial (2, deg);
    for (int i=0; i<deg; i++)
      P += Polynomial (C(i) * C(i), i); }

  TMPLW static void
  op (Polynomial_polynomial& P, int deg) {
    if (deg == -1) { P= Polynomial (); return; }
    P= Polynomial_polynomial (Polynomial(1), deg);
    for (int i=0; i<deg; i++) {
      Polynomial c; op (c, deg);
      P += Polynomial_polynomial (c, i); } }
};

/******************************************************************************
* Tests by signature
******************************************************************************/

template<typename Op, typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_unary_op (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  Polynomial0 p0, r0;
  Polynomial1 p1, r1;
  mmout << "Testing " << Op::name ()
        << " of polynomials over " << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  for (int i = -1; i <= max_deg; i++) {
    mmout << "                                   \r";
    mmout << "  up to degree " << i << "\r" << flush_now;
    S0::op (p0, i); p1 = as<Polynomial1> (p0);
    r0 = Op::op (p0); r1 = Op::op (p1);
    if (r1 != as<Polynomial1> (r0)) {
      mmout << "\nError with:\n";
      mmout << "p0 = " << p0 << "\n";
      mmout << "r0 = " << r0 << "\n";
      mmout << "p1 = " << p1 << "\n";
      mmout << "r1 = " << r1 << "\n";
      exit (-1);
    }
  }
  mmout << "\n";
}

template<typename Op, typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_binary_op (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing " << Op::name ()
        << " of polynomials over " << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, q0, r0;
  Polynomial1 p1, q1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = -1; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i); S1::op (q0, j);
      p1 = as<Polynomial1> (p0); q1 = as<Polynomial1> (q0);
      r0 = Op::op (p0, q0); r1 = Op::op (p1, q1);
      if (r1 != as<Polynomial1> (r0)) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "q0 = " << q0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "r1 = " << r1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_div (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing exact division of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, q0, r0;
  Polynomial1 p1, q1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = 0; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i); S1::op (q0, j); p0 *= q0;
      p1 = as<Polynomial1> (p0); q1 = as<Polynomial1> (q0);
      r0 = p0 / q0; r1 = p1 / q1;
      if (r1 != as<Polynomial1> (r0)) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "q0 = " << q0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "r1 = " << r1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_gcd (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing gcd of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, q0, r0;
  Polynomial1 p1, q1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = -1; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i); S1::op (q0, j);
      p1 = as<Polynomial1> (p0); q1 = as<Polynomial1> (q0);
      r0 = gcd (p0, q0); r1 = gcd (p1, q1);
      if (deg (r0) != deg (r1) || prem (r1, as<Polynomial1> (r0)) != 0) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "q0 = " << q0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "r1 = " << r1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_multi_gcd (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing simultaneous gcd of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0; vector<Polynomial0> q0, r0;
  Polynomial1 p1; vector<Polynomial1> q1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = -1; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i);
      q0 = sample_vector<S1>::template op<C,V0> (1 + (j >> 2));
      p1 = as<Polynomial1> (p0);
      q1 = as<vector<Polynomial1> > (q0);
      r0 = gcd (p0, q0); r1 = gcd (p1, q1);
      for (nat k= 0; k < N(q0); k++)
        if (deg (r0[k]) != deg (r1[k])
            || prem (r1[k], as<Polynomial1> (r0[k])) != 0) {
          mmout << "\nError with:\n";
          mmout << "p0 = " << p0 << "\n";
          mmout << "q0 = " << q0 << "\n";
          mmout << "r0 = " << r0 << "\n";
          mmout << "r1 = " << r1 << "\n";
          exit (-1);
        }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_xgcd (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing gcd of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, q0, r0, u0, v0;
  Polynomial1 p1, q1, r1, u1, v1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = -1; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i); S1::op (q0, j);
      p1 = as<Polynomial1> (p0); q1 = as<Polynomial1> (q0);
      r0 = gcd (p0, q0, u0, v0); r1 = gcd (p1, q1, u1, v1);
      if (deg (r0) != deg (r1) || prem (r1, as<Polynomial1> (r0)) != 0 ||
          deg (u0) != deg (u1) || prem (u1, as<Polynomial1> (u0)) != 0 ||
          deg (v0) != deg (v1) || prem (v1, as<Polynomial1> (v0)) != 0) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "q0 = " << q0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "u0 = " << u0 << "\n";
        mmout << "v0 = " << v0 << "\n";
        mmout << "r1 = " << r1 << "\n";
        mmout << "u1 = " << u1 << "\n";
        mmout << "v1 = " << v1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_lcm (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing lcm of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, q0, r0;
  Polynomial1 p1, q1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = -1; j <= max_deg; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i << " x " << j << "\r" << flush_now;
      S0::op (p0, i); S1::op (q0, j);
      p1 = as<Polynomial1> (p0); q1 = as<Polynomial1> (q0);
      r0 = lcm (p0, q0); r1 = lcm (p1, q1);
      if (deg (r1) != deg (r0) || prem (r1, as<Polynomial1> (r0)) != 0) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "q0 = " << q0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "r1 = " << r1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0> void
test_pade (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing pade approximants from polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, num0, den0;
  Polynomial1 p1, num1, den1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = 1; j <= i; j++) {
      mmout << "                                   \r";
      mmout << "  up to degree " << i
           << " with numerator of degree less than " << j << "\r";
      S0::op (p0, i); p1 = as<Polynomial1> (p0);
      pade (p0, i+1, j, num0, den0);
      pade (p1, i+1, j, num1, den1);
      if (deg (num1) != deg (num0) ||
          prem (num1, as<Polynomial1> (num0)) != 0 ||
          deg (den1) != deg (den0) ||
          prem (den1, as<Polynomial1> (den0)) != 0) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "num0 = " << num0 << "\n";
        mmout << "den0 = " << den0 << "\n";
        mmout << "num1 = " << num1 << "\n";
        mmout << "den1 = " << den1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1, typename S0> void
test_evaluate (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing evaluation of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 p0, r0; vector<C> x;
  Polynomial1 p1, r1;
  for (int i = -1; i <= max_deg; i++)
    for (int j = 0; j <= 2 * i; j++) {
      mmout << "  up to degree " << i << " at " << j << " points\r";
      mmout << flush_now;
      S0::op (p0, i); p1 = as<Polynomial1> (p0);
      x= vector<C> (C(0), (nat) j);
      for (int k= 0; k < j; k++) x[k]= k * k;
      r0 = evaluate (p0, x); r1 = evaluate (p1, x);
      if (r1 != as<Polynomial1> (r0)) {
        mmout << "\nError with:\n";
        mmout << "p0 = " << p0 << "\n";
        mmout << "r0 = " << r0 << "\n";
        mmout << "p1 = " << p1 << "\n";
        mmout << "r1 = " << r1 << "\n";
        exit (-1);
      }
    }
  mmout << "\n";
}

template<typename C, typename V0, typename V1> void
test_interpolate (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing interpolation of polynomials over "
        << cpp_demangle (typeid(C).name ())
        << " with " << cpp_demangle (typeid(V0).name ())
        << " versus " << cpp_demangle (typeid(V1).name ()) << "\n";
  Polynomial0 r0; Polynomial1 r1; vector<C> x, y;
  for (int i = -1; i <= max_deg; i++) {
    mmout << "  up to degree " << i << "\r" << flush_now;
    x= vector<C> (C(0), (nat) (i+1));
    y= vector<C> (C(0), (nat) (i+1));
    for (int k= 0; k < i+1; k++) x[k]= k;
    if (big_mul (x) == 0) continue;
    for (int k= 0; k < i+1; k++) y[k]= k * k;
    r0 = interpolate_bis<C,V0> (y, x);
    r1 = interpolate_bis<C,V1> (y, x);
    if (r1 != as<Polynomial1> (r0)) {
      mmout << "\nError with:\n";
      mmout << "x = " << x << "\n";
      mmout << "y = " << y << "\n";
      mmout << "r0 = " << r0 << "\n";
      mmout << "r1 = " << r1 << "\n";
      exit (-1);
    }
  }
  mmout << "\n";
}

template<typename C, typename V0, typename V1,
         typename S0, typename S1> void
test_subresultant (int max_deg) {
  typedef polynomial<C,V0> Polynomial0;
  typedef polynomial<C,V1> Polynomial1;
  mmout << "Testing subresultants over "
        << cpp_demangle (typeid (C).name ())
        << " with " << cpp_demangle (typeid (V0).name ())
        << " and " << cpp_demangle (typeid (V1).name ()) << "\n";
  Polynomial0 f0, g0, d0;
  Polynomial1 f1, g1, d1;
  for (int i = -1; i < max_deg; i++) {
    for (int j = -1; j < max_deg; j++) {
      mmout << "  in degrees " << i << " x " << j << "  \r" << flush_now;
      nat k= (nat) max (min (i,j), 0);
      Polynomial0 one0 (promote (1, CF(f0)));
      Polynomial1 one1 (promote (1, CF(f1)));
      vector<Polynomial0> s0 (one0, k), co0_f (one0, k), co0_g (one0, k);
      vector<Polynomial1> s1 (one1, k), co1_f (one1, k), co1_g (one1, k);
      S0::op (f0, i); S1::op (g0, j);
      f1 = as<Polynomial1> (f0); g1 = as<Polynomial1> (g0);
      s0= subresultants (f0, g0, co0_f, co0_g);
      s1= subresultants (f1, g1, co1_f, co1_g);
      if (s0 != as<vector<Polynomial0> > (s1) ||
          co0_f != as<vector<Polynomial0> > (co1_f) ||
          co0_g != as<vector<Polynomial0> > (co1_g) ||
          f0 * co0_f + g0 * co0_g != s0) { 
        mmout << "\nError: \n";
        mmout << "f= " << f0 << "\n";
        mmout << "g= " << g0 << "\n";
        mmout << "s0= " << s0 << "\n";
        mmout << "co0_f= " << co0_f << "\n";
        mmout << "co0_g= " << co0_g << "\n";
        mmout << "s1= " << s1 << "\n";
        mmout << "co1_f= " << co1_f << "\n";
        mmout << "co1_g= " << co1_g << "\n";
        exit (-1);
      }
    }
  }
  mmout << "\n";
}

/******************************************************************************
* Tests by category
******************************************************************************/

template<typename C, typename V0, typename V1> void
test_over_ring (int max_deg) {
  typedef sample S;
  test_unary_op<square_op  , C, V0, V1, S, S> (max_deg);
  test_unary_op<graeffe_op , C, V0, V1, S, S> (max_deg);
  test_binary_op<add_op    , C, V0, V1, S, S> (max_deg);
  test_binary_op<sub_op    , C, V0, V1, S, S> (max_deg);
  test_binary_op<mul_op    , C, V0, V1, S, S> (max_deg);
  test_binary_op<prem_op   , C, V0, V1, S, S> (max_deg);
  test_binary_op<pquo_op   , C, V0, V1, S, S> (max_deg);
  test_binary_op<compose_op, C, V0, V1, S, S> (max_deg);
  test_evaluate<C, V0, V1, S> (max_deg);
  test_subresultant<C, V0, V1,
    sample_subresultant0, sample_subresultant1> (max_deg);
}

template<typename C, typename V0, typename V1> void
test_over_gcd_ring (int max_deg) {
  typedef sample S;
  test_over_ring<C, V0, V1> (max_deg);
  test_div<C, V0, V1, S, S> (max_deg);
  test_gcd<C, V0, V1, S, S> (max_deg);
  test_multi_gcd<C, V0, V1, S, S> (max_deg);
  test_lcm<C, V0, V1, S, S> (max_deg);
  test_pade<C, V0, V1, S> (max_deg);
}

template<typename C, typename V0, typename V1> void
test_over_field (int max_deg) {
  typedef sample S;
  test_over_gcd_ring<C, V0, V1> (max_deg);
  test_xgcd<C, V0, V1, S, S> (max_deg);
  test_binary_op<rem_op, C, V0, V1, S, S> (max_deg);
  test_binary_op<quo_op, C, V0, V1, S, S> (max_deg);
  test_interpolate<C, V0, V1> (max_deg);
}

int
main () {
try {
  int m= 10;

  test_over_ring<unsigned char,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<short int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);
  /*
  test_over_ring<int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);
  */
  test_over_ring<long int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);
  test_over_ring<long long int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<unsigned char,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<unsigned short int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<unsigned int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<unsigned long int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);

  test_over_ring<unsigned long long int,
    polynomial_ring_naive<polynomial_naive>,
    polynomial_int> (m);
}
catch (exception& e) {
  mmout << e << "\n";
  return -1;
}
return 0;
}