Implemented modular GCD algorithm for univariate polynomials.
[ginac.git] / ginac / polynomial / normalize.tcc
1 #ifndef GINAC_UPOLY_NORMALIZE_TCC
2 #define GINAC_UPOLY_NORMALIZE_TCC
3 #include "upoly.hpp"
4 #include "ring_traits.hpp"
5 #include "debug.hpp"
6
7 namespace GiNaC
8 {
9
10 /// Make the univariate polynomial @a a \in F[x] unit normal.
11 /// F should be a field.
12 /// Returns true if the polynomial @x is already unit normal, and false
13 /// otherwise.
14 static bool normalize_in_field(umodpoly& a, cln::cl_MI* content_ = 0)
15 {
16         if (a.size() == 0)
17                 return true;
18         if (lcoeff(a) == the_one(a[0])) {
19                 if (content_)
20                         *content_ = the_one(a[0]);
21                 return true;
22         }
23
24         const cln::cl_MI lc_1 = recip(lcoeff(a));
25         for (std::size_t k = a.size(); k-- != 0; )
26                 a[k] = a[k]*lc_1;
27         if (content_)
28                 *content_ = lc_1;
29         return false;
30 }
31
32 /// Make the univariate polynomial @a x unit normal. This version is used
33 /// for rings which are not fields. 
34 /// Returns true if the polynomial @x is already unit normal, and false
35 /// otherwise.
36 template<typename T> bool
37 normalize_in_ring(T& x, typename T::value_type* content_ = 0, int* unit_ = 0)
38 {
39         typedef typename T::value_type ring_t;
40         static const ring_t one(1);
41         if (x.empty())
42                 return true;
43
44         bool something_changed = false;
45         if (minusp(lcoeff(x))) {
46                 something_changed = true;
47                 if (unit_)
48                         *unit_ = -1;
49                 for (std::size_t i = x.size(); i-- != 0; )
50                         x[i] = - x[i];
51         }
52
53         if (degree(x) == 0) {
54                 if (content_)
55                         *content_ = x[0];
56                 if (x[0] == one)
57                         return something_changed;
58                 x[0] = one;
59                 return false; // initial polynomial was unit normal
60         }
61                 
62         // Compute the gcd of coefficients
63         ring_t content = lcoeff(x);
64         // We want this function to be fast when applied to unit normal
65         // polynomials. Hence we start from the leading coefficient.
66         for (std::size_t i = x.size() - 1; i-- != 0; ) {
67                 if (content == one) {
68                         if (content_)
69                                 *content_ = one;
70                         return something_changed;
71                 }
72                 content = gcd(x[i], content);
73         }
74
75         if (content == one) {
76                 if (content_)
77                         *content_ = one;
78                 return something_changed;
79         }
80         
81         lcoeff(x) = one;
82         for (std::size_t i = x.size() - 1; i-- != 0; )
83                 x[i] = exquo(x[i], content);
84
85         if (content_)
86                 *content_ = content;
87         return false; // initial polynomial was not unit normal
88 }
89
90 } // namespace GiNaC
91         
92 #endif // GINAC_UPOLY_NORMALIZE_TCC
93