skip to main content.

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

by felix.

welcome to the third part of the multivariate polynomials in c++ using templates series. this part will explain how to compute partial derivatives using partial template specialization tricks.

the idea is very simple. given f \in R[X_0, \dots, X_{n-1}], written as f = \sum_{i=0}^k a_i X_0^i with a_i \in R[X_1, \dots, X_{n-1}], we have

\frac{\partial f}{\partial X_0} = \sum_{i=1}^k a_i X_0^{i-1}

and

\frac{\partial f}{\partial X_j} = \sum_{i=0}^k \biggl( \frac{\partial}{\partial X_j} a_i \biggr) X_0^i \quad \text{for } j > 0.

this allows us to implement partial derivatives iteratively. we assume that the number of the indeterminate by which one should differentiate is already known during compile time. the implementation will be done by a templated function

1template<int N, int n, class T>
2inline poly<n, T> derivative(const poly<n, T> & p);

the purpose of the template parameter N is to say that \frac{\partial p}{\partial X_N} should be computed. given a polynomial p, the function can be called simply by derivative<N>(p). to implement derivative<N>(p), we would like to again use partial template specialization to distinguish between the cases N == 0 and N > 0. unfortunately, partial template specialization is not available for template functions. hence, we need to introduce a helper template class:

1template<int N, int n, class T>
2class derivative_computer;

it should contain one static function, static inline void computeDerivative(const poly<n, T> & src, poly<n, T> & dest); which allows to compute the partial derivative \frac{\partial}{\partial X_N} of the polynomial in src and store it into dest. then derviative<N>(p) can simply be implemented by

1template<int N, int n, class T>
2inline poly<n, T, std::allocator<T> > derivative(const poly<n, T, std::allocator<T> > & p)
3{
4    poly<n, T> res;
5    derivative_computer<N, n, T>::computeDerivative(p, res);
6    return res;
7}

now let us focus on the class template derivative_computer<N, n, T>. the general implementation, for N > 0, uses the second formula from above and iteratively calls derivative_computer<N-1, n-1, T>::computeDerivative():

 1template<int N, int n, class T>
 2class derivative_computer
 3{
 4public:
 5    static inline void computeDerivative(const poly<n, T> & src, poly<n, T> & dest)
 6    {
 7        dest.d_value.resize(src.degree() + 1);
 8        for (int i = 0; i <= src.degree(); ++i)
 9            derivative_computer<N - 1, n - 1, T>::computeDerivative(src[i], dest[i]);
10        dest.normalize();
11    }
12};

now we are left with the case that N = 0. the implementation here is also simple:

 1template<int n, class T>
 2class derivative_computer<0, n, T>
 3{
 4public:
 5    static inline void computeDerivative(const poly<n, T> & src, poly<n, T> & dest)
 6    {
 7        dest.d_value.resize(src.degree() >= 0 ? src.degree() : 0);
 8        for (int i = 1; i <= src.degree(); ++i)
 9            dest[i - 1] = src[i] * (T)i;
10        dest.normalize();
11    }
12};

here, we cast i to an object of type T before multiplication. from this point on, no further iteration is needed.

we are left with the case that some (malicious or just confused) user calls derivative<N>(p) for a polynomial p of type poly<n, T> with n <= N. in that case, the compiler tries to instanciate derivative_computer<k, 0, T>::computeDerivative() and tries to treat poly<0, T> as a polynomial. this will result in strange error messages, and it would be nice if there is something more readable. for that reason, i added another two specializations:

 1template<class T>
 2class derivative_computer<0, 0, T>
 3{
 4public:
 5    class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
 6    
 7    static inline void computeDerivative(const poly<0, T> &, poly<0, T> &)
 8    {
 9        ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
10    }
11};
12
13template<int N, class T>
14class derivative_computer<N, 0, T>
15{
16public:
17    class ERROR_N_must_be_strictly_less_than_n_in_derivative_template;
18    
19    static inline void computeDerivative(const poly<0, T> &, poly<0, T> &)
20    {
21        ERROR_N_must_be_strictly_less_than_n_in_derivative_template();
22    }
23};

these two definitions compile fine, but as soon as one of them is instanciated by the compiler, it will try to find out something about the class ERROR_N_must_be_strictly_less_than_n_in_derivative_template and not succeed, and spit out some error message hopefully containing the string ERROR_N_must_be_strictly_less_than_n_in_derivative_template. for example, if i replace the line computing the partial derivative of the intersection equation in poly-test.cpp by std::cout << derivative<2>(l) << "\n";, g++ will complain as follows:

1In file included from poly-test.cpp:2:0:
2polynomial.hpp: In static member function ‘static void derivative_computer<N, 0, T>::computeDerivative(const poly<0, T>&, poly<0, T>&) [with int N = 1, T = double]’:
3polynomial.hpp:1200:13:   instantiated from ‘static void derivative_computer<N, n, T>::computeDerivative(const poly<n, T>&, poly<n, T>&) [with int N = 2, int n = 1, T = double]’
4polynomial.hpp:1247:5:   instantiated from ‘poly<n, T> derivative(const poly<n, T>&) [with int N = 2, int n = 1, T = double]’
5poly-test.cpp:19:33:   instantiated from here
6polynomial.hpp:1227:9: error: invalid use of incomplete type ‘struct derivative_computer<1, 0, double>::ERROR_N_must_be_strictly_less_than_n_in_derivative_template’
7polynomial.hpp:1223:11: error: declaration of ‘struct derivative_computer<1, 0, double>::ERROR_N_must_be_strictly_less_than_n_in_derivative_template’
8make: *** [poly-test.o] Error 1

this is not exactly readable, but the uppercase ERROR quickly lets you look at the important part of that message. (note that i removed the allocators from the error messages, so that it is in sync with the above writing. if you try it out yourself with the source code, you'll see a slightly more complicated error message.)

this completes part three of the series. in the next part, i will discuss how to evaluate polynomials efficiently using templates.

posted in: computer
tags:
places:

comments.

no comments.