利用模板元編程實現解循環優化

簡介

在《C++ Templates: The Complete Guide》一書中（以下簡稱書），提出了模板元編程最早的實際應用之一：在數值運算中進行解循環優化。

而本文的標題是噱頭！本文的真正目的是指出這種優化措施在增加復雜性的同時，並不一定能明顯改善效率。應當謹慎使用該技術——默認不使用該技術，在點積計算確實是效率瓶頸時考慮采用該技術，並認真測試該技術是否真能提高效率。

背景

數值運算庫中常常需要提供向量點積（dot_product）運算。其定義用C++代碼描述也許更清楚～

template<typename T>
T dot_product(int dim,const T v1[],const T v2[]) {
　　T result=0;
　　for (int i=0;i<dim;++i)
　　　　result += v1[i]*v2[i];
　　return result;
}我們可以使用這個函數，求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，... ？

無疑這是一個冗長乏味易出錯的方案。

書中提出了一種使用模板元編程解決該問題的方案，讓我們來看看他是怎麼做的。

模板元編程——解循環

// 首先是遞歸模板
template<int DIM,typename T>
struct DotProduct {
　　static T execute(const T v1[],const T v2[]);
}
template<int DIM,typename T>
T DotProduct<DIM,T>::execute(const T v1[],const T v2[]) {
　　return v1[0]*v2[0] + DotProduct<DIM-1,T>::execute(v1+1,v2+1);
}
// 遞歸邊界模板
template<typename T>
struct DotProduct<1,T> {
　　static T execute(const T v1[],const T v2[]);
}
template<typename T>
T 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函數。template<int DIM,typename T>
T dot_product(const T v1[],const T v2[]) {
　　return DotProduct<DIM,T>::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<<r1<<endl;
　　cout<<r2<<endl;
　　return 0;
}

程序輸出當然是2個32，讓我們進入反匯編看看實際工作情況：

令人吃驚的是，前4行代碼根本沒有生成任何機器碼

生成機器碼的語句只有2條輸出語句和return 0;cout<<r1<<endl; // 生成7條指令
mov　　　　eax,dword ptr [__imp_std::endl (402040h)
mov　　　　ecx,dword ptr [__imp_std::cout (402044h)]
push　　　　eax
push　　　　20h
call　　dword ptr [__imp_std::basic_ostream<char,std::char_traits<char> >::operator<< (40203Ch)]
mov　　ecx,eax
call　　dword ptr [__imp_std::basic_ostream<char,std::char_traits<char> >::operator<< (402038h)]這些指令的的含義是什麼呢？

我們知道 cout<<r1<<endl; 的完整形式是

(cout<<r1)<<endl; 即是 cout.operator<<(r1).operator<<(endl);

cout.opeator<<(...) 返回的是自身引用，所以可以進行鏈式輸出

前4條指令執行後：

堆棧應該是：

std::endl

20h

ECX 保存的是 std::cout

std::cout 的類型是std::basic_ostream<char,std::char_traits<char> >

所以第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<<r2<<endl; 生成了相似的7條指令。

這個原本就短小的程序，被編譯器“篡改”為更短小的程序：

int main() {
　　cout<<32<<endl;cout<<32<<endl;
　　return 0;
}

如果改用 printf 輸出：

printf(“%d\n”,r1); printf(“%d\n”,r2);

從反匯編可以更清楚的看出，編譯器實際做出的代碼是：

printf(“%d\n”,0x20); printf(“%d\n”,0x20);

（詳細代碼見 sample1）

比較合理的解釋是，編譯器知道了太多的上下文，做出了足夠的優化，將本來應該運行時得出的結果—— 一個不變的結果 ——放入的可執行文件中，直接輸出。

為了驗證模板元解循環技術，我們要另外想辦法。

驗證2：

我們需要“阻撓”編譯器，讓它不能猜測出運行結果，從而生成真正的計算代碼。

這裡使用了一個比較簡單的方案：

template<class OutIt>
void __stdcall MyGenerate(OutIt first,OutIt last,const char *name) {
　　cout<<"generating "<<name<<" with clock\n";
　　generate(first,last,clock);
　　ostream_iterator<int> oit(cout," ");
　　copy(first,last,oit);
　　cout<<endl;
}

用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<<r1<<endl;
cout<<r2<<endl;

仍然不夠， 編譯器會將r2的計算，推遲到 cout<<r2<<endl;中，使得r2的計算不夠清晰。

所以我們再加入一個函數

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<dim>(v1,v2);
　　ForceCalcResult_non_inline(r2);
　　cout<<r1<<endl;
　　cout<<r2<<endl;
　　return 0;
}

這樣，編譯器就屈服了，乖乖的把r1和r2的計算代碼展現出來。

//MyGenerate(v1,v1+dim,"v1");
mov　　　　eax,offset string "v1" (402134h)
lea　　　　edi,[esp+24h]
lea　　　　ebx,[esp+18h]
call　　　　MyGenerate<int *> (401220h)
//MyGenerate(v2,v2+dim,"v2");
mov　　　　eax,offset string "v2" (402138h)
lea　　　　edi,[esp+18h]
lea　　　　ebx,[esp+0Ch]
call　　　　MyGenerate<int *> (401220h)從這4條指令，我們可以看出v1和v2的地址：v1[0]:esp+18h, v1[1]:esp+1Ch, v1[2]:esp+20h, v1[dim]:esp+24h
v2[0]:esp+0Ch, v2[1]:esp+10h, v2[2]:esp+14h, v2[dim]:esp+18h讓我們先看模板解循環版本：//int r2 = dot_product<dim>(v1,v2);
mov　　　　edi,dword ptr [esp+10h]
mov　　　　edx,dword ptr [esp+14h]
imul　　　　edi,dword ptr [esp+1Ch]
imul　　　　edx,dword ptr [esp+20h]
// edi = v2[1]*v1[1]
// edx = v2[2]*v1[2]
mov　　　　eax,dword ptr [esp+0Ch]
imul　　　　eax,dword ptr [esp+18h]
// eax = v2[0]*v1[0]
add　　　　edi,edx
add　　　　edi,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);
push　　　　edi
call　　　　ForceCalcResult (401000h)結果確實是存放在edi中的。

我們接著看普通“循環”版本：//int r1 = dot_product(dim,v1,v2);
mov　　　　esi,dword ptr [esp+10h]
mov　　　　eax,dword ptr [esp+14h]
imul　　　　esi,dword ptr [esp+1Ch]
imul　　　　eax,dword ptr [esp+20h]
// esi = v2[1]*v1[1]
// eax = v2[2]*v1[2]
mov　　　　ecx,dword ptr [esp+0Ch]
imul　　　　ecx,dword ptr [esp+18h]
// ecx = v2[0]*v1[0]
add　　　　esi,eax
add　　　　esi,ecx
// esi = esi+eax+ecx = v2[1]*v1[1] + v2[2]*v1[2] + v2[0]*v1[0]
//ForceCalcResult_non_inline(r1);
push　　　　esi　
call　　　　ForceCalcResult (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　mov　　eax,offset string "v1" (402134h)
0040101C　lea　　edi,[esp+68h]
00401020　lea　　ebx,[esp+3Ch]
00401024　call　　MyGenerate<int *> (401280h)
//MyGenerate(v2,v2+dim,"v2");
00401029　mov　　eax,offset string "v2" (402138h)
0040102E　lea　　edi,[esp+3Ch]
00401032　lea　　ebx,[esp+10h]
00401036　call　　MyGenerate<int *> (401280h)
// v1[0]:esp+3Ch, v2[0]:esp+10h
//int r1 = dot_product(dim,v1,v2);
0040103B　xor　　ebp,ebp
// int result = ebp=0;
0040103D　xor　　eax,eax
// eax=0;
0040103F　nop
// align
// i=eax;
// loops:
00401040　mov　　ecx,dword ptr [esp+eax+10h]
// ecx = *(v2+i);
00401044　imul　　ecx,dword ptr [esp+eax+3Ch]
// ecx *= *(v1+i);
00401049　add　　eax,4
// ++i;
0040104C　add　　ebp,ecx
// result+=ecx
0040104E　cmp　　eax,2Ch
// if (i<11)
00401051　jl　　　　main+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;i+=step) {
　　　　results[0] += v1[i]*v2[i];
　　　　results[1] += v1[i+1]*v2[i+1];
　　　　...
　　　　results[step-1] += v1[i+step-1]*v2[i+step-1];
　　}
　　for (;i<DIM;++i)
　　　　results[0] += v1[i]*v2[i];
　　return results[0]+result[1]+...result[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條

本文配套源碼