check/exam_mod_gcd.cpp

   1 /** @file exam_misc.cpp
   2  *
   3  *  Testing modular GCD.
   4  */
   5
   6 /*
   7  *  GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany
   8  *
   9  *  This program is free software; you can redistribute it and/or modify
  10  *  it under the terms of the GNU General Public License as published by
  11  *  the Free Software Foundation; either version 2 of the License, or
  12  *  (at your option) any later version.
  13  *
  14  *  This program is distributed in the hope that it will be useful,
  15  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17  *  GNU General Public License for more details.
  18  *
  19  *  You should have received a copy of the GNU General Public License
  20  *  along with this program; if not, write to the Free Software
  21  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  22  */
  23
  24 #include <cln/random.h>
  25 #include <iostream>
  26 #include <map>
  27 #include <string>
  28
  29 #include "polynomial/upoly.h"
  30 #include "polynomial/upoly_io.h"
  31 #include "polynomial/mod_gcd.h"
  32 #include "ginac.h"
  33 using namespace GiNaC;
  34
  35 static upoly ex_to_upoly(const ex& e, const symbol& x);
  36 static ex upoly_to_ex(const upoly& p, const symbol& x);
  37
  38 // make a univariate polynomial \in Z[x] of degree deg
  39 static upoly make_random_upoly(const std::size_t deg);
  40
  41 static void run_test_once(const std::size_t deg)
  42 {
  43         static const symbol xsym("x");
  44
  45         const upoly a = make_random_upoly(deg);
  46         const upoly b = make_random_upoly(deg);
  47
  48         upoly g;
  49         mod_gcd(g, a, b);
  50
  51         ex ea = upoly_to_ex(a, xsym);
  52         ex eb = upoly_to_ex(b, xsym);
  53         ex eg = gcd(ea, eb);
  54
  55         const upoly g_check = ex_to_upoly(eg, xsym);
  56         if (g != g_check) {
  57                 std::cerr << "a = " << a << std::endl;
  58                 std::cerr << "b = " << b << std::endl;
  59                 std::cerr << "mod_gcd(a, b) = " << g << std::endl;
  60                 std::cerr << "sr_gcd(a, b) = " << g_check << std::endl;
  61                 throw std::logic_error("bug in mod_gcd");
  62         }
  63 }
  64
  65 int main(int argc, char** argv)
  66 {
  67         std::cout << "examining modular gcd. ";
  68         std::map<std::size_t, std::size_t> n_map;
  69         // run 256 tests with polynomials of degree 10
  70         n_map[10] = 256;
  71         // run 32 tests with polynomials of degree 100
  72         n_map[100] = 32;
  73         std::map<std::size_t, std::size_t>::const_iterator i = n_map.begin();
  74         for (; i != n_map.end(); ++i) {
  75                 for (std::size_t k = 0; k < i->second; ++k)
  76                         run_test_once(i->first);
  77         }
  78         return 0;
  79 }
  80
  81 static upoly ex_to_upoly(const ex& e, const symbol& x)
  82 {
  83         upoly p(e.degree(x) + 1);
  84         for (int i = 0; i <= e.degree(x); ++i)
  85                 p[i] = cln::the<cln::cl_I>(ex_to<numeric>(e.coeff(x, i)).to_cl_N());
  86         return p;
  87 }
  88
  89 static ex upoly_to_ex(const upoly& p, const symbol& x)
  90 {
  91         exvector tv(p.size());
  92         for (std::size_t i = 0; i < p.size(); ++i)
  93                 tv[i] = pow(x, i)*numeric(p[i]);
  94         return (new add(tv))->setflag(status_flags::dynallocated);
  95 }
  96
  97 static upoly make_random_upoly(const std::size_t deg)
  98 {
  99         static const cln::cl_I biggish("98765432109876543210");
 100         upoly p(deg + 1);
 101         for (std::size_t i = 0; i <= deg; ++i)
 102                 p[i] = cln::random_I(biggish);
 103
 104         // Make sure the leading coefficient is non-zero
 105         while (zerop(p[deg]))
 106                 p[deg] = cln::random_I(biggish);
 107         return p;
 108 }