welcome to the fifth (also last and easiest) part of the multivariate polynomials in c++ using templates series. this part will explain how i implemented long division and the euclidean algorithm.
long division.
long division in is defined as follows: given polynomials
, find polynomials
such that
and
. provided that the leading term of
is a unit in
(or at least divides most things),
and
exist and are unique (in fact, i'm implicitly assuming that
is commutative and unitary). this division makes
an euclidean ring in case
is a field (i.e. all non-zero elements are invertible). our
T is in general not a field, but we can still try to do long division, and let the user worry if this actually makes sense. for example, if T == int, one can divide two elements (assuming the second is not zero), but this is not the division a mathematician would expect. hence, if is not monic in case of
T == int, our implmentation will output <i>something</i> which does not necessarily make sense.
anyway. my implementation is rather straightforward, which is not very surprising:
1template<class T> 2void eucliddiv(const poly<1, T> & f, const poly<1, T> & g, poly<1, T> & q, poly<1, T> & r) 3{ 4 if (f.degree() < g.degree()) 5 { 6 // Ignore the trivial case 7 q = poly<1, T>(); 8 r = f; 9 return; 10 } 11 q = poly<1, T>(true, f.degree() - g.degree() + 1); 12 r = f; 13 for (int i = q.degree(); i >= 0; --i) 14 { 15 q[i] = r[i + g.degree()] / g.leading(); 16 for (int j = g.degree(); j >= 0; --j) 17 r[i + j] -= q[i] * g[j]; 18 } 19 r.normalize(); 20 // Note that the degree of q is already correct. 21}
here, we assume that f.leading() / g.leading() is non-zero; otherwise, q.normalize() would have to be called, but in that case the whole result will probably be broken anyway. we create the polynomial q using the special constructor mentioned here.
euclidean algorithm.
the implementation of the euclidean algorithm is not depending on the implementation anymore (long division made only use of the special constructor). first, we take care of special cases (one of the two arguments is zero), and then proceed with the main loop. we normalize the remainder in every iteration of the loop, and also the result, to (hopefully) increase numerical stability a bit. hence, for this to work, we really need that T is a field, or otherwise it will only work in very special cases.
1template<class T> 2poly<1, T> GCD(const poly<1, T> & f, const poly<1, T> & g) 3{ 4 if (f.isZero()) 5 { 6 if (g.isZero()) 7 return f; // both are zero 8 else 9 return g / g.leading(); // make g monic 10 } 11 poly<1, T> d1(f / f.leading()), d2(g / g.leading()), q, r; 12 while (!d2.isZero()) 13 { 14 eucliddiv(d1, d2, q, r); 15 d1.swap(d2); 16 d2 = r; 17 d2 /= r.leading(); // make r monic 18 } 19 return d1; 20}
i don't think there's much to write about this, which isn't already written somewhere else. this just shows how to implement a generic polynomial based algorithm.
comments.