基于表达式模版（expression template）的惰性求值（lazy evaluation）

读代码碰到惰性求值的模型编程技术，动手实践下。

• 目标：设计一个支持任意四则表达式运算的向量类（Vec）

定义一个类Vec

template<typename T>struct Vec {        Vec() : len(0), dptr(0), owned(false) {}    Vec(T* dptr, int len) : len(len), dptr(dptr), owned(false) {}    Vec(int len) : len(len), owned(true) {        dptr = new T[len];    }    ~Vec() {        if (owned) delete[] dptr;           dptr = 0;    }    inline int length() {        return len;    }    inline T operator[] (int i) const {        return dptr[i];    }    inline T & operator[] (int i) {        return dptr[i];    }    inline Vec & operator= (const Vec<T>& src) {        #pragma omp parallel for        for (int i = 0; i < len; i++) {            dptr[i] = src[i];        }        return *this;    }private:    int len;    T* dptr;    bool owned;};

支持任意长度四则运算，如

Vec<double> a, b, c, d, e, f;a = b + c - d / e * f;

• first try
  重定义运算符+、-、*、/。这是C++的一个常规解决方案

template<T>Vec<T> & operator+ (Vec<T> &lhs, Vec<T> &rhs) {    // assert(lhs.size() == rhs.size());    Vec<T> rs(lhs.size());    for (int i = 0; i < lhs.size(); i++) {        rs[i] = lhs[i] + rhs[i];    }    return rs;}

这一方案的问题在于，运算过程中需要分配临时空间。此外，还存在多次函数调用。

• 更好的方案——表达式模型
  这里， 我们使用表达式模版实现运算的惰性求值。不仅不需要额外的空间，也减少函数调用开销。

• 表达式模版
  表达式（等号右边部分）可以用一个表达式树抽象的表示。其中，叶子节点（终结符）是我们的向量，它也是一种特殊的表达式。树的根节点是运算符，左右子树是子表达式。

• 实现
  – 首先，我们定义表达式。它是所有可能表达式的父类。

template<typename RealType>struct Exp {    inline const RealType& self() const {        return *static_cast<const RealType*>(this);    }};

实质上，它只是一个wrapper。它的作用是，当我们需要将一个对象做为表达式传递是时，它将其他封装。在传递之后，通过self()函数再得到原来的对象。

例如，我们如下定义Vec：

template<T>struct Vec: Exp<Vec<T>> {    ... ..}

对比常规定义：

template<T>struct Vec {    ... ...}

Vec的完整定义如下：

template<typename T>struct Vec : public Exp < Vec<T> > {    typedef T value_type;    int len;    T* dptr;    bool owned;    Vec(T* dptr, int len) : len(len), dptr(dptr), owned(false) {}    Vec(int len) : len(len), owned(true) {        dptr = new T[len];    }    ~Vec() {        if (owned) delete[] dptr;           dptr = 0;    }    inline T operator[] (int i) const {        return dptr[i];    }    template<typename EType>    inline Vec & operator= (const Exp<EType>& src_) {        const EType &src = src_.self();#pragma omp parallel for        for (int i = 0; i < len; i++) {            dptr[i] = src[i];        }        return *this;    }};

唯一需要解释是赋值操作

template<typename EType>    inline Vec & operator= (const Exp<EType>& src_) {        const EType &src = src_.self();#pragma omp parallel for        for (int i = 0; i < len; i++) {            dptr[i] = src[i];        }        return *this;    }

Vec接受一个表达式，表达式必须提供operator[]函数，返回相应的值。正是由于[]的定义，使得惰性求值成为可能。

以上，我们已经有了叶子节点（Vec）。要构造表达式树，我们要定义每个中间节点和根节点。它们本质上是二元操作。

template<typename Op, typename TLhs, typename TRhs>struct BinaryOpExp : Exp < BinaryOpExp<Op, TLhs, TRhs> > {    const TLhs &lhs;    const TRhs &rhs;    typedef typename ReturnType<TLhs, TRhs>::value_type value_type;    BinaryOpExp(const TLhs &lhs, const TRhs &rhs) : lhs(lhs), rhs(rhs) {}    inline value_type operator[] (int i) const {        return Op::eval(lhs[i], rhs[i]);    }};

其中，ReturnType 只是一个简单的功能模版。

template<typename TLhs, typename TRhs>struct ReturnType {    typedef typename TLhs::value_type value_type;};

作为表达式，BinaryOpExp 重载了我们需要的 operator[]。

最后要做的是，重载+号等运算符

template<typename T>struct add {    inline static T eval(const T& lhs, const T& rhs) {        return lhs + rhs;    }};template<typename TLhs, typename TRhs>inline BinaryOpExp<add<typename ReturnType<TLhs, TRhs>::value_type>, TLhs, TRhs>operator+ (const Exp<TLhs> &lhs, const Exp<TRhs> &rhs) {    return BinaryOpExp<detail::add<typename ReturnType<TLhs, TRhs>::value_type>, TLhs, TRhs>(lhs.self(), rhs.self());}

一个简单的测试：

int main() {    const int n = 3;    double sa[n] = { 1, 2, 3 };    double sb[n] = { 2, 3, 4 };    double sc[n] = { 3, 4, 5 };    double sd[n] = { 4, 5, 6 };    double se[n] = { 5, 6, 7 };    double sf[n] = { 6, 7, 8 };    Vec<double> A(sa, n), B(sb, n), C(sc, n), D(sd, n), E(se, n), F(sf, n);    // run expression, this expression is longer:)    A = B + C - D * E / F;    for (int i = 0; i < n; ++i) {        printf("%d:%f == %f + %f - %f * %f / %f == %f\n", i,            A[i], B[i], C[i], D[i], E[i], F[i], B[i] + C[i] - D[i] * E[i] / F[i]);    }    return 0;}

输出结果：

0:1.666667 == 2.000000 + 3.000000 - 4.000000 * 5.000000 / 6.000000 == 1.6666671:2.714286 == 3.000000 + 4.000000 - 5.000000 * 6.000000 / 7.000000 == 2.7142862:3.750000 == 4.000000 + 5.000000 - 6.000000 * 7.000000 / 8.000000 == 3.750000

除基本的+、-、*、\之外，我们还可以自定义二元运算符。

template<typename Op, typename TLhs, typename TRhs>inline BinaryOpExp<Op, TLhs, TRhs> F(const Exp<TLhs> &lhs, const Exp<TRhs> &rhs) {    return BinaryOpExp<Op, TLhs, TRhs>(lhs.self(), rhs.self());}

类似的，我们可以定义一元操作运算

template<typename Op, typename T>struct UnaryOpExp : Exp < UnaryOpExp<Op, T> > {    const T &arg;    typedef typename T::value_type value_type;    UnaryOpExp(const T &arg) : arg(arg) {}    inline value_type operator[] (int i) const {        return Op::eval(arg[i]);    }};template<typename Op, typename T>inline UnaryOpExp<Op, T> F(const Exp<T> &arg) {    return UnaryOpExp<Op, T>(arg.self());}

我们重载sin函数

template<typename T>struct sinOp {    inline static T eval(const T& arg) {        return std::sin(arg);    }};template<typename T>UnaryOpExp<detail::sinOp<typename T::value_type>, T> sin(const Exp<T> &arg) {    return UnaryOpExp<detail::sinOp<typename T::value_type>, T>(arg.self());}

一个简单的测试：

int main() {    const int n = 3;    double sa[n] = { 1, 2, 3 };    double sb[n] = { 2, 3, 4 };    double sc[n] = { 3, 4, 5 };    Vec<double> A(sa, n), B(sb, n), C(sc, n);    A = sin(B) + sin(C);    for (int i = 0; i < n; ++i) {        printf("%d:%f == sin(%f) + sin(%f) == %f\n", i, A[i], B[i], C[i], sin(B[i]) + sin(C[i]));    }       return 0;}

输出结果如下：

0:1.050417 == sin(2.000000) + sin(3.000000) == 1.0504171:-0.615682 == sin(3.000000) + sin(4.000000) == -0.6156822:-1.715727 == sin(4.000000) + sin(5.000000) == -1.715727