Happy New Year!
[ginac.git] / check / exam_mod_gcd.cpp
1 /** @file exam_misc.cpp
2  *
3  *  Testing modular GCD.
4  */
5
6 /*
7  *  GiNaC Copyright (C) 1999-2019 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 dynallocate<add>(tv);
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 }