淺談最佳化程式效能(上)

前言

我們知道，多項式定義為：

在幾何學中，多項式是最簡單的平滑曲線。簡單是指它僅由乘法及加法構成，平滑是因為它類同口語中的平滑，以數學術語來說，它是無限可微，即它的所有高次微分都存在。事實上，多項式的微分也是多項式。簡單及平滑的特點，使多項式在數值分析、圖論，以及電腦繪圖等，都發揮極大的作用。多項式求值是解決許多問題的核心技術。以數值分析為例，多項式函數常常用作對數學庫中的三角函數求近似值。

現在，讓我們來用 C 語言寫一個對多項式求值的函數吧。

直接的演算法

直接按照多項式的定義使用迴圈求值：

double poly(double a[], double x){  double result = 0, p = 1;  for (int i = 0; i < N; i++, p *= x) result += a[i] * p;  return result;}

這個演算法需要執行 2N 個乘法和 N 個加法。

秦九韶演算法(Horner's method)

我們可以通過使用秦九韶演算法，或者被大多數外國人以及一些中國人稱為 Horner's method 的演算法，來減少乘法的數量：

相應的 C 語言程式如下所示：

double polyh(double a[], double x){  double result = 0;  for (int i = N - 1; i >= 0; i--) result = result * x + a[i];  return result;}

這個演算法需要執行 N 個乘法和 N 個加法。乘法的數量是原始演算法的一半。

看來，通過使用秦九韶演算法，我們大大地最佳化了程式的效能。

且慢，還是以事實說話吧，讓我們來作一些測試吧。

測試程式

下面就是比較這兩個演算法效能優劣的測試程式 poly.c ：

#include <stdio.h>#include <time.h>#define N 23456789void initialize(double a[]){  for (int i = 0; i < N; i++) a[i] = i - 12345678.9012345;}double poly(double a[], double x){  double result = 0, p = 1;  for (int i = 0; i < N; i++, p *= x) result += a[i] * p;  return result;}double polyh(double a[], double x){  double result = 0;  for (int i = N - 1; i >= 0; i--) result = result * x + a[i];  return result;}int main(int argc, char *argv[]){  static double a[N];  initialize(a);  double result = 0, (*func)(double*, double) = (argc > 1) ? polyh : poly;  clock_t elapsed = clock();  for (int i = 0; i < 1234; i++) result += func(a, 2.34 / (i - 1234567) - 1);  elapsed = clock() - elapsed;  printf("%p %g %11f\n", func, result, (double)elapsed / CLOCKS_PER_SEC);}

這個測試程式通過對一個二十多萬項的多項式用不同的自變數值求值一千多遍來比較兩個演算法的優劣。測試程式中的各項參數經過精心選擇，不會在求值過程中造成浮點溢出。在 C 語言中，每個雙精確度浮點數佔用八個位元組，所以這個測試程式大約需要一百八十MB的記憶體空間來運行。

測試結果

我們在 32-bit Windows 作業系統和 64-bit Linux 作業系統下均進行了測試，結果如下所示：

D:\work> cl /O2 poly.cpp用於 80x86 的 Microsoft (R) 32 位 C/C++ 最佳化編譯器 16.00.40219.01 版著作權(C) Microsoft Corporation。著作權所有，並保留一切權利。poly.cppMicrosoft (R) Incremental Linker Version 10.00.40219.01Copyright (C) Microsoft Corporation.  All rights reserved./out:poly.exepoly.objD:\work> poly00871000 1.4271e+029  107.531000D:\work> poly h00871080 1.4271e+029  127.686000D:\work> poly h00871080 1.4271e+029  127.686000D:\work> poly00871000 1.4271e+029  106.860000D:\work> poly00871000 1.4271e+029  107.935000D:\work> poly h00871080 1.4271e+029  128.662000

ben@vbox:~/work> gcc -std=c99 -O2 poly.cben@vbox:~/work> ./a.out0x400640 1.4271e+29  125.240000ben@vbox:~/work> ./a.out0x400640 1.4271e+29  124.820000ben@vbox:~/work> ./a.out h0x400680 1.4271e+29  160.890000ben@vbox:~/work> ./a.out h0x400680 1.4271e+29  162.210000ben@vbox:~/work> ./a.out0x400640 1.4271e+29  125.450000ben@vbox:~/work> ./a.out h0x400680 1.4271e+29  162.230000

上述運行結果中，第一欄是 func 變數的值，表示指向 poly 或 polyh 函數的指標。第二欄是對多項式求值結果的總和。第三欄是已耗用時間，單位是秒。總結一下，如下表所示：

作業系統	Windows (32-bit)		Linux (64-bit)
所用演算法	直接的	秦九韶	直接的	秦九韶
函數地址	00871000	00871080	0x400640	0x400680
計算結果	1.4271e+029	1.4271e+029	1.4271e+29	1.4271e+29
1	107.531	127.686	125.24	160.89
2	106.860	127.686	124.82	162.21
3	107.935	128.662	125.45	162.23
平均(秒)	107.442	128.011	125.17	161.78

測試說明

是不是非常意外？無論是在 32-bit Windows 作業系統，還是在 64-bit Linux 作業系統中測試，多次測試結果都表明秦九韶演算法比原始的演算法慢，同一演算法在同一作業系統中的多次啟動並執行所需的時間也很接近。由於我們使用的 Linux 作業系統是 64-bit 的，同一演算法比在 32-bit 的 Windows 作業系統運行要稍微慢一點。所有的測試結果中，無論是用哪種演算法，還是在不同的作業系統中，對多項式求值結果的總和都是相同的，這是預料之中，如果不同就有問題了。

是測試程式運行時的記憶體和 CPU 使用方式，由於是雙核的 CPU，測試程式運行時 CPU 使用率為百分之五十左右：

可以看出進程數為83個，記憶體使用量965MB。測試程式運行完畢，進程數減少一個，記憶體也降到795MB：

是在 Linux 作業系統中測試程式的健全狀態，由於 VirtualBox 只給 Linux 作業系統分配一個 CPU，所以測試程式運行時 CPU 佔用接近百分之百，記憶體佔用大約 179 MB，如所示：

測試環境

本次測試是在 Dell inspiron 1520 本本上啟動並執行，該本本只有一個 CPU，是雙核心的 Intel Core2 Duo。安裝了 Windows Vista Hoem Premium SP2 (32-bit) 作業系統。

在 Oracle VM VirtualBox 上啟動並執行 openSuSE 12.1 (64-bit) 作業系統：

結論

通過以上論述，說明了最小化一個計算中的運算元量不一定會提高它的效能。至於出現這個情況的原因，敬請期待本文的下篇。

參考資料

1. Wikipedia: Polynomial
2. Wikipedia: Horner's method
3. 維基百科: 秦九韶演算法
4. C++ Reference: clock