数学与算法的艺术
為了更好體現數學在計算機算法中的重要作用，我來介紹一些數學知識以及它們的應用。
?
Contens
?
??? 1. 秦九韶算法???????????????? 9.? 最小二乘法
??? 2. 斯特林公式?????????????????10. 自守數
??? 3. 外觀數列
??? 4. 整數拆分問題
??? 5. 阿貝爾變換
??? 6. 二項式反演
??? 7. 馬青公式
??? 8. 艾森斯坦判別法
?
?
1. 秦九韶算法
?
???秦九韶算法是中國南宋時期數學家秦九韶提出的一種計算多項式的優化算法。在西方又被稱為Horner算法。
?
?? 給定一個多項式，計算這個多項式的值。
?
???可以對這個多項式進行如下改寫
?
???
?
??? 那么由內到外就是
???
????
???
????時間復雜度為。
?
代碼：
int Horner(int a[], int n, int x) {int v1 = a[n];for(int i = n - 1; i >= 0; i--){int v2 = v1 * x + a[i];v1 = v2;}return v1; }
?
2. 斯特林公式
?
?? 斯特林公式是用來取階乘近似值的數學公式。尤其是在比較大的時候，用斯特林公式非常高效。公式如下
?
??
?
?? 在維基百科上有關于斯特林數更詳細的介紹，如下
?
?? 鏈接：http://zh.wikipedia.org/wiki/%E6%96%AF%E7%89%B9%E9%9D%88%E5%85%AC%E5%BC%8F
?
?? 題目：http://acm.hdu.edu.cn/showproblem.php?pid=1018
?
?? 分析：題意就是求階乘的位數，由于給定數比較大，直接做很容易TLE，采用斯特林數易過之。
??
?? 拓展：給定一個正整數，比較大，求如下表達式結果的位數。
?
???????
??
???分析：作如下變換，然后取對數即可。
?
????????
?
3. 外觀數列
??
???外觀數列是這樣一個數列，首項為1，以后的每一項都是對前一項的描述，數列如下
?
??
?
?? 比如312211是對111221的描述，即3個1，2個2，1個1。
?
?? 1987年，康威（John Conway）發現，在這個數列中，相鄰兩數的長度之比越來越接近一個固定的常數，當數列
?? 項數趨近無窮大的時候，該比值約為1.303577，此常數叫做康威常數。外觀數列還有一個性質，4永遠不出現。
?
?? 更多有關康威常數的了解請關注Matrix 67的博客，如下
?
?? 鏈接：http://www.matrix67.com/blog/archives/3870
?
?? 題目：http://acm.hdu.edu.cn/showproblem.php?pid=4148
?
?? 題意：求外觀數列指定某一項的長度。
?
?? 分析：我們只需考慮前一項與后一項的關系，然后遞推即可。
?
代碼：
#include <iostream> #include <string.h> #include <stdio.h>using namespace std; char list[32][10005];void Init() {for(int i = 0; i < 32; i++)memset(list, 0, sizeof(list));list[1][0] = '1';for(int i = 2; i < 32; i++){int k = 0;int cnt = 1;int len = strlen(list[i - 1]);for(int j = 0; j < len; j++){if(list[i - 1][j] == list[i - 1][j + 1]) cnt++;else{list[i][k++] = cnt + '0';list[i][k++] = list[i - 1][j];cnt = 1;}}} }int main() {int n;Init();while(cin>>n && n)cout<<strlen(list[n])<<endl;return 0; }
?
4. 整數拆分問題
?
?? 問題：給定一個自然數，把它拆分為若干個數的和，記這若干個數的積為，求的最大值。
?
?? 分析：假設把拆分為個數之和，那么有和，由不等式
?
????????
?
??????? 進一步得到如下結論
?
???????
?
??????? 我們需要求的最大值，等價于求如下函數的最大值
?
???????
?
??????? 求導后得到
?
???????
?
??????? 那么函數在處取得最大值，因為與3最接近，所以可得到如下結論
?
??????? 結論：，使乘積最大的拆分方案是，先拆分出盡多可能的3，如果剩余
???????????? 的數為2或0，則拆分結束。如果剩余的數為1，則將拆分好的3拿一個和剩余的1組合拆分為兩個2，
?????????????拆分結束。
?
????????而本題也可以用dp的思路做，只不過時間復雜度比較高，僅適合比較小的。狀態轉移方程如下
?
??????? ?
?
??????? 其中表示把自然數拆分為若干數之和，得到的最大乘積。
?
?
5. 阿貝爾變換
?
?? 給定兩個數列和，并記
?
??
?
?? 并且，那么有如下結論
?
??
?
?? 則稱上式為阿貝爾變換或者阿貝爾部分求和公式。證明步驟如下
?
?? 由于，則進一步有
?
??
?
?? 證明完畢！
?
?? 應用阿貝爾變換可以解決一些比較復雜的數學問題，接下來看看一個經典例題。
??
?? 例題：已知，且滿足如下兩個條件
?
???????
?
???????證明如下不等式
?
??????
?
?? 證明：記中全體非負實數之和為，全體負數之和為，由條件可知且。所以得到
?
???????
?
????????記，并且
?
???????
?
??????? 很明顯，對于所有的，都有
?
???????
???????
????????所以根據阿貝爾變換得到
?
???????
???????
????????證明完畢！
?
?
6. 二項式反演
?
?? 設和是定義在非負整數集上的兩個函數，并且滿足
?
???
?
???那么得到
?
???
?
???這就是二項式反演的主要內容。
?
?? 證明：首先我們要知道
?
???
?
???有了這個，就很容易證明了，步驟如下
?
???
?
?? 證明完畢！
?
?? 接下來，我會用二項式反演來求錯排公式。
?
?? 個不同的元素排列的方法數為，那么恰有個元素錯排的方法數為，那么有
?
???
?
?? 根據二項式反演，進而得到
?
???
?
?
7. 馬青公式
??
?? 馬青公式是英國天文學教授約翰·馬青于1706年發現的，用它來計算圓周率的值非常高效。公式內容如下
?
???
??
???關于這個公式的推導參考WiKi，鏈接如下
?
?? 鏈接：http://zh.wikipedia.org/wiki/%E6%A2%85%E6%AC%BD%E9%A1%9E%E5%85%AC%E5%BC%8F
?
?? 對于的冪級數展開式推導步驟如下
?
?? 設，那么有
?
??
??
???兩邊分別積分，最終得到
?
???
?
???也就是說
?
???
?
???然后得到的計算方法如下
?
??
?
?? 通過這個公式，可以很容易計算出小數點后精確到某一位的值。當然這里計算有技巧，由于
?
??
?
?? 設，并且
?
??
?
???對于可以用秦九韶算法求之，最終得到，為了提高效率，乘法運算可用FFT。
?
代碼：
#include <iostream> #include <string.h> #include <stdio.h>using namespace std; typedef long long LL;const int BASE = 10000; const int DIGIT = 10000; const int NUM = DIGIT >> 2;int res[NUM], A[NUM], B[NUM], C[NUM];void set(int res[], int x) {for(int i = 0; i < NUM; i++)res[i] = 0;res[0] = x; }void copy(int res[], int x[]) {for(int i = 0; i < NUM; i++)res[i] = x[i]; }bool zero(int res[]) {for(int i = 0; i < NUM; i++)if(res[i]) return false;return true; }void add(int res[], int x[]) {for(int i = NUM - 1; i >= 0; i--){res[i] += x[i];if(res[i] >= BASE && i > 0){res[i] -= BASE;res[i - 1]++;}} }void sub(int res[], int x[]) {for(int i = NUM - 1; i >= 0; i--){res[i] -= x[i];if(res[i] < 0 && i > 0){res[i] += BASE;res[i - 1]--;}} }void multi(int res[], int x) {int t = 0;for(int i = NUM - 1; i >= 0; i--){res[i] *= x;res[i] += t;t = res[i] / BASE;res[i] %= BASE;} }void div(int res[], int x) {int t = 0;for(int i = 0; i < NUM; i++){res[i] += t * BASE;t = res[i] % x;res[i] /= x;} }void arctan(int res[], int A[], int B[], int x) {int m = x * x;int cnt = 1;set(res, 1);div(res, x);copy(A, res);do{div(A, m);copy(B, A);div(B, 2 * cnt + 1);if(cnt & 1)sub(res, B);elseadd(res, B);cnt++;}while(!zero(B)); }void print(int res[]) {printf("%d.\n", res[0]);for(int i = 1; i < NUM; i++){printf("%04d", res[i]);if(i % 15 == 0) puts("");}puts(""); }void Machin() {arctan(res, A, B, 5);multi(res, 4);arctan(C, A, B, 239);sub(res, C);multi(res, 4);print(res); }int main() {Machin();return 0; }
?
8. 艾森斯坦判別法
?
?? 艾森斯坦判別法是目前為止用來判斷內一個多項式可約與否的最好的方法。內容如下
?
???給定一個次本原多項式，并且這是一個整系數多項式，若存在一
?? 個素數，使得下面3個條件都成立
?
??
?
???那么為不可約多項式。
?
?? 上面提到了本原多項式，其實本原多項式就是滿足多項式所有系數的最大公約數為1的多項式。
?
?
9. 最小二乘法
?
?? 最小二乘法又叫最小平方法，是一種數學優化技術，它通過最小化誤差的平方和尋找數據的最佳函數匹配。
?
?? 通常情況，最小二乘法用于求回歸問題。以簡單的線性最小二乘為例，在二維平面上給定個點的坐標，確定
?? 一條直線，要求盡量吻合這個點的走向。
?
?? 首先設要找的直線方程為，這個點的坐標為，那么現在要
?? 確定和的值，使得
?
??
?
???最小，根據這種方法來求和的值就是典型的最小二乘法。
?
?? 可以看出是關于和的一個二元函數，要求的最小值，先求偏導，并且偏導值均為零，即
?
??
?
???進一步得到
?
???
?
?? 然后聯立兩式可以解出和，如果方程個數比較多，可以通過高斯消元來做，最小二乘法的解是唯一的。
?
?? 問題：在一個平面直角坐標系中給定個點，找出一條直線，要求這些點到這條直線距離
??????? 的平方和最小，求出這條直線的方程。
?
?? 問題：在一個三維空間中給定個點，求一過原點的平面，使得這些點到這個平面的距
??????? 離的平方和最小，求出這個平面的方程。（2014年編程之美復賽）
?
?
10. 自守數
?
?? 一個位的自然數，如果它的平方后得到的最后位跟原數相同，那么稱為自守數。其數學表達式為
?
??
?
?? 一位的自守數有3個，即1，5，6。一位數以上的自守數分為兩類
?
? （1）個位數為5，形式為
? （2）個位數為6，形式為
?
?? 可以看出兩個位數相同的自守數之和為。自守數還有兩個重要的性質
?
??（1）一個數為自守數當且僅當它為另一個自守數的后綴
??（2）位（1位數除外）的自守數僅有兩個，其位數包括前導零
?
?? 自守數的計算方法，一個位的自然數可以由求得，具體如下
?
??（1）對于個位數為5的自守數
?
????? 求的平方，取最后位，如果第位為零，則取最后位。例如
?
?????
?
??（2）對于個位數為6的自守數
?
????? 求的平方，取最后位，把最后第位的數用10減之代替，若第位是零，則取最后
??????第位來減，第位保持為零。例如
?
?????
?
?? 題目：http://acm.timus.ru/problem.aspx?space=1&num=1698
?
?? 題意：求位數不大于的自守數的個數。
?
代碼：
#include <iostream> #include <string.h> #include <stdio.h>using namespace std; typedef long long LL; const int N = 2005;int a[N], b[N]; int ans, n;bool check(int k) {b[k] = 0;for(int i = 0; i <= k; i++)b[k] += a[i] * a[k-i];if(k) b[k] += b[k-1] / 10;return b[k] % 10 == a[k]; }void dfs(int k) {if(k >= n) return;for(int i=9; i >= 0; i--){a[k] = i;if(check(k)){if(a[k]) ans++;dfs(k+1);}} }int main() {while(cin>>n){ans = 0;dfs(0);cout<<ans<<endl;}return 0; }
?
?
?
總結
