check/exam_mod_gcd.cpp

   1 #include "polynomial/upoly.hpp"
   2 #include "polynomial/upoly_io.hpp"
   3 #include "polynomial/mod_gcd.hpp"
   4 #include "ginac.h"
   5 #include <cln/random.h>
   6 #include <string>
   7 #include <iostream>
   8 #include <map>
   9 using namespace GiNaC;
  10
  11 static upoly ex_to_upoly(const ex& e, const symbol& x);
  12 static ex upoly_to_ex(const upoly& p, const symbol& x);
  13
  14 // make a univariate polynomial \in Z[x] of degree deg
  15 static upoly make_random_upoly(const std::size_t deg);
  16
  17 static void run_test_once(const std::size_t deg)
  18 {
  19         static const symbol xsym("x");
  20
  21         const upoly a = make_random_upoly(deg);
  22         const upoly b = make_random_upoly(deg);
  23
  24         upoly g;
  25         mod_gcd(g, a, b);
  26
  27         ex ea = upoly_to_ex(a, xsym);
  28         ex eb = upoly_to_ex(b, xsym);
  29         ex eg = gcd(ea, eb);
  30
  31         const upoly g_check = ex_to_upoly(eg, xsym);
  32         if (g != g_check) {
  33                 std::cerr << "a = " << a << std::endl;
  34                 std::cerr << "b = " << b << std::endl;
  35                 std::cerr << "mod_gcd(a, b) = " << g << std::endl;
  36                 std::cerr << "sr_gcd(a, b) = " << g_check << std::endl;
  37                 throw std::logic_error("bug in mod_gcd");
  38         }
  39 }
  40
  41 int main(int argc, char** argv)
  42 {
  43         std::cout << "examining modular gcd. ";
  44         std::map<std::size_t, std::size_t> n_map;
  45         // run 256 tests with polynomials of degree 10
  46         n_map[10] = 256;
  47         // run 32 tests with polynomials of degree 100
  48         n_map[100] = 32;
  49         std::map<std::size_t, std::size_t>::const_iterator i = n_map.begin();
  50         for (; i != n_map.end(); ++i) {
  51                 for (std::size_t k = 0; k < i->second; ++k)
  52                         run_test_once(i->first);
  53         }
  54         return 0;
  55 }
  56
  57 static upoly ex_to_upoly(const ex& e, const symbol& x)
  58 {
  59         upoly p(e.degree(x) + 1);
  60         for (int i = 0; i <= e.degree(x); ++i)
  61                 p[i] = cln::the<cln::cl_I>(ex_to<numeric>(e.coeff(x, i)).to_cl_N());
  62         return p;
  63 }
  64
  65 static ex upoly_to_ex(const upoly& p, const symbol& x)
  66 {
  67         exvector tv(p.size());
  68         for (std::size_t i = 0; i < p.size(); ++i)
  69                 tv[i] = pow(x, i)*numeric(p[i]);
  70         return (new add(tv))->setflag(status_flags::dynallocated);
  71 }
  72
  73 static upoly make_random_upoly(const std::size_t deg)
  74 {
  75         static const cln::cl_I biggish("98765432109876543210");
  76         upoly p(deg + 1);
  77         for (std::size_t i = 0; i <= deg; ++i)
  78                 p[i] = cln::random_I(biggish);
  79
  80         // Make sure the leading coefficient is non-zero
  81         while (zerop(p[deg]))
  82                 p[deg] = cln::random_I(biggish);
  83         return p;
  84 }
  85