代码优化实例:分数化小数

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的倍速。