当我用不超过100ms的时间解决了，用时要14s 的这个问题的时候，我泪流满面!!!!!!!:lol
公式编辑器不太熟，还是百度贴吧的形式吧。
上面的情形对S9的情形处理的很好了。通过用欧拉公式，我想没有比这更好的了。只降到了3秒的运行时间，不错。但我还是不太满意，今天上课的时间突然想起了这事，也就思考了一下，发现还有潜力。
并且很优美。
当然是对不含9的情形，是不是觉得n=8,n=9时计算其仍然要不少的时间呢，所以从这里动刀，看能否将这两项的不含9的计算时间砍到0。
我首先想到的仍然是欧拉公式，这个的确没得说。细看#89楼的search_x()系列的循环，觉得可以用欧拉公式，并且很一致。
想法：
当n,m均较大时，有估计式：ln(n) - ln(m)=1/n+1/(n-1) + ... +1/(m-1)
比如：对于n=8时的情形有：
search_8()
=(ln10000008 - ln9999999)+ (ln10000018 -ln10000009) + ...+ (ln88888888 -ln 88888879)
=ln(10000008/9999999)+ ln(10000018/10000009) +... + ln(88888888/88888879)
=ln(1+9/9999999) + ln(1+9/10000009) +.... +ln(1+9/88888879)
到这里，你一定会有想法了，呵呵。没错，我们知道当x趋于0时，有估计式：ln (1+x) =x
好了，问题得到了简化：
=9/9999999 +9/10000009 + .... + 9/88888879
=9* (1/9999999+1/10000009+... +1/88888879)
是不是有点爽，发现这样一来，分母每次增量变成了10，几乎节省了一个数量级的除法次数。
但不太妙，这还是不好算啊，你很快会想到精度要求不太高的情形下，1/9999999=1/10000000于是有:
=9* (1/10000000+1/10000010 + ... +1/88888880)
嗯，不错了，.....，不是说了分母是以10为增量吗？所以有
=0.9 * (1/1000000 +10000001 +... + 1/88888888)
到此你有何感想，我反正想到了：
=0.9 * search_7()
所以要算：s0和s0，只要算一次search_7():)
s0=0.9*search_7(); //如果精度不太高的要求，还可以往前推（比如推到search6需要75ms）。
s0=0.9*s0 ;
---------------------------------------------------------
注意：不要用ln(10000008 /10000000)来估计ln(10000008/9999999) 这样本身就有误差，还会传递给ln(1+x)=x，误差更大。
---------------------------------------------------------
差不多了，我也没什么更好的办法了,此事，我可以安心结束了。#include <stdio.h>
#include <time.h>
void disp(double s0,double s9,int n)
{
printf( "s0[%d] = %.6f\ns9[%d] = %.6f\ns0/s9 = %.6f\n\n\n", n,s0,n,s9,s0/s9);
}
double calcS(int n)//n=1,2,3,...,
{
int e={1,10,100,1000,10000,100000,1000000,10000000,100000000,1000000000};
int k;
double s=0.0;
for(k=e;k<e;k++)
s+=1.0/k;
return s;
}
//------------------------------
double search_1()
{
double s9,s0=0.0;
int k;
for (k=1;k<=8;k++)
s0+=1.0/k;
return s0;
}
double search_2()
{
double s9,s0=0.0;
int a,b,k=10;
for (b=1;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
double search_3()
{
double s9,s0=0.0;
int a,b,c,k=100;
for (c=1;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
double search_4()
{
double s9,s0=0.0;
int a,b,c,d,k=1000;
for (d=1;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
double search_5()
{
double s9,s0=0.0;
int a,b,c,d,e,k=10000;
for (e=1;e<=8;e++,k+=1000)
for (d=0;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
double search_6()
{
double s9,s0=0.0;
int a,b,c,d,e,f,k=100000;
for (f=1;f<=8;f++,k+=10000)
for (e=0;e<=8;e++,k+=1000)
for (d=0;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
double search_7()
{
double s9,s0=0.0;
int a,b,c,d,e,f,g,k=1000000;
for (g=1;g<=8;g++,k+=100000)
for (f=0;f<=8;f++,k+=10000)
for (e=0;e<=8;e++,k+=1000)
for (d=0;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
//下面的可以省了哦，^_^，哈哈
//---------------------------
/*double search_8()
{
double s9,s0=0.0;
int a,b,c,d,e,f,g,h,k=10000000;
for (h=1;h<=8;h++,k+=1000000)
for (g=0;g<=8;g++,k+=100000)
for (f=0;f<=8;f++,k+=10000)
for (e=0;e<=8;e++,k+=1000)
for (d=0;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}*/
/*double search_9()
{
double s9,s0=0.0;
int a,b,c,d,e,f,g,h,i,k=100000000;
for (i=1;i<=8;i++,k+=10000000)
for (h=0;h<=8;h++,k+=1000000)
for (g=0;g<=8;g++,k+=100000)
for (f=0;f<=8;f++,k+=10000)
for (e=0;e<=8;e++,k+=1000)
for (d=0;d<=8;d++,k+=100)
for (c=0;c<=8;c++,k+=10)
for (b=0;b<=8;b++,k+=1)
for (a=0;a<=8;a++,k++)
s0+=1.0/k;
return s0;
}
*/
//----------------------
int main()
{
doubles0={0};
doubles9={0};
//由欧拉公式lim(1+1/2+...+1/n)=lnn+C(n->无穷,c为常数)
//从第n=8开始，n值以达到其极限的精度要求
//完全可以用估计式ln( n/m )=1/n +1/(n-1)+...+1/(m+1)
//这样明显减少计算时间
doubles_8=2.3025851829940506340183244547094;
doubles_9=2.3025851019940457335179917876844;
doublet0=clock();
s0=search_1();
s0=search_2();
s0=search_3();
s0=search_4();
s0=search_5();
s0=search_6();
s0=search_7();//如果想得到更粗糙一点的近似值时，可以用s0=0.9*s0;输出前几位
s0=0.9*s0;//search_8();
s0=0.9*s0;//search_9();
s9=calcS(1)-s0;
disp(s0,s9,1);
s9=calcS(2)-s0;
disp(s0,s9,2);
s9=calcS(3)-s0;
disp(s0,s9,3);
s9=calcS(4)-s0;
disp(s0,s9,4);
s9=calcS(5)-s0;
disp(s0,s9,5);
s9=calcS(6)-s0;
disp(s0,s9,6);
s9=calcS(7)-s0;
disp(s0,s9,7);
//==================
s9=s_8-s0;
disp(s0,s9,8);
s9=s_9-s0;
disp(s0,s9,9);
printf("Elapsed time: %f s\n",(clock()-t0)/CLOCKS_PER_SEC);
system("Pause");
return 0;
}部分结果：
92# G-Spider
\(^o^)/~，优化无极限啊。看截图你的win7 DPI比较高
本帖最后由 chyanog 于 2013-11-4 16:37 编辑
又写了两个Mathematica版，第一个算9位的需要比较大的内存（2G内存不够）；第二个接近了C的速度，算法和86楼的一致，不同的是只算不含9的数，如果把86楼的calcS去掉（这样也只算不含9的），速度几乎一致，注意：需要Mathematica8或以上版本且有Mathematica能识别的C编译器
f := 1/N@Flatten@Outer & @@ Array, 8] &, n] // Tr;
h := H - H /. H -> HarmonicNumber
Table[{Framed@n, h@n - #, #} &@f@n, {n, 8}] // Grid // AbsoluteTiming
len = 9;
f1 :=Array[{Unique["x"], Boole[# == 1], 8} &, n] /.
it_ :> {1/FromDigits, it} /. {e_, {it__}} :> Hold@Sum
cf = Join @@ Table, {n, len}] /.
Hold :> Compile[{}, {x}, RuntimeOptions -> "Speed",CompilationTarget -> "C"];
h := H - H /. H -> HarmonicNumber
MapIndexed[{Framed @@ #2, #} &,Thread@{Array - #, #} &@cf[]] // Grid // AbsoluteTiming