skip to main content.

posts about euclidean algorithm.

feeds:

atom feed. (euclidean algorithm, en) rss feed. (euclidean algorithm, en)

implementing multivariate polynomials in c++ using templates [part 5/5].

november 25, 2011 at 20:02 by felix. no comments.

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 R[X] is defined as follows: given polynomials f, g \in R[X], find polynomials q, r \in R[X] such that f = q \cdot g + r and \deg r < \deg g. provided that the leading term of g is a unit in R (or at least divides most things), q and r exist and are unique (in fact, i'm implicitly assuming that R is commutative and unitary). this division makes R[X] an euclidean ring in case R 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 g 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.

posted in: computer
tags:
places: