Currently we define polynomials over doubles which are approximations of real numbers from mathematics. However, from the mathematical point of view the underlying structure for a polynomial is a ring. The polynomial requires the operations "*" and "+" to be defined for the underlying type. Real numbers fulfill the axioms of a ring but much more general objects do likewise. So it would be nice if our code for polynomials worked for arbitrary types with a "*" and "+" defined. Exactly this is possible with "templates" which are parametric types. What is a template or parametric type? Consider the following function declaration:
1 double f(double x);It declares a function getting a double and returning a double. x is called a variable because it is variable. The type double however is not variable. Templates are variable types. Defining a grammar which allows to define variable types is not really trivial (think about it). In C++ the parametric type is declared before its actual usage in a variable declaration:
1 template <class T> T f(T x);This tells the compiler: Attention, T is going to be a parametric type. So this line reads: There is a function f getting any type of argument returning the same type of argument. If f is now called say with a string it will return a string. Nevertheless, the whole template thing is much more flexible than can be described in this short tutorial. Analogously to the template function declaration just shown we may declare a template class. We shall apply this to our Polynomial class:
1 #include <iostream> // use the iostream library (for input and output)
2 #include <vector> // use the STL type vector (a container type)
3
4 using namespace std; // add the namespace "std" to standard scope resolution
5
6 //-------------------------------------------------------------
7 // declaration:
8 template <class T>
9 class Polynomial : public vector<T> { // declare a new type "Polynomial" which inheritas from a vector of Ts
10 public: // the following is publicly visibel
11 // constructor
12 Polynomial(int n); // gets an integer as degree
13
14 // function evaluation
15 T operator () (T x) const;
16
17 // declare addition operator
18 Polynomial operator + (const Polynomial & p) const;
19
20 }; // the class declaration of "Polynomial" ends here
21
22 // define output operator
23 template <class T>
24 ostream & operator << (ostream & o, const Polynomial<T> & p);
25
26 //-------------------------------------------------------------
27 // implementation:
28 template <class T>
29 Polynomial<T>::Polynomial(int n) :
30 vector<T>(n+1, 0) // allocate space for the coefficients and initialize them to 0
31 {} // nothing else to do
32
33 template <class T>
34 T Polynomial<T>::operator () (T x) const
35 {
36 T xx = 1; // temporary for powers of x
37 T y = 0; // will contain the result
38 for ( unsigned int i=0 ; i<this->size() ; ++i )
39 {
40 y += (*this)[i] * xx; // add coefficient * x^n
41 xx *= x; // next power
42 }
43 return y; // return result
44 }
45
46
47 template <class T>
48 Polynomial<T> Polynomial<T>::operator + (const Polynomial<T> & p) const
49 {
50 Polynomial<T> r(*this); // copy left polynomial into result
51 for ( unsigned int i=0 ; i<p.size() ; ++i ) // loop over right polynomial
52 if ( i<this->size() ) // check if coefficient index exist in result polysnomial
53 r[i] += p[i]; // yes: add coefficients from right polynomial
54 else
55 r.push_back(p[i]); // no: append coefficients from right polynomial
56 return r; // return result
57 }
58
59 template <class T>
60 ostream & operator << (ostream & o, const Polynomial<T> & p)
61 {
62 for ( int i=p.size()-1 ; i>=0 ; i-- )
63 {
64 if ( p[i]!=0 )
65 {
66 if ( i!=p.size()-1 )
67 o << " + ";
68 o << p[i];
69 if ( i>0 )
70 o << "*x^" << i;
71 }
72 }
73 return o;
74 }
75
76
77 int main()
78 {
79 Polynomial<double> p(4); // create a polynomial of Ts of 4th degree
80 p[0] = 1; // 1*x^0
81 p[4] = 2; // 2*x^4
82
83 Polynomial<double> q(2); // create a polynomial of Ts of 2nd degree
84 q[0] = 1; // 1*x^0
85 q[1] = 2; // 2*x^1
86 q[2] = 3; // 3*x^2
87
88 cout << "p = " << p << endl; // print it
89 cout << "q = " << q << endl; // print it
90 cout << "p(3.4) = " << p(3.4) << endl; // evaluate and print it
91 cout << "q(3.4) = " << q(3.4) << endl; // evaluate and print it
92 cout << "p(3)+q(3) = " << p(3)+q(3) << endl;
93 cout << "p+q = " << p+q << endl;
94 cout << "(p+q)(3) = " << (p+q)(3) << endl;
95 return 0;
96 }
Download Full Source
This Polynomial will work exactly as before but is now pretty general. Consequently, we may do things like the following:
1 Polynomial<Matrix<double> > p(3); 2 Matrix<Polynomial<double> > q(3);The first line generates a polynomial of matrices of doubles. The second line generates a matrix of polynomials of doubles. This is expressive power! There is no runtime penalty. Think about what you would have to do to achieve this with FORTRAN! You might say: "I never use polynomials of matrices." You may be right. But if you write a quantum chemistry program package you could certainly use something like the following:
1 Basis<SlaterDeterminant<SpinOrbital> >