数值计算优化方法C/C++(七)——奇异递归模板

奇异递归模板简介

奇异递归模板是模板的一种使用技巧，通常这种技巧和一种叫做静态多态特性一起出现。是一种继承时将子类型作为模板参数传给基类的一种模板使用方法。

奇异递归模板应用

对象计数

有时我们可能会需要对一些类型的对象计数，通常我们是通过一个static int成员来实现计数功能。这个实现本身并不复杂，但是如果我们想要实现多个类型都能够计数，我们就需要在每个类型中都添加一个static int成员，同时在它们的构造函数中添加一个计数操作。如果需要计数的类型比较多就可能会有较多的重复性操作。那么我们自然希望能通过某种方式来把计数的功能单独分离出来，形成一个可以复用的模块。面向对象的思想中最容易想到的就是继承，做一个基类实现计数功能，所有需要计数的类型只需要继承它就可以。直觉上似乎并没有问题，但是仔细思考就会发现由于计数本身是利用的静态成员，这样所有继承都会公用一个计数，多个子类的构造都会在记在同一个变量中，而无法做到每个子类单独计数。利用奇异递归模板可以很容易实现不同子类拥有各自的静态计数成员。代码如下

template <typename T>
struct Count {
  static size_t count;
  Count() {count++;}
  ~Count() {count--;}
};

struct A : public Count<A> {
  A() {cout << Count<A>::count << "\n";}
  ~A() {cout << Count<A>::count << '\n';}
};

template<>
size_t Count<A>::count = 0;

struct B : public Count<B> {
  B() {cout << Count<B>::count << "\n";}
  ~B() {cout << Count<B>::count << '\n';}
};

template<>
size_t Count<B>::count = 0;
int main() {
  A *p[5];
  B *q[5];
  for (int i=0;i<5;i++) {
    p[i] = new A;
    q[i] = new B;
  }
  for (int i=0;i<5;i++) {
    delete p[i];
    delete q[i];
  }
  return 0;
}

可以看到继承Count的子类A和B可以实现单独计数，互相不会干扰。这就是奇异递归模板的一个作用，利用子类作为模板参数传给基类，从而让基类能够相互独立。在我们想在基类中设置静态成员实现一些功能，却又希望不同子类之间不会因为使用基类中的静态成员而互相干扰时，奇异递归模板就是一个很好的选择。

静态多态

静态多态就是一个比较神奇的操作了，如下代码

template <typename T>
struct Base {
  void show() {
    static_cast<T*>(this)->show();
  }
};
struct A : public Base<A> {
  void show() {
    cout << "this is A\n";
  }
};
struct B : public Base<B> {
  void show() {
    cout << "this is B\n";
  }
};
int main() {
  A a;
  B b;
  a.show();
  b.show();
}

首先这个例子并不是多态，尽管有继承，有奇异递归模板，但是这确实不是多态。a，b的不同表现完全来自于a，b本身就是不同的类型。所以我们需要利用模板本身就具有的一种表达接口多态形式的特性，配合模板我们可以实现一种接近真正多态的表达形式。
在上述实现的基础上我们添加如下函数

template <typename T>
void show(Base<T> *in) {
  in->show();
}
int main() {
  A a;
  B b;
  show(a);
  show(b);
}

这样是不是就有点像真正的多态了。又有基类指针，又有继承。这回就是真的调用的基类的show接口，然后实现了不同表现形式。不过其实往深了考虑一层就会发现这个多态其实依然不是真正的多态，因为这次虽然调用了基类的同名接口实现了子类的不同行为，但是前面计数的实现我们说了，这两个基类本就不同，而且所谓多态的实现是在基类中做了一次指针强转，把基类指针强转成了子类指针，本质上这依旧是重写而不是多态。但毕竟我们这个特性叫做静态多态，是一种编译期的特性和真正的运行期多态肯定不可能完全相同，二者处理的问题也不一样。不过这种方法确实在一定程度上可以像真正的多态一样，保持相同的接口，却可以使用不同的实现。下面介绍一种编译期多态的真正应用。

表达式模板

在之前的文章中我介绍过表达式模板这种元编程技术，它利用模板的特性可以实现延迟计算，表达式展开的功能，使得vector的四则运算可以节省一个临时vector对象的构造。但是那个文章里面表达式模板的实现及其复杂，需要大量特化模板，而且不易扩展，我当时写那篇文章时其实是测试了加减乘以及数乘四种运算的，但是考虑到文章篇幅(其实就是懒)，我并没有贴上全部代码，因为那确实太长了。这次我们有了奇异递归模板，我们就可以把那时候实现中的一些公共操作抽出来，作为基类，利用继承来实现代码复用，从而使得扩展变得容易。
首先我们需要做一个描述运算的基类

template <typename DT, typename Ch>
struct BaseOp {
  typedef DT Type;
  Type operator[](size_t i) {
    return (*static_cast<Ch*>(this))[i];
  }
};

然后通过这个基类派生出不同类型的运算，这里就做一个简单矢量加法和置零的例子，所以需要定义一元运算和二元运算

template <typename T1, typename Ch, typename DT = typename Trait<T1>::Type>
struct UnaryOp : public BaseOp<DT, UnaryOp<T1, Ch>> {
  typedef DT Type;
  Type op(Type &x) {
    return static_cast<Ch*>(this)->op(x);
  }
  Type operator[](size_t i) {
    return op(operand1[i]);
  }
  UnaryOp(const T1 &t1):operand1(t1) {}
  private:
  T1 operand1;
};
template <typename T1, typename T2, typename Ch,
          typename F = typename enable_if<
          is_same<typename Trait<T1>::Type,
                  typename Trait<T2>::Type>::value,
                  void>::type>
struct BinaryOp : public BaseOp<typename Trait<T1>::Type, BinaryOp<T1, T2, Ch>> {
  typedef typename Trait<T1>::Type Type;
  Type op(Type &x, Type y) {
    return static_cast<Ch*>(this)->op(x, y);
  }
  Type operator[](size_t i) {
    return op(operand1[i], operand2[i]);
  }
  BinaryOp(const T1& t1, const T2& t2):
    operand1(t1), operand2(t2) {}
  private:
  T1 operand1;
  T2 operand2;
};

然后我们还需要一个基础的用于计算的数据类型

template <typename T>
struct Vector : public BaseOp<typename Trait<T>::Type, Vector<T>> {
  typedef typename Trait<T>::Type Type;
  Type& operator[](size_t i) {
    return data[i];
  }
  Vector(T& t):data(t){}
  private:
  T& data;
};

这样我们基础设施就完善了，接下来就可以加真正的计算函数了(这里已经不是c++中的那个函数了，而是一系列计算用的模板类)

// 置零操作
template <typename T>
struct Zero : public UnaryOp<T, Zero<T>> {
  typedef typename T::Type Type;
  Type op(Type &x) {
    return (x = 0);
  }
  Zero(T& t):UnaryOp<T, Zero<T>>(t) {};
};
// 矢量加
template <typename T1, typename T2>
struct Add : public BinaryOp<T1, T2, Add<T1, T2>> {
  typedef typename Trait<T1>::Type Type;
  Type op(Type x, Type y) {
    return x+y;
  }
  Add(T1 &t1, T2 &t2):BinaryOp<T1, T2, Add<T1, T2>>(t1, t2) {}
};
// 赋值操作
template <typename T1, typename T2>
struct Assign : public BinaryOp<T1, T2, Assign<T1, T2>> {
  typedef typename Trait<T1>::Type Type;
  Type op(Type &x, Type y) {
    return x = y;
  }
  Assign(T1 &t1, T2 &t2):BinaryOp<T1, T2, Assign<T1, T2>>(t1, t2) {}
};
// 表达式模板的作用是延迟计算，Calculate用于表示该数据需要使用了，进行真正的计算
template <typename T>
struct Calculate : public BaseOp<typename Trait<T>::Type, Calculate<T>> {
  typedef typename Trait<T>::Type Type;
  Type operator[](size_t i) {
    return operand[i];
  }
  Calculate(T& t, size_t len):operand(t) {
    for (int i = 0; i < len; i++) (*this)[i];
  }
  private:
  T operand;
};

尽管这样就加好了所有运算操作，但是这显然不是一个可以用的接口，因此还要封装成一个函数得到一个易用的接口

template <typename T>
Zero<T> zero(T& x) {
  return Zero<T>(x);
}

template <typename T1, typename T2>
Add<T1, T2> add(T1 &t1, T2 &t2) {
  return Add<T1,T2>(t1,t2);
}

template <typename T1, typename T2>
Assign<T1, T2> assign(T1 &t1, T2 &t2) {
  return Assign<T1,T2>(t1,t2);
}

template <typename T>
Calculate<T> calculate(T& x, size_t len) {
  return Calculate<T>(x, len);
}

测试一下

int main(){
  vector<int> x({1, 2, 3});
  Vector<vector<int>> a(x);
  cout << a[0] << '\t' << a[1] << '\t' << a[2] << '\n';
  auto t1 = add(a, a); 
  auto t12 = assign(a, t1);
  auto r1 = calculate(t12, x.size());
  cout << a[0] << '\t' << a[1] << '\t' << a[2] << '\n';
  auto t2 = zero(a);
  auto r2 = calculate(t2, x.size());
  cout << a[0] << '\t' << a[1] << '\t' << a[2] << '\n';
}

计算结果如下

1       2       3
2       4       6
0       0       0

我们测试一下这个计算和通过一个接收C++原生vector的函数来完成加法的操作的性能比较。

std::vector<int> origin_add(const std::vector<int> &a, const std::vector<int> &b) {
	std::vector<int> ret(a.size());
	for (int i = 0; i < ret.size(); i++)
		ret[i] = a[i] + b[i];
	return ret;
}
int main(){
  int N = 500000;
  vector<int> x(N, 1), y(N, 0);
  Vector<vector<int>> a(x);
  Vector<vector<int>> b(y);

  double time = 0;
  clock_t s = clock();
  auto op1 = add(a, a); 
  auto op2 = assign(b, op1);
  auto ret = calculate(op2, x.size());
  time = static_cast<double>(clock() - s)/CLOCKS_PER_SEC;
  cout << y[5000] << "\t|\ttime:\t" << time <<"\n";

  s = clock();
  y = origin_add(x, x);
  time = static_cast<double>(clock() - s)/CLOCKS_PER_SEC;
  cout << y[5000] << "\t|\ttime:\t" << time <<"\n";
}

使用g++的O3优化结果为

2       |       time:   0.000306
2       |       time:   0.001324

可以看到使用奇异递归模板时性能更好。