welcome to the fourth part of the multivariate polynomials in c++ using templates series. in this part, i want to explain how to implement efficient polynomial evaluation using a collection of interacting templates. the output will be quite optimized code (assuming the compiler has a decent enough optimizer).
evaluating multivariate polynomials is more complicated. i started with a naive implementation of operator() in poly<n, T>:
1 template<class S> 2 poly<n - 1, S> operator() (const S & x) const 3 { 4 poly<n - 1, S> res; 5 S xx = (S)1; 6 for (unsigned i = 0; i < d_value.size(); ++i) 7 { 8 res += d_value[i] * xx; 9 xx *= x; 10 } 11 return res; 12 }
if f is of type poly<2, int>, then f(4) will be of type poly<1, int>, whence f(4)(2) will call poly<1, int>::operator()<int>(const int &) to yield a polynomial of type poly<0, int>, which automatically casts to an int. unfortunately, in this process, first f(4) will be created, which requires arithmetic of polynomials of type poly<1, int>, and then the resulting polynomial will be evaluated again. in particular, if n is larger, this is far from being optimal. hence, this solution is fine if evaluations are seldomly done, but if they are done more often, it is too slow.
to stay with our example of f of type poly<2, int>. one could evaluate , written as
f(4)(2), directly as follows:
1int result = 0, xx1 = 1; 2for (int i = 0; i <= f.degree(); ++i) 3{ 4 int result2 = 0, xx2 = 1; 5 for (int j = 0; j <= f[i].degree(); ++j) 6 { 7 result2 += f[i][j] * xx2; 8 xx2 *= 2; 9 } 10 result += result2 * xx1; 11 xx1 *= 4; 12}
clearly, this is tiresome (and prone to typos) to write every time. moreover, this is also not optimal, as our operator[] will do more than just returning an element. in particular, if f is not const in this context, the compiler will insert code at every operation f[i] and f[i][j] to check whether the index is out of range (to resize d_value in that case). writing more directly
1int result = 0, xx1 = 1; 2for (int i = 0; i < f.d_value.size(); ++i) 3{ 4 int result2 = 0, xx2 = 1; 5 for (int j = 0; j < f.d_value[i].d_value.size(); ++j) 6 { 7 result2 += f.d_value[i].d_value[j] * xx2; 8 xx2 *= 2; 9 } 10 result += result2 * xx1; 11 xx1 *= 4; 12}
would result in faster code, if it would compile – d_value is private. for this reason, i came up with two templates:
1template<int n, class T, class S = T> 2class poly_evaluator; 3template<int n, class T, class HL, class S> 4class poly_evaluator_impl;
the template poly_evaluator<n, T, S> is instanciated in poly<n, T>::operator()<S>(const S &), which is defined as follows:
1 template<class S> 2 inline poly_evaluator<n, T, S> operator() (const S & x) const 3 { 4 return poly_evaluator<n, T, S>(*this, x); 5 }
the idea is that if operator() of poly_evaluator<n, T, S>(*this, x) is called, it will spawn an object of type poly_evaluator_impl<n-1, T, poly_evaluator<n, T, S>, S>. if operator() of this new object is called, it will create an object of type poly_evaluator_impl<n-2, T, poly_evaluator_impl<n-1, T, poly_evaluator<n, T, S>, S>, S>, and so on. the purpose of carrying the type of the caller around is that the outer loop of the evaluation function is the loop for the inner-most object (of type poly_evaluator<n, T, S>). if operator() would be right-associative instead of left-associative (as it is), this trouble would not be necessary.
the idea is that both templates poly_evaluator<n, T, S> and poly_evaluator<n, T, HL, S> provide an internal function template<class SS, class Fun> evaluate(SS & res, Fun evalfun), which in case of poly_evaluator<n, T, HL, S> calls the corresponding evaluate function of its owner (with evalfun replaced), and in the case of poly_evaluator<n, T, S> implements the outermost loop, which loops over the coefficients of and uses
evalfun to evaluate the coefficients. hence, the evalful handed upwards must be extended in every step to include the next variable.
assume that we have a polynomial f of type poly<3, int>, and we want to evaluate f(10)(20)(30). this is an object of type poly_evaluator_impl<1, int, poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>, int>, whose cast operator operator int() will start evaluation by defining a local class poly_evaluator_impl<1, int, poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>, int>::eval_fun, defined as follows:
1 class eval_fun 2 { 3 const poly_evaluator_impl<1, int, poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>, int> & d_owner; 4 5 public: 6 inline eval_fun(const poly_evaluator_impl<1, int, poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>, int> & owner) 7 : d_owner(owner) 8 { 9 } 10 11 inline int operator() (const poly<1, int> & p) const 12 { 13 int res = 0; 14 int xx = 1; 15 for (int i = 0; i < (int)p.d_value.size(); ++i) 16 { 17 res += p.d_value[i] * xx; 18 xx = xx * d_owner.d_evalpoint; 19 } 20 return res; 21 } 22 };
an object of it is created by operator int() as follows:
1 inline operator int() const 2 { 3 int res = 0; 4 d_owner.evaluate(res, eval_fun(*this)); 5 return res; 6 }
the object eval_fun(*this) can now be used to evaluate polynomials of type poly<1, int> at the specified evaluation point d_owner.d_evalpoint, which equals 30 in this case.
the operator int() calls poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>::evaluate<int, Fun>(int &, Fun), which in turn is implemented as follows:
1 template<class Fun> 2 inline void evaluate(int & res, const Fun & evalfun) const 3 { 4 d_owner.evaluate(res, eval_fun<int, Fun>(*this, evalfun)); 5 }
the newly created object of type eval_fun<int, Fun> will now evaluate polynomials of type poly<2, int> by using the provided functor evalfun, which is an object of the class poly_evaluator_impl<1, int, poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int>, int>::eval_fun defined above. the class eval_fun in poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int> is defined as follows:
1 template<class Fun> 2 class eval_fun 3 { 4 const poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int> & d_owner; 5 const Fun & d_evalfun; 6 7 public: 8 inline eval_fun(const poly_evaluator_impl<2, int, poly_evaluator<3, int, int>, int> & owner, const Fun & evalfun) 9 : d_owner(owner), d_evalfun(evalfun) 10 { 11 } 12 13 inline int operator() (const poly<2, int> & p) const 14 { 15 int res = 0; 16 int xx = 1; 17 for (int i = 0; i < (int)p.d_value.size(); ++i) 18 { 19 res += d_evalfun(p.d_value[i]) * xx; 20 xx = xx * d_owner.d_evalpoint; 21 } 22 return res; 23 } 24 };
note that d_owner.d_evalpoint is 20 in this case. finally, we are at the top level, namely in the function poly_evaluator<3, int, int>::evaluate<int, Fun>(), defined as follows:
1 template<class Fun> 2 inline void evaluate(int & res, const Fun & evalfun) const 3 { 4 int xx = 1; 5 for (int i = 0; i < (int)d_poly.d_value.size(); ++i) 6 { 7 res += evalfun(d_poly.d_value[i]) * xx; 8 xx = xx * d_evalpoint; 9 } 10 }
(here, d_evalpoint is 10.) this is the outer loop. the middle loop is inserted when calling evalfun(d_poly.d_value[i]), and into that the innermost loop is inserted when it calles d_evalfun(p.d_value[i]). this all essentially boils down to the following combination:
1 { 2 int res = 0; 3 int xx = 1; 4 for (int i = 0; i < (int)d_poly.d_value.size(); ++i) 5 { 6 int tmp; 7 const poly<2, int> & p = d_poly.d_value[i]; 8 { 9 int res2 = 0; 10 int xx2 = 1; 11 for (int i2 = 0; i2 < (int)p.d_value.size(); ++i2) 12 { 13 int tmp2; 14 const poly<2, int> & p2 = p.d_value[i2]; 15 { 16 int res3 = 0; 17 int xx3 = 1; 18 for (int i3 = 0; i3 < (int)p2.d_value.size(); ++i3) 19 { 20 res3 += p2.d_value[i3] * xx3; 21 xx3 = xx3 * d_owner.d_evalpoint; 22 } 23 tmp2 = res3; 24 } 25 res2 += tmp2 * xx2; 26 xx2 = xx2 * d_owner.d_evalpoint; 27 } 28 tmp = res2; 29 } 30 res += tmp * xx; 31 xx = xx * d_evalpoint; 32 } 33 return res; 34 }
of course, this is just an intermediate step (to highlight the similiarities to all the code snippets from above); the compiler will optimize this to something more compact.
this “example” already presented most of the concept. the implementation itself is more complicated, since different types can be used for evaluation (at different stages, even), and since the implementation has support for allocators.
an important point is that both templates poly_evaluator<n, T, S> and poly_evaluator<n, T, HL, S> are specialized for n == 1, as in that case, there is no lower level, and hence these specializations provide no evaluate template and try to pass evaluation further down, but right away pass it higher up to their owner, or in case of poly_evaluator<1, T, S>, it will just do the evaluation.
this concludes the fourth part. the last part will concentrate something less complicated: implementing long division and the euclidean algorithm.
comments.