代碼最佳化執行個體:分數化小數

http://community.csdn.net/Expert/topic/5563/5563568.xml

題目描述：

將分數轉化為小數，相信很多人都會吧．在電腦中並能直接進行分數運算，需要將分數轉換化為浮點數或雙精確度數才能運算，但這樣會導致結果的不精確，那麼，這裡給定一個分數N/D,N為分子，D為分母(N,D均為整數)，請給出分數精確運算的方法並編程求出N/D的精確小數形式，當然如果這個小數為無限迴圈小數，則把迴圈的部分用括弧括起來，接著迴圈的部分則省略不寫。比如：
1/3    =0.(3)
22/5=4.4
1/7    =0.(142857)
2/2    =1.0
3/8    =0.375
45/56    =0.803(571428)

比如計算N/D=1/1000000007

分析：

這個我覺得沒有什麼可最佳化的，輸出結果花費的時間比例太高了.

先給以個最簡單的程式：

#include <stdio.h>
#ifdef WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif

void output_intpart(int intpart, int shift){
    int i, v;
    v=1;
    for(i=0;i<shift;i++)v*=10;
    printf("%d.",intpart/v);
    printf("%0*d",shift,intpart);
}

int main(int argc, char *argv[]){
    int N = atoi(argv[1]);
    int D = atoi(argv[2]);
    longlong LN;
    int count2=0,count5=0;
    int shift,i;
    int intpart;
    while(D%2==0){count2++;D/=2;}
    while(D%5==0){count5++;D/=5;}
    if(count2>count5){
        shift=count2;
        LN = N;
        for(i=0;i<count2-count5;i++)LN*=5;
    }else if(count2<count5){
        shift=count5;
        LN = N;
        for(i=0;i<count5-count2;i++)LN*=2;
    }else{
        shift=count2;
        LN = N;
    }
    intpart=LN/D;
    N = LN%D;
    output_intpart(intpart, shift);
    if(N>0){
        int cur=N;
        printf("(");
        do{
            int d;
            cur*=10;
            d=cur/D;
            printf("%d",d);
            cur=cur%D;
        }while(cur!=N);
        printf(")");
    }
    printf("/n");
    return 0;
}

先不考慮輸入輸出，可以看出：

代碼中主要部分在N>0後面的那個迴圈
我們知道除法效率比較低，一個可以考慮的是將除法轉化為乘法。
由於這裡D相對於這個迴圈是常數，這個是可以做到的，（參考http://blog.csdn.net/mathe/archive/2006/09/01/1153575.aspx）
不過產生的代碼還是需要用彙編代碼來寫才行，不然效率還是不行。

還有一種可行的方法是這個迴圈中本身用到的就是很多類似
  10*a=u*D+v (1<=v<D)
的這種分解過程。
  如果我們對所以的a (1<=a<D)都事先分解出來，那麼就不需要使用乘除法了（只需要加減運算就可以了）。不過問題在於對於大的D,我們需要很大的記憶體空間。對於太大的D，2G記憶體的地址空間就不夠用了。

還有一種可行方法是將迴圈中每次
   cur*=10
改成
   cur*=10^k (1<=k<=10)
這樣每次可以得到多位而不是一位結果，（當然，在每個迴圈中，我們需要檢查其每一位來判斷是否已經出現迴圈）。不過這種方法最多也就提高10倍速度。

 medie2005(阿諾) ：不過呢，我已經用數學方法求得小數的迴圈節長度，因此mathe所說的：“當然，在每個迴圈中，我們需要檢查其每一位來判斷是否已經出現迴圈”在我的方法中是不需要的。通過測試，正如mathe所說，效率確實只能提高10倍左右。至於是否還能再最佳化，我降低一下標準，將原來的10秒改成30秒，如果有人達到30秒之內，就可以了

 事先計算迴圈節長度不難，需要對整數D做因子分解。由於D的範圍不大，我們只需要事先將不超過2^16的所有素數計算出來就可以了。
比如
D=2^a*3^b*5^c*p1^d1*p2^d2*...*pk^dk
其中p1,p2,...,pk是不小於7的素數。
那麼如果(N,D)=1, N/D的迴圈節長度為（結果同N無關）
  L(D) = 3^s(b,2)*(p1-1)p1^(d1-1)*(p2-1)p2^(d2-1)*....*(pk-1)pk^(dk-1)
其中s(b,2)在b>=2時為b-2,不然為0。
不過這個L(D)不一定是N/D的最小迴圈節長度，也可能是L(D)的一個因子。
  只要我們找到一個數x使得10^x=1(mod D),那麼x就是N/D的迴圈節長度了。
當然是否去計算最小迴圈節可能不重要，比如0.(3)的結果寫成0.(33)也是可以接受的。
  而如果一定要計算最小迴圈節，考慮到實際中L(D)的因子數目不會太多，我們窮舉L(D)的因子x判斷是否10^x=1(mod D)也不難。當然這個計算過程中還有很多技巧，這個不細說了。

做到這些後，我不認為餘下還有多少機會，唯一還可以實驗的可能就是每步計算
cur*=10^k之後
我們需要計算
cur/=D;
這一步有可能可以將除法轉化為乘法。
我給一個對於unsigned類型的數，將除法轉化為乘法的程式：
#include <stdio.h>
unsigned int_inv(unsigned x){
    unsigned long long L;
    unsigned long long M=1ULL<<32;
    unsigned w;
    int bits=0;
    L=1ULL;
    while(L<x){L*=2;bits++;}
    do{
      w=(L/x)*x+x-L;
      if(L/w>=M)return (unsigned)(L/x+1);
      bits++;L<<=1;
    }while(bits<64);
    return 0;
}

int main(int argc, char *argv[]){
    unsigned n=atoi(argv[1]);
    printf("%u/n",int_inv(n));
}
上面的程式就是對於任何輸入的n,將輸出一個數字m
那麼對於任意的計算a/n可以轉化為(a*m)>>(bits-32). 其中bits是函數int_inv中最後用的的bits.
不過這個計算最後還是要用彙編實現效率才高。

我試了一下Divid by constant最佳化的作用
為了簡單起見，我準備變更一下演算法，讓最後的迴圈逆序輸出，這樣，我們就可以設計一個演算法讓編譯器自己做最佳化了，
首先我們需要一個能夠計算任何數關於10的離散倒數，這個很簡單：
unsigned inv(unsigned a, unsigned b){
    int s,t;
    a=a%b;
    if(a==1)return 1;
    s=inv(b,a);
    t=(s*b-1)/a;
    return b-t;
}

unsigned inv10(unsigned p){
    return inv(p,10);
}

而且這個代碼效能不重要，我就不最佳化了。
然後我們可以將if(N>0)裡面的代碼替換為：
        int u=inv10(D);
        int w;
        int index[10];
        longlong cur=N;
        for(w=1;w<10;w++){
            index[w]=((10-w)*u)%10;
        }
        printf("(");
        do{
            w = index[cur%10];
            cur=(cur+w*D)/10;
//            printf("%d",w);
        }while(cur!=N);
        printf(")%d",w);

這樣，所有的除法運算就已經變成除以常數10了。
很遺憾，我發現這個代碼計算速度還是同原先代碼幾乎一樣。
稍微分析一下，就可以知道了，主要原因在於cur被聲明為long long,是64為整數，躍出了最佳化的範圍了。
將cur聲明改成int後（這樣，對於D=100000007就不能使用了，乘法要越界了）
結果對於50000017,計算只需要0.7s了（不輸出結果）。同原先4.6秒相比，提高了很多倍。
所以這種最佳化是非常有效，只是對數字範圍有點要求。

上面代碼可以稍微修改一下，就不會有越界問題了：
代碼如下：唯一問題是迴圈節內部的資料是顛倒輸出的：
#include <stdio.h>
#ifdef WIN32
typedef __int64 longlong;
#else
typedef long long longlong;
#endif

unsigned inv(unsigned a, unsigned b){
    int s,t;
    a=a%b;
    if(a==1)return 1;
    s=inv(b,a);
    t=(s*b-1)/a;
    return b-t;
}

unsigned inv10(unsigned p){
    return inv(p,10);
}

void output_intpart(int intpart, int shift){
    int i, v;
    v=1;
    for(i=0;i<shift;i++)v*=10;
    printf("%d.",intpart/v);
    if(shift>0){
       printf("%0*d",shift,intpart);
    }
}

int main(int argc, char *argv[]){
    int N = atoi(argv[1]);
    int D = atoi(argv[2]);
    longlong LN;
    int count2=0,count5=0;
    int shift,i;
    int intpart;
    while(D%2==0){count2++;D/=2;}
    while(D%5==0){count5++;D/=5;}
    if(count2>count5){
        shift=count2;
        LN = N;
        for(i=0;i<count2-count5;i++)LN*=5;
    }else if(count2<count5){
        shift=count5;
        LN = N;
        for(i=0;i<count5-count2;i++)LN*=2;
    }else{
        shift=count2;
        LN = N;
    }
    intpart=LN/D;
    N = LN%D;
    output_intpart(intpart, shift);
    if(N>0){
        int u=inv10(D);
        int w;
        int index[10];
        int cur=N;
        for(w=1;w<10;w++){
            index[w]=((10-w)*u)%10;
        }
        printf("(");
        do{
            w = index[cur%10];
            cur=(cur+w*D)/10;
//            printf("%d",w);
        }while(cur!=N);
        printf(")%d",w);
    }
    printf("/n");
    return 0;
}

最後我們來解決檔案輸入輸出問題：
通過將檔案資料分塊輸出，可以顯著加速寫檔案速度：
比如：
#define BUFFER_LEN (4096)
...
    char buf[BUFFER_LEN];
    int used=0;
    FILE *fout=fopen("out.txt","wb");
    ...
        do{
            int curm10=cur%10;
            int curd10=cur/10;
            w = index[curm10];
            cur=curd10+Dd10*w+(curm10+Dm10*w)/10;
            buf[used++]=(char)('0'+w);
            if(used==BUFFER_LEN){
                fwrite(buf,BUFFER_LEN,1,fout);
                used=0;
            }
//            printf("%d",w);
        }while(cur!=N);
        if(used>0){
            fwrite(buf,used,1,fout);
        }
...
通過這種修改，包含檔案輸入輸出的df3版本對於1000000007的時間從
 2m42s降低到39s.
而我發現在我的電腦上，我簡單把輸出結果用cp命令複製一份也需要花費同樣長的時間，這說明已經沒有任何最佳化機會了。

現在我們就可以設計一個演算法，基本上可以達到極限速度，而且按正常順序輸出結果。方法很簡單。
i)事先計算出迴圈節長度。
ii)使用上面演算法算出迴圈節中每位結果。結果是逆序的，我們可以按逆序順序將每一段(比如4k)資料寫入數組buf中。 根據迴圈節長度可以計算出這組資料應該在檔案中儲存的位置，通過調用fseek將資料寫入檔案指定位置。需要注意的需要調整使得寫入的資料的檔案中起始位移量是4k的倍速。