擂台：分数化小数

[medie2005]

题目描述：

将分数转化为小数，相信很多人都会吧．在计算机中并能直接进行分数运算，需要将分数转换化为浮点数或双精度数才能运算，但这样会导致结果的不精确，那么，这里给定一个分数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/100003 time : 0ms

当N/D=1/1000003 time : 16ms

当N/D=1/10000019 time : 1515ms

当N/D=1/50000017 time : 7500ms

当N/D=1/100000007 time : 15094ms

当N/D=1/1000000007 time : 201594ms

[medie2005]

多谢各位赏脸。

首先，很佩服mathe。其实我的方法就是他说的：

“还有一种可行方法是将循环中每次
cur*=10
改成
cur*=10^k (1<=k<=10)
这样每次可以得到多位而不是一位结果，”
==============================================================

不过呢，我已经用数学方法求得小数的循环节长度，因此mathe所说的：

“当然，在每个循环中，我们需要检查其每一位来判断是否已经出现循环”

在我的方法中是不需要的。

通过测试，正如mathe所说，效率确实只能提高10倍左右。


至于是否还能再优化，我降低一下标准，将原来的10秒改成30秒，如果有人达到30秒之内，就可以了。

[medie2005]

谢谢二位的回复。

不过，我还是得再说明一下：我并不是没想到算法，也不是没写出代码，只是我想看看是否能再降低运行时间（结果的输出占的时间也包括在内）。

这样吧，如果有那位高人的程序在N/D=1/1000000007时，运行时间(包含结果输出时间)低于10秒，100分送上!

[ahuige]

无限不循环小数怎么办?

[mathe]

事先计算循环节长度不难，需要对整数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.
不过这个计算最后还是要用汇编实现效率才高。

[JTZY]

循环节长度即为“指数”：先求欧拉函数，而后求出其约数子集，满足 a^r mod m = 1 最小的 r 即为所求。

[mathe]

To fire_woods
包含printf的结果：
[newtest]$ time ./df 1 50000017 >out

real 0m11.584s
user 0m10.740s
sys 0m0.850s
[newtest]$ time ./df 1 100000007 >out

real 0m23.513s
user 0m22.660s
sys 0m0.510s
[newtest]$ time ./df 1 1000000007 >out

real 3m41.662s
user 3m33.530s
sys 0m7.640s

不包含的时间：
[newtest]$ time ./df2 1 50000017 >out2

real 0m4.605s
user 0m4.580s
sys 0m0.000s
[newtest]$ time ./df2 1 100000007 >out2

real 0m9.235s
user 0m9.220s
sys 0m0.000s
[newtest]$ time ./df2 1 1000000007 >out2

real 1m31.499s
user 1m31.470s
sys 0m0.000s

上面结果用的代码就是上面的代码，只是其中cur和d的声明改成long long(不然最后两个数据计算中要溢出）。而且对于循环中printf被注释的那个版本，d的声明移到循环外面，并且在外面的printf中输出最后的结果d,以防止编译器的优化去掉整个循环（不然循环的输出没有被使用）

[mathe]

我试了一下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秒相比，提高了很多倍。
所以这种优化是非常有效的，只是对数字范围有点要求。

[mathe]

上面代码可以稍微修改一下，就不会有越界问题了，结果如下：
[newtest]$ time ./df3 1 50000017 >out

real 0m0.717s
user 0m0.710s
sys 0m0.010s

[newtest]$ time ./df3 1 100000007 >out

real 0m1.422s
user 0m1.420s
sys 0m0.000s
[newtest]$ time ./df3 1 1000000007 >out

real 0m14.274s
user 0m14.270s
sys 0m0.000s

代码如下：唯一问题是循环节内部的数据是颠倒输出的：
#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;
}

[mathe]

注意上面代码用gcc编译至少要打开-O3编译选项。Windows下我没有试验过，如果大家性能没有达到，可能是因为优化选项打开不够。
上面代码如果将10改成100，10000等，应该还可以提高速度，所以纯粹计算部分提高到10s以内没有任何问题。（当然我计算的循环节是逆序输出的，但是这不是问题。至少证明将除法改变成乘法可以极大提高计算速度）。
问题在于输出部分占用时间太多，我们可能可以通过缓存部分数据在能存然后一次性输出提高部分速度，但是不能改变本质问题。

[mathe]

在贴一个这个代码中如果将所有输出打开所花费的时间：
[newtest]$ time ./df3 1 50000017 >out

real 0m7.668s
user 0m7.390s
sys 0m0.280s
[newtest]$ time ./df3 1 100000007 >out

real 0m15.427s
user 0m14.840s
sys 0m0.580s
[newtest]$ time ./df3 1 1000000007 >out

real 2m42.223s
user 2m35.300s
sys 0m6.630s

[mathe]

试着用100的版本来替换10的版本，结果发现如果用gcc编译，结果同10的版本新能相差无几（还稍微降低了一些）。

我试着用Intel的C/C++编译器来优化（对于10的版本，我试验过速度同gcc差不多，没有优势），
结果发现果然可以提高将近一倍（不包含输入输出）
[newtest]$ time ./df4 1 10000017
0.()00

real 0m0.006s
user 0m0.000s
sys 0m0.000s
[newtest]$ time ./df4 1 100000007
0.()00

real 0m0.756s
user 0m0.770s
sys 0m0.000s
[newtest]$ time ./df4 1 1000000007
0.()00

real 0m7.545s
user 0m7.540s
sys 0m0.000s

[fire_woods]

我测试了一下.
winXp
P4 3G, 512M
VC6.0 O2
1000000007 只要2s左右的时间,当然不包括printf.
输出到文件的话,大概要20多秒.

[mathe]

你的机器输出这么快？难道是我的硬盘太慢了？:(

[mathe]

To fire_woods,你用的是否是100000007,而不是1000000007。我觉得你的计算机不可能这么快。
2s左右。我的可是P43.6G，2G内存（内存不重要）

[fire_woods]

我少了个0.哈哈

[mathe]

通过将文件数据分块输出，可以显著加速写文件速度：
比如：
#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命令复制一份也需要花费同样长的时间，这说明已经没有任何优化机会了。

[mathe]

现在我们就可以设计一个算法，基本上可以达到极限速度，而且按正常顺序输出结果。方法很简单。
i)事先计算出循环节长度。
ii)使用df3中算法算出循环节中每位结果。结果是逆序的，我们可以按逆序顺序将每一段(比如4k)数据写入数组buf中。 根据循环节长度可以计算出这组数据应该在文件中保存的位置，通过调用fseek将数据写入文件指定位置。需要注意的需要调整使得写入的数据的文件中起始偏移量是4k的倍速。

[chiyq118]

同上

[weartoby]

k

[xiaoluomeixi]

这是一道很有意思的题目