注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

fanhq666的博客

Fan-Fun

 
 
 

日志

 
 

Sqrt(2)前10000位的各种常数优化  

2010-08-18 11:08:32|  分类: 程序 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
我曾听说有人的程序能在1s内计算出2的算术平方根的前10000位,当时感到相当佩服。
现在,我决定挑战一下。
说一下结果吧:在我的机器上,用时109ms计算sqrt(2)的前一万位

先是算法的选择。
1.牛顿迭代法。这个需要计算高精度除法,如果朴素实现,复杂度是O(n*n)(n是位数),如果有更高级的
实现,那么复杂度可以降下来。不过我没有学过高精度除法高级的实现......只能放弃了。
2.逐位尝试法。其实就是手工开平方的方法。
算法如下:
1.   a=1,c=1,n=0
2.   while n<所需位数
3.        c=c*100
4.        b=max{b|b*b+20ab<c}
5.        输出b
6.        c=c-b*b-20ab
7.        a=10a+b
难点在于第4步的实现。其实也没什么难的,b的范围是0~9,挨个试就行了。
我实现了上面的算法,测得运行时间:
4734ms
嗯,惨了点。怎么办呢?优化!

最先想到的优化:转化为10000进制。
这样理论上可以快4倍。
不过,这时第4步的实现就很困难了。
用二分查找来计算b如何?
实践证明:好一点了,但依然很慢。

这时我想了一个非常自然,但是非常重要的优化的思路:
其实当c,a都很大的时候,b约等于c/20000a(是20000乘以a,记住现在是10000进制了!)
而如果知道c和20000a的前几位,就可以把c/20000a计算得八九不离十了。
举例:
c= 1021002...(一共101位,十进制)
20000a=    215698... (一共100位,十进制)
不用知道其他位,我就可以很负责任地说:b=4
这样,基本的框架出炉了:
A.如果c和20000a都至少有5位(10000进制意义下的5位)
用两个long long分别存c和20000a的前3位(10000进制意义下的3位,如果c和20000a的位数不同,就补0)
设它们是c2,a2
那么有(c2-1)/(a2+1)<=b<=(c2-1)/(a2+1)+1
这个证明其实很简单,就不说了。于是,计算(c2-1)/(a2+1),检查b+1是否符合第4步中的等式,如果符合,就b++
B.那如果c和20000a都不够5位呢?
更简单了,把c和20000a转换为long long,二分查找即可。
C.那如果c和20000a一个够5位,一个不够5位怎么办?
呵呵,如果c够20000a不够,那么照套情况A就行,如果20000a够c不够,b一定就是0(当然)

加入这个优化,程序运行如飞:458ms即跑完!

但常数优化向来是让人贪婪的。我还想继续优化,怎么办呢?
只能上“见不得人”的方法了!

1.编译选项加入-O2
-O2可是一个很牛的选项,它能让系统库里的sort加快一倍,让自己写的代码加快数倍
这时,运行时间325ms,并且大部分时间都在输出。
2.将输出所用的时间刨除
方法就是测量时间的时候只计算结果不输出!
这样,运行时间变为180ms。
这和我的结果109ms依然有距离。我是怎么做到的呢?
3.将IE、chrome、360杀毒、DEV C++、资源浏览器统统关掉,再跑一次
哈哈,运行时间是109ms了!

附:代码
#include <cstdio>
//#include <ctime>
using namespace std;
struct Int{
int n;
int a[10100];
};
Int a,c,d,e;
int n;
int Greater(Int &u,Int &v){
if (u.n>v.n)return 1;
if (u.n<v.n)return 0;
for (int i=u.n-1;i>=0;i--){
if (u.a[i]>v.a[i])return 1;
if (u.a[i]<v.a[i])return 0;
}
return 0;
}
void Add(Int &u,int v){
u.a[u.n++]=0;
for (int i=0;i<u.n && v;i++){
v+=u.a[i];
u.a[i]=v%10000;
v/=10000;
}
while (u.n>1 && u.a[u.n-1]==0)u.n--;
}
void Subtract(Int &u,Int &v){
int j=0;
for (int i=0;i<v.n;i++){
j+=u.a[i]-v.a[i];
if (j<0)u.a[i]=j+10000,j=-1;
else u.a[i]=j,j=0;
}
if (j){
for (int i=v.n;i<u.n;i++){
j+=u.a[i];
if (j<0)u.a[i]=j+10000,j=-1;
else{
u.a[i]=j;
break;
}
}
}
while (u.n>1 && u.a[u.n-1]==0)u.n--;
}
void Multi(Int &u,int v){
u.a[u.n++]=0;
int j=0;
for (int i=0;i<u.n;i++){
j+=u.a[i]*v;
u.a[i]=j%10000;
j/=10000;
}
while (u.n>1 && u.a[u.n-1]==0)u.n--;
}
void output(Int u){
if (u.n)printf("%d",u.a[u.n-1]);
for (int i=u.n-2;i>=0;i--)printf("%04d",u.a[i]);
puts("");
}
void Sqrt2(){
int b,j;
a.n=1;a.a[0]=1;n=0;
c.n=1;c.a[0]=1;
while (n<2500){
//c=100000000c
c.n+=2;
for (int i=c.n-1;i>=2;i--)c.a[i]=c.a[i-2];
c.a[1]=c.a[0]=0;
//b=max{b|b^2+20ab<c}
// d=20000a
d.n=a.n+1;
for (int i=1;i<d.n;i++)d.a[i]=a.a[i-1]+a.a[i-1];
d.a[0]=0;
d.a[d.n++]=0;
j=0;
for (int i=0;i<d.n;i++){
j+=d.a[i];
if (j>=10000)d.a[i]=j-10000,j=1;
else d.a[i]=j,j=0;
}
while (d.n>1 && d.a[d.n-1]==0)d.n--;
// find b, and calculate c=c-b^2-db
if (c.n<=4 && d.n<=4){
long long c2=0,d2=0;
for (int i=c.n-1;i>=0;i--)
c2=c2*10000+c.a[i];
for (int i=d.n-1;i>=0;i--)
d2=d2*10000+d.a[i];
int low=0,high=10000;
if (d2 && high>c2/d2+1)high=c2/d2+1;
while (high-low>1){
b=((high+low)>>1);
if (c2>d2*b+b*b)low=b;
else high=b;
}
b=low;
c2=c2-b*b-d2*b;
for (int i=0;i<c.n;i++)
c.a[i]=int(c2%10000),c2/=10000;
while (c.n>1 && c.a[c.n-1]==0)c.n--;
}else if (c.n<d.n){
b=0;
}else{
long long c2=c.a[c.n-1]*100000000LL+
c.a[c.n-2]*10000LL+
c.a[c.n-3];
long long d2=0;
for (int i=c.n-1;i>=c.n-3;i--){
if (i<d.n)d2=d2*10000+d.a[i];
}
b=(c2-1)/(d2+1);
e.n=d.n;
for (int i=0;i<d.n;i++)e.a[i]=d.a[i];
Multi(e,b);
Add(e,b*b);
Subtract(c,e);
Add(d,b+b+1);
if (Greater(c,d)){
b++;
Subtract(c,d);
}
}
//a=10000a+b
a.n++;
for (int i=a.n-1;i>0;i--)a.a[i]=a.a[i-1];
a.a[0]=b;
n++;
}
}
int main()
{
Sqrt2();
output(a);
//printf("time=%d\n",(int)clock());
getchar();
return 0;
}

  评论这张
 
阅读(1569)| 评论(4)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017