網友在知乎的一個提問帖：
有沒有一段程式碼，讓你覺得人類的智慧也可以璀璨無比？
不一定要是完整演算法，就是那種看著看著就覺得嗨爆了，驚為天人的結構或語句。
下面是【燒茄子】引用了知名博主 Matrix67 的一篇博文：
《用三段 140 字元以內的程式碼生成一張 1024×1024 的圖片》
Kyle McCormick 在 StackExchange 上發起了一個叫做 Tweetable Mathematical Art 的比賽，參賽者需要用三條推這麼長的程式碼來生成一張圖片。
具體地說，參賽者需要用 C++ 語言編寫 RD 、 GR 、 BL 三個函式，每個函式都不能超過 140 個字元。每個函式都會接到 i 和 j 兩個整型引數（0 ≤ i, j ≤ 1023），然後需要返回一個 0 到 255 之間的整數，表示位於 (i, j) 的畫素點的顏色值。舉個例子，如果 RD(0, 0) 和 GR(0, 0) 返回的都是 0 ，但 BL(0, 0) 返回的是 255 ，那麼影象的最左上角那個畫素就是藍色。
參賽者編寫的程式碼會被插進下面這段程式當中（我做了一些細微的改動），最終會生成一個大小為 1024×1024 的圖片。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
// NOTE: compile with g++ filename.cpp -std=c++11 #include <iostream> #include <cmath> #include <cstdlib> #define DIM 1024 #define DM1 (DIM-1) #define _sq(x) ((x)*(x)) // square #define _cb(x) abs((x)*(x)*(x)) // absolute value of cube #define _cr(x) (unsigned char)(pow((x),1.0/3.0)) // cube root unsigned char GR(int,int); unsigned char BL(int,int); unsigned char RD(int i,int j){ // YOUR CODE HERE } unsigned char GR(int i,int j){ // YOUR CODE HERE } unsigned char BL(int i,int j){ // YOUR CODE HERE } void pixel_write(int,int); FILE *fp; int main(){ fp = fopen("MathPic.ppm","wb"); fprintf(fp, "P6\n%d %d\n255\n", DIM, DIM); for(int j=0;j<DIM;j++) for(int i=0;i<DIM;i++) pixel_write(i,j); fclose(fp); return 0; } void pixel_write(int i, int j){ static unsigned char color[3]; color[0] = RD(i,j)&255; color[1] = GR(i,j)&255; color[2] = BL(i,j)&255; fwrite(color, 1, 3, fp); } |
我選了一些自己比較喜歡的作品，放在下面和大家分享。
首先是一個來自 Martin Büttner 的作品：
它的程式碼如下：
|
1 2 3 4 5 6 7 8 9 10 11 |
unsigned char RD(int i,int j){ return (char)(_sq(cos(atan2(j-512,i-512)/2))*255); } unsigned char GR(int i,int j){ return (char)(_sq(cos(atan2(j-512,i-512)/2-2*acos(-1)/3))*255); } unsigned char BL(int i,int j){ return (char)(_sq(cos(atan2(j-512,i-512)/2+2*acos(-1)/3))*255); } |
同樣是來自 Martin Büttner 的作品：
這是目前暫時排名第一的作品。它的程式碼如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 |
unsigned char RD(int i,int j){ #define r(n)(rand()%n) static char c[1024][1024];return!c[i][j]?c[i][j]=!r(999)?r(256):RD((i+r(2))%1024,(j+r(2))%1024):c[i][j]; } unsigned char GR(int i,int j){ static char c[1024][1024];return!c[i][j]?c[i][j]=!r(999)?r(256):GR((i+r(2))%1024,(j+r(2))%1024):c[i][j]; } unsigned char BL(int i,int j){ static char c[1024][1024];return!c[i][j]?c[i][j]=!r(999)?r(256):BL((i+r(2))%1024,(j+r(2))%1024):c[i][j]; } |
下面這張圖片仍然出自 Martin Büttner 之手：
難以想象， Mandelbrot 分形圖形居然可以只用這麼一點程式碼畫出：
|
1 2 3 4 5 6 7 8 9 10 11 |
unsigned char RD(int i,int j){ float x=0,y=0;int k;for(k=0;k++<256;){float a=x*x-y*y+(i-768.0)/512;y=2*x*y+(j-512.0)/512;x=a;if(x*x+y*y>4)break;}return log(k)*47; } unsigned char GR(int i,int j){ float x=0,y=0;int k;for(k=0;k++<256;){float a=x*x-y*y+(i-768.0)/512;y=2*x*y+(j-512.0)/512;x=a;if(x*x+y*y>4)break;}return log(k)*47; } unsigned char BL(int i,int j){ float x=0,y=0;int k;for(k=0;k++<256;){float a=x*x-y*y+(i-768.0)/512;y=2*x*y+(j-512.0)/512;x=a;if(x*x+y*y>4)break;}return 128-log(k)*23; } |
Manuel Kasten 也製作了一個 Mandelbrot 集的圖片，與剛才不同的是，該圖描繪的是 Mandelbrot 集在某處區域性放大後的結果：
它的程式碼如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
unsigned char RD(int i,int j){ double a=0,b=0,c,d,n=0; while((c=a*a)+(d=b*b)<4&&n++<880) {b=2*a*b+j*8e-9-.645411;a=c-d+i*8e-9+.356888;} return 255*pow((n-80)/800,3.); } unsigned char GR(int i,int j){ double a=0,b=0,c,d,n=0; while((c=a*a)+(d=b*b)<4&&n++<880) {b=2*a*b+j*8e-9-.645411;a=c-d+i*8e-9+.356888;} return 255*pow((n-80)/800,.7); } unsigned char BL(int i,int j){ double a=0,b=0,c,d,n=0; while((c=a*a)+(d=b*b)<4&&n++<880) {b=2*a*b+j*8e-9-.645411;a=c-d+i*8e-9+.356888;} return 255*pow((n-80)/800,.5); } |
這是 Manuel Kasten 的另一作品：
生成這張圖片的程式碼很有意思：函式依靠 static 變數來控制繪畫的程式，完全沒有用到 i 和 j 這兩個引數！
|
1 2 3 4 5 6 7 8 9 10 11 |
unsigned char RD(int i,int j){ static double k;k+=rand()/1./RAND_MAX;int l=k;l%=512;return l>255?511-l:l; } unsigned char GR(int i,int j){ static double k;k+=rand()/1./RAND_MAX;int l=k;l%=512;return l>255?511-l:l; } unsigned char BL(int i,int j){ static double k;k+=rand()/1./RAND_MAX;int l=k;l%=512;return l>255?511-l:l; } |
這是來自 githubphagocyte 的作品：
它的程式碼如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
unsigned char RD(int i,int j){ float s=3./(j+99); float y=(j+sin((i*i+_sq(j-700)*5)/100./DIM)*35)*s; return (int((i+DIM)*s+y)%2+int((DIM*2-i)*s+y)%2)*127; } unsigned char GR(int i,int j){ float s=3./(j+99); float y=(j+sin((i*i+_sq(j-700)*5)/100./DIM)*35)*s; return (int(5*((i+DIM)*s+y))%2+int(5*((DIM*2-i)*s+y))%2)*127; } unsigned char BL(int i,int j){ float s=3./(j+99); float y=(j+sin((i*i+_sq(j-700)*5)/100./DIM)*35)*s; return (int(29*((i+DIM)*s+y))%2+int(29*((DIM*2-i)*s+y))%2)*127; } |
這是來自 githubphagocyte 的另一個作品：
這是一張使用 diffusion-limited aggregation 模型得到的圖片，程式執行起來要耗費不少時間。程式碼很有意思：巧妙地利用巨集定義，打破了函式與函式之間的界限，三段程式碼的字數限制便能合在一起使用了。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
unsigned char RD(int i,int j){ #define D DIM #define M m[(x+D+(d==0)-(d==2))%D][(y+D+(d==1)-(d==3))%D] #define R rand()%D #define B m[x][y] return(i+j)?256-(BL(i,j))/2:0; } unsigned char GR(int i,int j){ #define A static int m[D][D],e,x,y,d,c[4],f,n;if(i+j<1){for(d=D*D;d;d--){m[d%D][d/D]=d%6?0:rand()%2000?1:255;}for(n=1 return RD(i,j); } unsigned char BL(int i,int j){ A;n;n++){x=R;y=R;if(B==1){f=1;for(d=0;d<4;d++){c[d]=M;f=f<c[d]?c[d]:f;}if(f>2){B=f-1;}else{++e%=4;d=e;if(!c[e]){B=0;M=1;}}}}}return m[i][j]; } |
最後這張圖來自 Eric Tressler ：
這是由 logistic 對映得到的 Feigenbaum 分岔圖。和剛才一樣，對應的程式碼也巧妙地利用了巨集定義來節省字元：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
unsigned char RD(int i,int j){ #define A float a=0,b,k,r,x #define B int e,o #define C(x) x>255?255:x #define R return #define D DIM R BL(i,j)*(D-i)/D; } unsigned char GR(int i,int j){ #define E DM1 #define F static float #define G for( #define H r=a*1.6/D+2.4;x=1.0001*b/D R BL(i,j)*(D-j/2)/D; } unsigned char BL(int i,int j){ F c[D][D];if(i+j<1){A;B;G;a<D;a+=0.1){G b=0;b<D;b++){H;G k=0;k<D;k++){x=r*x*(1-x);if(k>D/2){e=a;o=(E*x);c[e][o]+=0.01;}}}}}R C(c[j][i])*i/D; } |
下面是 高城 的分享，伯樂線上已徵得許可
我在這裡提供我見識到的三個精彩演算法的解析，強烈地推薦給初學的演算法愛好者，它們可能會令你眼界大開，同時堅定你在演算法大道上勇往直前的信念。
#3. 二進位制是人類的好朋友，線上的樹的最近公共祖先（LCA）演算法：
利用數的二進位制表示可以產生很多加速演算法，online-LCA是其中之一。許多演算法的加速是對數率的，就是利用了數的二進位制表示。
首先定義二維陣列：prede[N+1][B+1], N表示樹的結點的數量，結點以數字1到N代指，B滿足條件：2^(B)>=N
令fa[i]表示結點i的父結點，那麼prede[i][b]的含義是：
|
1 2 |
prede[i][0] = fa[i]; prede[i][b] = prede[prede[i][b-1]][b-1]; // b >= 1 |
也就是說，prede[i][b]指的是從結點i往上走2^(b)步，所到達的結點。如果走到了盡頭，就令prede[i][b]為0。
我們只需要O(NlogN)的複雜度，就可以完成prede的初始化。此外，我們還需要預處理出所有結點的高度，也就是depth[i]，定義為：
|
1 2 |
depth[root] = 0; depth[i] = depth[fa[i]] + 1; |
當遇到詢問LCA(x, y)，我們只需要採取如下行動，就可以O(logN)的代價獲得答案：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
int lca(int x, int y) { if (depth[x] > depth[y]) swap(x, y); for(int i = B; i >= 0; i --){ //令x和y高度一致 if (depth[prede[y][i]] >= depth[x]) y = prede[y][i]; } //注意此時有可能出現x == y，那麼LCA(x,y) == x，下方的for //就不起作用了。 for(int i = B; i >= 0; i --){ //如果prede[x][i]和prede[y][i]不相同，說明這兩者的高度 //都大於所求的LCA(x,y)，也就是在LCA(x,y)的下方，此時令 //x和y一同往根部以2^(i)的步數爬升 if (prede[x][i] != prede[y][i]) x = prede[x][i], y = prede[y][i]; } if (x == y) return x; //此時LCA（x,y) = x return prede[x][0]; //此時x和y有共同的父結點 } |
上述程式碼的精髓在於兩個for(int i = B; i >= 0; i –)，這裡利用了數的二進位制表示。可以證明，對於任何嚴格小於2^(B+1)的非負整數t，下面的程式碼執行之後可以令a == t，
|
1 2 3 4 |
int a = 0; for(int i = B; i >= 0; i --){ if(a + (1<<i) <= t) a += (1<<i); } |
#2. 集合之交，樹狀陣列，動態更改、查詢陣列字首和演算法。
實現樹狀陣列所需的程式碼極為簡易，實際上它是一棵殘缺的線段樹，它可以實現一部分線段樹的功能（但凡可以化為區間求和的問題基本上都能解決），但是畢竟不如線段樹功能完整，有興趣的讀者應該學習一下線段樹的知識。
問題描述：利用預處理的字首和陣列pre[N + 1]，我們可以O(1)的代價對靜態的陣列A[N + 1]求取區間和：
pre[i] = A[0] + A[1] + A[2] + … + A[i];
A[a] + A[a+1] + A[a+2] + … + A[b] = pre[b] – pre[a-1];
但是當需要對陣列A進行動態的更改時，上述程式碼就失效了。我們需要一種演算法，可以動態地更改以及查詢字首和陣列pre[N+1]。下面首先展示樹狀陣列的程式碼，然後解釋其數學原理，它的插入和查詢的代價都是O(logN)：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 |
int Count[BiggestN+1], N; //使用前令Count所有元素為0，規定A[0]沒有數 //據，也就是說資料從A[1]開始存，pre[0]總為零 //實現功能A[i] += add void insert(int i, int add) { while( i <= N ) { Count[i] += add; i += i&(-i); } } //返回pre[i]的值 int query(int i) { int num = 0; while( i > 0 ) { num += Count[i]; i -= i&(-i); } return num; } |
演算法中最關鍵的語句是位操作i&(-i)，讀者在稿紙上算一算就可以知道：
i -= i&(-i)的功能是令i的最低的非0位變為0；
i += i&(-i)的功能是令i的最低的非0位變為0，並往更高一位進一。
理解樹狀陣列的行為，需要構造兩個集合：
Define lowbit(i) = i&(-i);
up(a) = {a, a1, a2, …}, ai = a(i-1) + lowbit(a(i-1));
down(a) = {a, a1, a2, …, 0}, ai = a(i-1) – lowbit(a(i-1));
可以證明，對於任何a <= b的正整數對(a,b)，up(a)和down(b)的交集都有且僅有一個元素。對這個定理進行含糊的說明是很容易的，a == b的情況不必考慮，a < b時，總有一個最大的i，使得b的第i位大於a的第i位（也就是b的第i位是1，而a的第i位是0）,那麼對b產生down(b)，對a產生up(a)，它們的唯一交集就是(1<<i)。注意這裡討論的第i位的i是從0開始索引的。讀者可以在稿紙上找若干數對進行實驗來加深印象。
有了上述定理，我們就不難意會insert函式和query函式的作用了。
#1. 機器浮點數的祕密，”巧奪天工”的完美例項，基於標準浮點數的快速開平方倒數演算法
這是一個公開的祕密，這是一個所有程式設計師得以欣賞的智慧之美。她在許多程式設計師的心目中高居“最美程式碼”的第一位，所有溢美之詞都無法表達他們所感受到的震撼。
一定會有許多人想在這裡貼這段程式碼，少年，來我這裡，我幫你揭開她神祕的面紗。
公式太多，貼圖講解。
