為何不同標準庫實現的三角函數的執行效率差別如此巨大？

01-11

現在看，似乎就是標準庫實現的不同。改成double在gcc里居然快了一倍。再也不相信標準庫了。。。。
============================================
代碼非常簡單，四個標準庫的三角函數反覆調用N次：

#include & #include & #include & #include &
float result[1024];
int main(int argc, char** argv) { FILE* fh = std::fopen("simple_tri.txt", "wb"); for (int iter = 0; iter &< 50000; iter++) { for (int i = 0; i &< 1024; i++) { float tmp = i; result[i] = std::tan(tmp) + std::tan(2*tmp) + std::sin(tmp) + std::cos(tmp); } std::fwrite(result, sizeof(float), 1024, fh); // reset result for (int i = 0; i &< 1024; i++) { result[i] = 0.0f; } // write again std::fwrite(result, sizeof(float), 1024, fh); } fclose(fh); }
兩次write到文件以防任何可能的優化，實際上去掉write對執行時間也沒有本質的影響。分別在Linux下面使用gcc 4.9.2、Windows里使用VS2013、MinGW 4.8編譯。VS2013使用release模式，gcc和MinGW開-O3、-mavx、-ffast-math。
進程執行時間如下：

gcc 4.9.2 @ Debian Jessie：25.816s

VS2013 @ Windows 10：1.3041494s

MinGW 4.8 @ Windows 10：13.0604687s

執行效率居然差出了將近20倍。這是為什麼？

====================================
GCC換了下更激進的優化參數，然並卵，依然是那麼長時間。
-march=native -O3 -ffast-math -g

VC提供了可以一次算兩個值的版本，充分利用SIMD特性，而glibc裡面只有__kernel_XXXf(sysdeps/ieee754flt-32/k_XXXf.c)，雖然編譯後也用了SSE，但一次只能算一個值，這裡差了兩倍。即便源碼稍做調整，讓VC也調用一次算一個值的版本，性能差距仍有將近10倍，所以還有其它原因。

為了進一步研究造成這種性能差距的原因，答主將glibc的代碼移植到了windows上，實現了零外部依賴的三角函數，最終使用VC編譯，這樣可以排除操作系統和硬體的干擾，也方便做profile。簡單起見，下文中僅比較兩者tanf的實現，核心代碼如下。

for (int i = 0; i &< total; i++) { float tmp = i; result[i] = tanf(tmp); }

做profile之前，答主對比了GCC和VC生成的可執行文件中的機器(彙編)代碼，發現基本相似，有可比性。

profile結果顯示，glibc版本三角函數只有10.1%的時間花在__kernel_XXXf上，所以三角函數核心演算法中SSE指令集的區別對性能的影響是有限的，這點答主原先的理解是有誤的，已經贊同過的知友請注意。

真正的瓶頸在一個名為__ieee754_rem_pio2f函數上，佔用了88.4%的時間。三角函數是周期函數，對於超過周期或半周期的值要先進行規約再計算，__ieee754_rem_pio2f實現的就是這個功能。但pi是無理數，無法用浮點數精確表示，簡單的乘除法對pi求余肯定是有誤差的。為了提高精度，減少不必要的誤差，__ieee754_rem_pio2f中要做的工作很多。

那麼問題來了，VC的實現究竟是怎樣的，以至於有10倍以上的性能差距？是否犧牲了計算精度？經過比較，前1024個數中，VC實現所得的結果與glibc實現所得的結果有33個不同，其中一共60位不同，也就是說二者精度非常接近。

通過VC編譯時生成的pdb可以發現源文件是手寫的彙編代碼，我們只能通過反彙編可執行文件來觀察其實現邏輯。這部分答主仍然在做，目前來看，應該是直接使用更高精度的double來計算簡化了很多精度控制邏輯。

鑒於glibc中這部分代碼編寫在1993年，那時intel pentium剛剛問世，SSE還不知為何物，有差距也是正常的。double版本可能最早成形於2001年，對SIMD會更友善些，相應的性能有提升也是正常的。

有更詳細結論後再更新。

閑的無聊，用release編譯一下看了看反彙編代碼，確認了幾件事情：

MSVC沒有把三角函數運算摺疊為常量
最外層的50000次循環也沒有被優化掉
三角函數的內部實現使用了CPU的SSE指令集，速度更快
MSVC Release的默認設置是會生成使用SSE2指令集的代碼的

我的cpu是core i7 4810MQ，支持SSE指令集

VS會很智能地根據目標平台的硬體做一些優化，如果題主想要用GCC達到同樣效果的話，推薦在編譯參數後面加上-march=corei7 （如果你自己的cpu不是i7，那麼需要查一下march支持的cpu代號，填上你自己的那個）

最後放反彙編代碼：

#define _CRT_SECURE_NO_WARNINGS #include & #include & #include & #include &


float result[1024];

int main(int argc, char** argv) { 00031002 DC 83 EC 08 83 E4 fadd qword ptr [ebx-1B7CF714h] 00031008 F0 83 C4 04 lock add esp,4 0003100C 55 push ebp 0003100D 8B 6B 04 mov ebp,dword ptr [ebx+4] 00031010 89 6C 24 04 mov dword ptr [esp+4],ebp 00031014 8B EC mov ebp,esp 00031016 83 EC 38 sub esp,38h 00031019 56 push esi 0003101A 57 push edi FILE* fh = std::fopen("simple_tri.txt", "wb"); 0003101B 68 08 31 03 00 push 33108h 00031020 68 0C 31 03 00 push 3310Ch 00031025 FF 15 B8 30 03 00 call dword ptr ds:[330B8h] 0003102B 83 C4 08 add esp,8 0003102E 89 45 F8 mov dword ptr [ebp-8],eax 00031031 C7 45 FC 50 C3 00 00 mov dword ptr [ebp-4],0C350h 00031038 0F 1F 84 00 00 00 00 00 nop dword ptr [eax+eax] for (int iter = 0; iter &< 50000; iter++) { for (int i = 0; i &< 1024; i++) 00031040 33 FF xor edi,edi 00031042 BE 98 43 03 00 mov esi,34398h { float tmp = i; 00031047 66 0F 6E C7 movd xmm0,edi 0003104B 66 0F 70 C0 00 pshufd xmm0,xmm0,0 00031050 66 0F FE 05 20 31 03 00 paddd xmm0,xmmword ptr ds:[33120h] 00031058 0F 5B C8 cvtdq2ps xmm1,xmm0 result[i] = std::tan(tmp) + std::tan(2 * tmp) + std::sin(tmp) + std::cos(tmp); 0003105B 0F 5A C1 cvtps2pd xmm0,xmm1 0003105E 0F 29 45 D0 movaps xmmword ptr [ebp-30h],xmm0 00031062 0F 28 05 30 31 03 00 movaps xmm0,xmmword ptr ds:[33130h] 00031069 0F 59 C1 mulps xmm0,xmm1 0003106C 0F 5A C0 cvtps2pd xmm0,xmm0 0003106F E8 BC 0C 00 00 call ___vdecl_tan2 (031D30h) 00031074 66 0F 5A C0 cvtpd2ps xmm0,xmm0 00031078 0F 29 45 E0 movaps xmmword ptr [ebp-20h],xmm0 0003107C 0F 28 45 D0 movaps xmm0,xmmword ptr [ebp-30h] 00031080 E8 AB 0C 00 00 call ___vdecl_tan2 (031D30h) 00031085 66 0F 5A C8 cvtpd2ps xmm1,xmm0 00031089 0F 28 45 E0 movaps xmm0,xmmword ptr [ebp-20h] 0003108D 0F 58 C1 addps xmm0,xmm1 00031090 0F 29 45 E0 movaps xmmword ptr [ebp-20h],xmm0 00031094 0F 28 45 D0 movaps xmm0,xmmword ptr [ebp-30h] 00031098 E8 83 0C 00 00 call ___vdecl_sin2 (031D20h) 0003109D 66 0F 5A C8 cvtpd2ps xmm1,xmm0 000310A1 0F 28 45 E0 movaps xmm0,xmmword ptr [ebp-20h] 000310A5 0F 58 C1 addps xmm0,xmm1 000310A8 0F 29 45 E0 movaps xmmword ptr [ebp-20h],xmm0 000310AC 0F 28 45 D0 movaps xmm0,xmmword ptr [ebp-30h] 000310B0 E8 5B 0C 00 00 call ___vdecl_cos2 (031D10h) 000310B5 0F 28 4D E0 movaps xmm1,xmmword ptr [ebp-20h] 000310B9 83 C7 02 add edi,2 000310BC 66 0F 5A C0 cvtpd2ps xmm0,xmm0 000310C0 0F 58 C8 addps xmm1,xmm0 000310C3 F2 0F 11 0E movsd mmword ptr [esi],xmm1 000310C7 83 C6 08 add esi,8 000310CA 81 FE 98 53 03 00 cmp esi,35398h 000310D0 0F 8C 71 FF FF FF jl main+47h (031047h) } std::fwrite(result, sizeof(float), 1024, fh); 000310D6 8B 75 F8 mov esi,dword ptr [fh] 000310D9 56 push esi 000310DA 68 00 04 00 00 push 400h 000310DF 6A 04 push 4 000310E1 68 98 43 03 00 push 34398h 000310E6 FF 15 B0 30 03 00 call dword ptr ds:[330B0h] // reset result for (int i = 0; i &< 1024; i++) { result[i] = 0.0f; } // write again std::fwrite(result, sizeof(float), 1024, fh); 000310EC 56 push esi 000310ED 68 00 04 00 00 push 400h 000310F2 33 C0 xor eax,eax 000310F4 B9 00 04 00 00 mov ecx,400h 000310F9 BF 98 43 03 00 mov edi,34398h 000310FE 6A 04 push 4 00031100 F3 AB rep stos dword ptr es:[edi] 00031102 68 98 43 03 00 push 34398h // reset result for (int i = 0; i &< 1024; i++) { result[i] = 0.0f; } // write again std::fwrite(result, sizeof(float), 1024, fh); 00031107 FF 15 B0 30 03 00 call dword ptr ds:[330B0h] 0003110D 83 C4 20 add esp,20h 00031110 83 6D FC 01 sub dword ptr [ebp-4],1 00031114 0F 85 26 FF FF FF jne main+40h (031040h) } fclose(fh); 0003111A 56 push esi 0003111B FF 15 B4 30 03 00 call dword ptr ds:[330B4h] 00031121 83 C4 04 add esp,4 } 00031124 33 C0 xor eax,eax 00031126 5F pop edi 00031127 5E pop esi 00031128 8B E5 mov esp,ebp 0003112A 5D pop ebp 0003112B 8B E3 mov esp,ebx 0003112D 5B pop ebx 0003112E C3 ret

按照 @孫明琦提供的彙編代碼，應該沒有做優化，所以不要考慮下面的答案。

============================

看一下彙編吧，我覺得可能是VC把三角函數的調用都解析為常量表達式了。所以根本沒有去算而是直接輸出了數據。

建議題主把輸入放到文件里，程序啟動時一次性讀進來，然後再開始計時計算。

可以試試手動按照級數展開的方法計算，看看速度是否會加快！

不用進程測，改為直接用clock()測呢？