利用模板元编程实现解循环优化

利用模板元编程实现解循环优化
Unroll Loops Optimization with Template Meta Programming

作者：王鹏

简介

在《C++ Templates: The Complete Guide》一书中（以下简称书），提出了模板元编程最早的实际应用之一：在数值运算中进行解循环优化。
而本文的标题是噱头！本文的真正目的是指出这种优化措施在增加复杂性的同时，并不一定能明显改善效率。应当谨慎使用该技术——默认不使用该技术，在点积计算确实是效率瓶颈时考虑采用该技术，并认真测试该技术是否真能提高效率。

背景

数值运算库中常常需要提供向量点积（dot_product）运算。其定义用C++代码描述也许更清楚～

templateT dot_product(int dim,const T v1[],const T v2[]) {T result=0;for (int i=0;i我们可以使用这个函数，求2个向量的点积
int v1[] = {1,2,3};int v2[] = {4,5,6};int r1 = dot_product(3,v1,v2);
得到r1=32

书中指出：“这个结果是正确的，但是在要求高效率的应用中，它耗费了太多的时间”，对于这类特殊的问题，“简单的把循环展开”，如：r1=v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]
“反而会好得多”。
如何便捷的展开循环？将每个dot_product(dim,...) 手工重写成展开式？还是设计一整套函数: dot_product_dim_2，dot_product_dim_3，... ？
无疑这是一个冗长乏味易出错的方案。
书中提出了一种使用模板元编程解决该问题的方案，让我们来看看他是怎么做的。
模板元编程——解循环
// 首先是递归模板templatestruct DotProduct {static T execute(const T v1[],const T v2[]);}templateT DotProduct::execute(const T v1[],const T v2[]) {return v1[0]*v2[0] + DotProduct::execute(v1+1,v2+1);}// 递归边界模板templatestruct DotProduct<1,T> {static T execute(const T v1[],const T v2[]);}templateT DotProduct<1,T>::execute(const T v1[],const T v2[]) {return v1[0]*v2[0];}
我们可以使用DotProduct 来进行点积计算了：
int v1[] = {1,2,3}; int v2[] = {4,5,6};int r2 = DotProduct<3,int>::execute(v1,v2);

计算r2的函数，在DotProduct<3,int>::execute 被实例化时就确定了。
编译器将会实例化
int DotProduct<3,int>::execute(const int v1[],const int v2[]) {return v1[0]*v2[0] + DotProduct<2,int>::execute(v1+1,v2+1);}
然后编译器继续实例化
int DotProduct<2,int>::execute(const int v1[],const int v2[]) {return v1[0]*v2[0] + DotProduct<1,int>::execute(v1+1,v2+1);}
这里，我们有一个 DotProduct<1,T> 的偏特化版本，计算函数将被实例化为
int DotProduct<1,int>::execute(const int v1[],const int v2[]) {return v1[0]*v2[0];}
而这3个函数都足够短小，编译器通常会为它们做inline工作，使得r2的计算被展开为
r2 = v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];

DotProduct<3,int>::execute(v1,v2) 的调用语法不是很友好。
书中还设计了一个helper函数。
templateT dot_product(const T v1[],const T v2[]) {return DotProduct::execute(v1,v2);}
如同STL的make_pair、bind2nd，这个模板函数将模板参数T的推导工作交给编译器，使DotProduct的界面更友好，更不容易出错。当然，维度DIM必须指定。
现在可以按如下方式使用 ：
int r2 = dot_product<3>(v1,v2); // 是不是和 int r1 = dot_product(3,v1,v2); 十分相似？
本文以下部分将从编译器生成的机器码层次来验证这一技术，同时与普通循环点积计算做比较。
测试环境 ：Windows Xp sp2，Microsoft Visual Studio 2005 Team Suite，release默认设置
验证1：
该技术真的能和预想的一样，实例化多个函数，并将它们inline，达到解循环的效果吗？
让我们写个短小程序验证一下：
int main() {int v1[] = {1,2,3};int v2[] = {4,5,6};int r1 = dot_product(3,v1,v2);//循环版本int r2 = dot_product<3>(v1,v2);//解循环版本cout<
程序输出当然是2个32，让我们进入反汇编看看实际工作情况：
令人吃惊的是，前4行代码根本没有生成任何机器码
生成机器码的语句只有2条输出语句和return 0;
cout< >::operator<< (40203Ch)]movecx,eaxcalldword ptr [__imp_std::basic_ostream >::operator<< (402038h)]
这些指令的的含义是什么呢？

我们知道 cout<(cout<cout.opeator<<(...) 返回的是自身引用，所以可以进行链式输出
前4条指令执行后：
堆栈应该是：
std::endl
20h
ECX 保存的是 std::cout
std::cout 的类型是std::basic_ostream > 
所以第5条指令正是对cout的成员函数“ opeator<<”进行调用
(this指针 std::cout,已经存入ECX，而且可以猜想这个成员函数正是 operator(int) 重载)
而参数，对的，是call指令前，堆栈栈顶“20h”—— 32的16进制。
第5条指令将执行 cout.operator<<(20h); ！
让我们把 r1被编译器偷偷替换成20h的事放在一边，继续解析后2条指令：
cout.operator<< 使用__thiscall调用约定，函数自己弹出堆栈中参数20h， 返回自身引用，返回值存入EAX。
所以第5条指令执行后：
堆栈应该是：
std::endl;
EAX 保存的是返回值，即是std::cout
第6条指令 mov ecx,eax 将返回值重新放入ECX，设置好this指针。
第7条指令 继续调用 operator<< 的另一个重载版本。（这时 std::endl 在栈顶）
完成语句 cout.operator<<(endl);
cout<
这个原本就短小的程序，被编译器“篡改”为更短小的程序：
int main() {cout<<32<如果改用 printf 输出：
printf(“%d\n”,r1); printf(“%d\n”,r2);
从反汇编可以更清楚的看出，编译器实际做出的代码是：
printf(“%d\n”,0x20); printf(“%d\n”,0x20);
（详细代码见 sample1）
比较合理的解释是，编译器知道了太多的上下文，做出了足够的优化，将本来应该运行时得出的结果—— 一个不变的结果 ——放入的可执行文件中，直接输出。
为了验证模板元解循环技术，我们要另外想办法。
验证2：
我们需要“阻挠”编译器，让它不能猜测出运行结果，从而生成真正的计算代码。
这里使用了一个比较简单的方案：
templatevoid __stdcall MyGenerate(OutIt first,OutIt last,const char *name) {cout<<"generating "< oit(cout," ");copy(first,last,oit);cout<
用clock 函数返回值填充一个序列，即使编译器“胆大包天”，也不敢猜测clock返回结果了～
多余的 name 参数和一些输出语句是因为在release模式下，watch窗口中看不到值。
而__stdcall 是让 MyGenerate后不生成多余的平衡堆栈代码。
按照如下方式使用 ：
 
MyGenerate(v1,v1+3,”v1”);MyGenerate(v2,v2+3,”v2”);int r1 = dot_product(3,v1,v2);MyGenerate(v1,v1+3,”v1”);MyGenerate(v2,v2+3,”v2”);int r2 = dot_product<3>(v1,v2);cout<
仍然不够， 编译器会将r2的计算，推迟到 cout<
所以我们再加入一个函数
void ForceCalcResult(int result) { }
强制编译器立刻计算点积。
同时，使用函数指针调用，避免编译器对该函数inline。
程序主干如下 ：
int main() {const int dim = 3;int v1[dim];int v2[dim];void ( *const ForceCalcResult_non_inline)(int) = ForceCalcResult;MyGenerate(v1,v1+dim,"v1");MyGenerate(v2,v2+dim,"v2");int r1 = dot_product(dim,v1,v2);ForceCalcResult_non_inline(r1);MyGenerate(v1,v1+dim,"v1");MyGenerate(v2,v2+dim,"v2");int r2 = dot_product(v1,v2);ForceCalcResult_non_inline(r2);cout<这样，编译器就屈服了，乖乖的把r1和r2的计算代码展现出来。
//MyGenerate(v1,v1+dim,"v1");moveax,offset string "v1" (402134h) leaedi,[esp+24h] leaebx,[esp+18h] callMyGenerate (401220h) //MyGenerate(v2,v2+dim,"v2");moveax,offset string "v2" (402138h) leaedi,[esp+18h] leaebx,[esp+0Ch] callMyGenerate (401220h)
从这4条指令，我们可以看出v1和v2的地址：
v1[0]:esp+18h, v1[1]:esp+1Ch, v1[2]:esp+20h, v1[dim]:esp+24hv2[0]:esp+0Ch, v2[1]:esp+10h, v2[2]:esp+14h, v2[dim]:esp+18h
让我们先看模板解循环版本：
//int r2 = dot_product(v1,v2);movedi,dword ptr [esp+10h] movedx,dword ptr [esp+14h] imuledi,dword ptr [esp+1Ch] imuledx,dword ptr [esp+20h] // edi = v2[1]*v1[1]// edx = v2[2]*v1[2]moveax,dword ptr [esp+0Ch] imuleax,dword ptr [esp+18h] // eax = v2[0]*v1[0]addedi,edx addedi,eax // edi = edi+edx+eax = v2[1]*v1[1]+v2[2]*v1[2]+v2[0]*v1[0]
循环被解开了！利用3个寄存器edi,edx,eax，最终结果应该是保存在edi中。
再看接下来的代码
//ForceCalcResult_non_inline(r2);pushedicallForceCalcResult (401000h)
结果确实是存放在edi中的。 
我们接着看普通“循环”版本：
//int r1 = dot_product(dim,v1,v2);movesi,dword ptr [esp+10h] moveax,dword ptr [esp+14h] imulesi,dword ptr [esp+1Ch] imuleax,dword ptr [esp+20h] // esi = v2[1]*v1[1]// eax = v2[2]*v1[2]movecx,dword ptr [esp+0Ch] imulecx,dword ptr [esp+18h] // ecx = v2[0]*v1[0]addesi,eax addesi,ecx // esi = esi+eax+ecx = v2[1]*v1[1] + v2[2]*v1[2] + v2[0]*v1[0]//ForceCalcResult_non_inline(r1);pushesi  callForceCalcResult (401000h)
几乎相同的代码，使用的是esi,ecx,eax，结果保存在esi中。
循环同样被解开了！
编译器在我们的重重算计下，被迫在运行时计算结果，同时将运算方法展现出来；但同时，它依然不屈不饶的向我们展示它的强大优化能力！
（详细代码见sample2）
验证2+：
在验证2中，编译器知道了 const int dim = 3; 这一上下文。从而将普通循环版本的点积函数进行展开。
让我们来看看，dim取更大值时，编译器会如何。
dim=10：
普通“循环”点积计算被展开，使用esi,eax,ecx,edx，10条mov与imul指令，9条add指令，最终将计算结果存入esi
模板元点击计算也被展开，使用edi,eax,ecx,edx，10条mov与imul指令，9条add指令，最终将计算结果存入edi
dim=11：
循环点积计算发生了变化，代码如下：
//MyGenerate(v1,v1+dim,"v1");00401017  moveax,offset string "v1" (402134h) 0040101C  leaedi,[esp+68h] 00401020  leaebx,[esp+3Ch] 00401024  callMyGenerate (401280h) //MyGenerate(v2,v2+dim,"v2");00401029  moveax,offset string "v2" (402138h) 0040102E  leaedi,[esp+3Ch] 00401032  leaebx,[esp+10h] 00401036  callMyGenerate (401280h) // v1[0]:esp+3Ch, v2[0]:esp+10h//int r1 = dot_product(dim,v1,v2);0040103B  xorebp,ebp // int result = ebp=0;0040103D  xoreax,eax // eax=0;0040103F  nop// align// i=eax;// loops:00401040  movecx,dword ptr [esp+eax+10h]// ecx = *(v2+i);00401044  imulecx,dword ptr [esp+eax+3Ch]// ecx *= *(v1+i);00401049  addeax,4// ++i;0040104C  addebp,ecx // result+=ecx0040104E  cmpeax,2Ch // if (i<11)00401051  jlmain+30h (401040h) // goto loops;
编译器已经不能“忍受”这样长度的循环展开，将其老老实实的翻译成真正的循环。
但是，编译器仍然对函数做了inline工作。
对于模板元解循环，因为代码要求11层函数调用，编译器能做的只有将11层调用inline，得到的是就地展开的计算。
dim=12：
循环点积计算： 编译器甚至连inline也不做了，进行正式的函数调用。
模板元点击计算：依然就地展开。
dim=32：
循环点积计算：不展开，调用
模板元点积计算：
比较有趣的是，编译器也没有将函数就地展开，而是分成3步。就地计算前9个乘法然后与DotProduct<23,int>的返回值相加。DP<23>也没有就地计算全部，而是计算前11个，然后与DP<12>的返回值相加。
3步计算当中都没有循环。
值得注意的是，循环点积计算函数也没有很老实的按照源代码的方式进行计算，而是进行了不完全的循环展开，类似于如下代码：
template< typename T>T dot_product(int DIM,const T v1[],const T v2[]) {const int step = complierKnow;T results[step] = {0};int i=0;for (;i+step-1
DIM和step似乎没有什么固定的规律。
（详细代码见sample2，取不同的dim值即可。）
验证3：
sample2中，编译器对被计算的向量的上下文依然知情：v1,v2在main函数的栈中。
sample3做了小小的改动：dot_product_loop和dot_product_unloop将无从得知传递给它们的向量在何处。（当然，在main函数中，仍然使用函数指针来调用这2个函数，使得编译器不做inline）
sample3得到的结果和sample2基本相同，只是对向量寻址由esp+v1_offset,esp+v2_offset变为[esp+4],[esp+8]
（详细代码见sampl3）
验证4：
在sample3的基础上，通过控制台输入读取dim，使得编译对其不知情。
当然，在这种情况下，是无法使用模板元解循环的，因为它要求dim是编译时常量。
在这种情况下，循环点积计算终于被编译器翻译成真正的“循环”了。
总结：
循环版本有更大的灵活性：在dim为编译时常量时，编译器会根据其大小，进行完全解循环或部分解循环。同时也支持dim为运行时的值。
模板元版本必须使用编译时常量作为dim，而且总是完全解开循环。
循环版本与模板元版本最终都会使用相同次数的乘法与运算，区别在于乘法指令和跳转指令的数目。
模板元版本在dim较大的时候，可能会使用跳转指令，将部分工作交给更小维度的点积计算。但是乘法指令数目与维数一样。
循环版本可能会使用更多的跳转指令，使得乘法指令数目大大减少。
多出的跳转指令会占用运行时间，但也能减少目标代码体积。权衡应该使用那一种版本与权衡某个函数是否应该inline十分相似。
书中17章最后部分也提醒读者，不计后果的解循环并不总是能优化运行时间。
借用大师的观点来结束本文 —— 优化的第一原则是：不要优化。优化的第二原则（仅适用于专家）是：还是不要优化。再三测试，而后优化。
《C++编程规范》 第8条