T. Veldhuizen, "Expression Templates," C++ Report, Vol. 7 No. 5 June 1995, pp. 26-31

C++ Report (ISSN 1040-6042) is published nine times per year by SIGS Publications Inc., 71 West 23rd St., 3rd Floor, New York, NY 10010. © Copyright 1995 SIGS Publications Inc.

Back to the Blitz++ Home Page
Back to Todd's Home Page

Introduction

The technique of expression templates allows expressions to be passed to functions as an argument and inlined into the function body. Here are some examples:

          // Integrate a function from 0 to 10
          DoublePlaceholder x;
          double result = integrate(x/(1.0+x), 0.0, 10.0);
          
          // Count the elements between 0 and 100 in a
          collection of int
          List<int> myCollection;
          IntPlaceholder y;
          int n = myCollection.count(y >= 0 && y <= 100);
          
          // Fill a vector with a sampled sine function
          Vector<double> myVec(100);
          Vector<double>::placeholderType i;
          myVec.fill(sin(2*M_PI*i/100));  // eqn in terms of element #

          Vector y(100), a(100), b(100), c(100),
          d(100);
          y = (a+b)/(c-d);
          

So how does it work? The trick is that the expression is parsed at compile time, and stored as nested template arguments of an "expression type. Let's work through a simplified version of the first example, which passes a mathematical function to an integration routine. We'll invoke a function evaluate() which will evaluate f(x) at a range of points:

          
          int main()
          {
              DExpr<DExprIdentity> x;         // Placeholder
              evaluate(x/(1.0+x), 0.0, 10.0);
              return 0;
          }

Figure 1. (A) Parse tree of "1.0+x", showing the similarity to the expression type generated by the compiler (B). The type DApAdd is an applicative template class which represents an addition operation. In (C), the type name is shown in a non-expanded view. For this implementation, type names are similar to a postfix traversal of the expression tree.
To write a function which accepts an expression as an argument, we include a template parameter for the expression type. Here is the code for the evaluate() routine. The parenthesis operator of expr is overloaded to evaluate the expression:

template<class Expr>
void evaluate(DExpr<Expr> expr, double start, double end)
{
    const double step = 1.0;

    for (double i=start; i < end; i += step)
        cout << expr(i) << endl;
}

              cout << i/(1.0+i) >> endl;

              evaluate(x/(1.0+x), 0.0, 10.0);

Figure 2. Template instance generated by x/(1.0+x): DExpr<DBinExprOp< DExpr <DExprIdentity>, DExpr<DBinExprOp<DExpr<DExprLiteral>, DExpr<DExprIdentity>, DApAdd>>, DApDivide>>
As an example of how the return types of operators cause the compiler to infer expression types, here is an operator/() which produces a type representing the division of two subexpressions DExpr<A> and DExpr<B>:

          template<class A, class B>
          DExpr<DBinExprOp<DExpr<A>, DExpr<B>, DApDivide> >
          operator/(const DExpr<A>& a, const DExpr<B>& b)
          {
              typedef DBinExprOp<DExpr<A>, DExpr<B>,DApDivide> ExprT;
              return DExpr<ExprT>(ExprT(a,b));
          }

The typedef ExprT is the type DBinExprOp< ... , ... , DApDivide> which is to be wrapped by a DExpr<...>. Expression classes take constructor arguments of the embedded types; for example, DExpr<A> takes a const A& as a constructor argument. So DExpr<ExprT> takes a constructor argument of type const ExprT&. These arguments are needed so that instance data (such as literals which appear in the expression) are preserved.

The idea of applicative template classes has been borrowed from the Standard Template Library (STL) developed by Alexander Stepanov and Meng Lee, of Hewlett Packard Laboratories [2]. The STL has been accepted as part of the ISO/ANSI Standard C++ Library. In the STL, an applicative template class provides an inline operator() which applies an operation to its arguments and returns the result. For expression templates, the application function can (and should) be a static member function. Since operator() cannot be declared static, an apply() member function is used instead:

          // DApDivide -- divide two doubles
          class DApDivide {
          public:
              static inline double apply(double a, double b)
              { return a/b; }
          };

          template<class A, class B, class Op>
          class DBinExprOp {
              A a_;
              B b_;
          
          public:
              DBinExprOp(const A& a, const B& b)
                : a_(a), b_(b)
              { }
          
              double operator()(double x) const
              { return Op::apply(a_(x), b_(x)); }
          };

          double mean=5.0, sigma=2.0;   // mathematical
          constants
          DExpr<DExprIdentity> x;
          evaluate( 1.0/(sqrt(2*M_PI)*sigma) * exp(sqr(x-mean)/
              (-2*sigma*sigma)), 0.0, 10.0);

          double __t1 = 1.0/(sqrt(2*M_PI)*sigma);
          double __t2 = (-2*sigma*sigma);
          evaluate( __t1 * exp(sqr(x-mean)/__t2), 0.0,
          10.0);

Figure 3 shows a benchmark of a simple integration function, comparing callback functions and expression templates to manually inlined code. Expression template versions ran at 95% the efficiency of hand-coded C. With larger expressions, the advantage over pointers-to-functions shrinks. This advantage is also diminished if the function using the expression (in this example, integrate()) does significant computations other than evaluating the expression.

It is possible to generalize the classes presented here to expressions involving arbitrary types (instead of just doubles). Expression templates can also be used in situations where multiple variables are needed -- such as filling out a matrix according to an equation. Since most compilers do not yet handle member function templates, a virtual function trick has to be used to pass expressions to member functions of a class. This technique is illustrated in the next example, which solves a significant outstanding problem in the design of matrix and vector class libraries.

Optimizing Vector Expressions

          DoubleVec y(1000), a(1000), b(1000), c(1000),
          d(1000);
          y = (a+b)/(c-d);

          DoubleVec __t1 = a+b;
          DoubleVec __t2 = c-d;
          DoubleVec __t3 = __t1/__t2;
          y = __t3;

          // double *yp, *ap, *bp, *cp, *dp all point into the vectors.
          // double *yend is the end of the y vector.
          do {
              *yp = (*ap + *bp)/(*cp - *dp);
              ++ap; ++bp; ++cp; ++dp;
          } while (++yp != yend);
          

          class DVec {
          
          private:
              double* data_;
              int length_;
          
          public:
              typedef double* iterT;
          
              DVec(int n)
                  : length_(n)
              { data_ = new double[n]; }
          
              ~DVec()
              { delete [] data_; }
          
              iterT begin() const
              { return data_; }
          
              iterT end() const
              { return data_ + length_; }
          
              template<class A>
              DVec& operator=(DVExpr<A>);
          };

          template<class A>
          class DVExpr {
          
          private:
              A iter_;
          
          public:
              DVExpr(const A& a)
                  : iter_(a)
              { }
          
              double operator*() const
              { return *iter_; }
          
              void operator++()
              { ++iter_; }
          };

          template<class A>
          DVec& DVec::operator=(DVExpr<A> result)
          {
              // Get a beginning and end iterator
          for the vector
              iterT iter = begin(), endIter = end();
          
              // Store the result in the vector
              do {
                  *iter = *result;   // Inlined expression
                  ++result;
              } while (++iter != endIter);
          
              return *this;
          }

As before, operators build expression types using nested template arguments. Several versions of the each operator are required to handle the combinations in which DVec, DVExpr<A>, and numeric constants can appear in expressions. For a simple example of how expressions are built, consider the code snippet

          DVec a(10), b(10);
          y = a+b;
          

DVExpr<DVBinExprOp<DVec::iterT,DVec::iterT,DApAdd> >
operator+(const DVec& a, const DVec& b)
{
    typedef DVBinExprOp<DVec::iterT,DVec::iterT,DApAdd> ExprT;
    return DVExpr<ExprT>(ExprT(a.begin(),b.begin()));
}

Figure 4. Type name for the expression (a+b)/(c-d): DVExpr<DVBinExprOp<DVExpr< DVBinExprOp<DVec::iterT, DVec::iterT, DApAdd>>, DVExpr<DVBinExprOp< DVec::iterT, DVec::iterT, DApSubtract>>, DApDivide>>

Benchmark Results

How to reach me: