Implemented modular GCD algorithm for univariate polynomials.
[ginac.git] / 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