阶乘和圆周率

无心人

我在考虑用256位整数计算的可能，按道理，256位整数，连续乘2^64次，误差接近64位吧~~~~

uk702

我个人的看法 ，当 n>119027672349942 时，逐个 n 进行搜索，如果没有几百个核的话，可行性极低。

而本命题的实质是，如何求得一个 n，使得 n!/pi 最大限度地接近 10^k，求对数的话，即 (log(n!) - log(pi))/log(10) 最大限度地接近某个整数 k。

注意到当 n 比较小时，如 n = 10044，和 n 充分大时，如 n>119027672349942 时，这时 log(n!) 的行为，特别是 log(n!) 的一阶差分 △log(n!) = log(n+1) 呈现较大的差异，
近邻的 n1 和 n2，它们的 log(n1+1)、log(n2 + 1) 之间的差异并不能忽视，也就是说，log(n!) 的每次增长，不能当然地被视作在某个区域内相对稳定。

实际上，当 n 比较小时， log(n!) 的二阶差分，在某个受限区域内，大体上可视作接近常数，因此，这时将 log(n!) 视作关于 n 某个二次函数（而不是线性函数），才有合理的成分。
但如何求解一个二次函数的最接近的整数点，以我的知识，除了遍历搜索外，似没有太好的办法。

uk702

但当 n 充分大时，如 当 n = 119027672349942 时，log(n!) 的每次增长为 log(n+1) = 32.410377122762490395465656202840610407，现在令 n=n+10^10，log(n+1) = 32.410461133310253975650287289788756489，两者之差为 8.401054776358*10-5，也就是说，即使 n 增长了 10^10，log(n!) 的每次增长（一阶差分），也小于 10^-5，
也就是说，它在相对大的空间内极为稳定。

至于 log(n!) 的二阶差分，当  n = 119027672349942 时，是一个 10^-15 的量，完全可以忽略不见。

因此，当 n > 119027672349942 时，如果我们将校准步长设置为 10^5，那么，
在第一个 0 < i <= 10^5 内搜索时，log((n+i)!) 可视作 c11 + c12 * i + o(10^-5)，其中 o < 1.0，当 i < 10^5 时。
在第二个 10^5 < i <= 2*10^5 内搜索时，log((n+i)!) 可视作 c21 + c22 * i + o(10^-5)，其中 o < 1.0，当 i < 10^5 时。
在第三个 2*10^5 < i <= 3*10^5 内搜索时，log((n+i)!) 可视作 c31 + c32 * i + o(10^-5)，其中 o < 1.0，当 i < 10^5 时。
...

大体上，我们可以每 10^5个间距，作一次校准，完全可以用 log((n+i)!)= c1+c2*i 的近似公式进行搜索，在搜索范围内，其误差 < 10^-5。

uk702

在逐个 n 进行搜索时，每个 n 基本上只需要做一次乘法，就能基本判断出它是否接近 pi，每个 n 的计算量是极小的。

而通过上述方法，需要通过每间隔 10^5 进行一次校准，以避免误差的积累，校准的时候，必须有高精度的计算才能准确判断误差，同时，一旦我们发现某个候选的 n，也需要通过高精度的计算才能准确判断误差。

如果我们通过斯特林公式或者 loggamma 来计算高精度误差的话，那就需要进行1次的高精度开平方，1次还是2次的求对数，多次乘法和除法，多次加法，每次都要对几十位数进行计算，其计算量可想而知。

因此，如果搜索步长不能做到在很大范围内（我想至少需要 10^4）保持误差稳定可预测的话，这种方法肯定是赔本的买卖，并不能给我们带来收益。

uk702

使用 julia 进行初步测试，发现 2 个问题，1 是速度还是太慢不具备实用性，2 是计算过程漏解，参数敏感。
问题1 正在改成 c++，问题2 不好说，搞不好是方案性的错误，如果是这样，这种方案将演变成一种概率解。

uk702

julia 的参考代码如下，主要思路如下：
1）在 10^5 内及每找到一个候选者只作校准，也就是指定新的偏差值 θ（代码中是变量 delta ）
2）按照新的 偏差值 θ，查找候选者表，如果这时  θ 对应的候选者表还没有建立，则建立之，
     否则，直接从候选者表获取候选者，再对该候选者计算偏差，如果偏差极小，则这个候选者就是所需的解
     否则，重新在 θ （这时 θ 是一个新值）表中查找下一个候选者。
3）若 n 离当初建立 θ 表的时候已经偏离超过 10^10，原  θ 表应当视作老化，应当重新建立  θ 表。
（这一步 julia 代码没有实现）

这里的 θ 相当 于 #133楼的 c11, c21, c31，代码中的 dt 相当于 #133楼的 c12，c22, c32，...，
但由于 c12，c22, c32 在 10^10 内几乎不变（偏差小于 10^-5），因此，这个参数只用于计算候选者，θ 表认为与之独立。

#现假设分多段来逼近，这样，我们首先选取 s1，使得 (log(u1)/log(10))*s1 接近一个整数
#但实际上肯定不会等于整数，这时有一个偏差  θ
#下一步假设偏差是 θ，u2=u1+s1，x=log(u2)/log(10)，这相当于求解 qx=p-θ
#若将 x 视作常数，可以考虑将 qx+θ 最接近整数的 q 制作成一张 [θ, q] 的表。
# θ 本身误差导致的偏差，公式表明 qx+θ 的偏差只取决于 θ 本身的分辨率，
#因此，取 θ 为 10^-5 的精度就足够了

1. 
   
2. lgpi=log(big(pi));lg10=log(big(10.));
   
3. h(n) = (n=big(n);log(sqrt(2.0*pi*n)) + n*log(n)) - n * (1-1/(0.4+12*n^2));
   
4. diff(n) = (n=big(n);delta = h(n) - lgpi; d10f = delta/lg10; d10i = round(BigInt, d10f); d10f-d10i);
   
5. 
   
6. M=10^5; N=M; ps=big(0.0002); nt = zeros(Int64, M*10);
   
7. function getnext(n, delta, dt, dd)
   
8.   id = rem(floor(Int, delta*M), M, RoundDown);
   
9. 
   
10.   # 偏差为 θ 的表已经建立，直接使用
    
11.   if nt[id*10+1] != 0
    
12.     return nt[id*10+1:id*10+10]
    
13.   end
    
14. 
    
15.   # 为 θ 建立候选者表
    
16.   i = 0; j = 1; dta = big(0.0);
    
17.   while true
    
18.     delta += dt
    
19.     dt += dd
    
20.     i += 1
    
21. 
    
22.     dta = delta - round(BigInt, delta)
    
23.     #println("delta = $delta, dta = $dta, abs(dta) < ps = $(abs(dta) < ps)")
    
24. 
    
25.     if abs(dta) < ps || i > N
    
26.       nt[id*10+j] = i
    
27.       j = j+1
    
28.       if j>10 || i > N
    
29.         nt[id*10+10] = i
    
30.         break
    
31.       end
    
32.     end
    
33.   end
    
34. 
    
35.   return nt[id*10+1:id*10+10]
    
36. end
    
37. 
    
38. let n=big(119027672349942); i = 0; dta = big(0.0); mindiff = big(100.0)
    
39.   while true
    
40.     delta = diff(n); absdt = abs(delta);
    
41.     # println("test $(n), est. = $dta, real = $delta")
    
42. 
    
43.     if absdt < 10^-8
    
44.       println((n, absdt))
    
45.     end
    
46. 
    
47.     if absdt < mindiff
    
48.       mindiff = absdt
    
49.     end
    
50. 
    
51.     #如果 delta\dt1\dt2 被约成 0 了，下面的计算就没有意义了
    
52.     #当 n 充分大时，二次差分的累积应该相当小了，
    
53.     #如果我们将步长限制在一个合理的范围，应该可以不考虑二次差分的影响
    
54.     dt1=(h(n+1)-h(n))/lg10
    
55.     dt2=(h(n+2)-h(n+1))/lg10
    
56.     dd = dt2-dt1
    
57.     #println((dt1, dt2, dd))
    
58. 
    
59.     ds = getnext(n, delta, dt1, dd)
    
60.     for k=1:10
    
61.       if ds[k] > 0
    
62.                delta = diff(n+ds[k]); absdt = abs(delta);
    
63.                if absdt < 10^-8
    
64.             println((n+ds[k], absdt))
    
65.          end
    
66. 
    
67.          if absdt < mindiff
    
68.             mindiff = absdt
    
69.          end
    
70.       end
    
71.     end
    
72. 
    
73.     n += ds[10]
    
74.   end
    
75. end
    
76. 
    
77. println("use time = $(time()-t1)")

复制代码

uk702

虽没找到 d14，总算找到超过 Ickiverar 的更优值：
n=119027672349942，3.141 592 653 589 7 087155083108662556094834692416520748510                Ickiverar 的 d13
n=523195175599421   3.141 592 653 589 7 2751919714105771491041459                                       

uk702

重大成果！！！找到了候选的 d15！！！

原 Ickiverar 找到了 d13，也就是满足小数点后 13 位的 n，现在不仅找到了满足小数点后 14 位的 n，还找到了第 15 位的 n. （但无法肯定这个 n 是满足要求的最小值）
Pi                          = 3.141 592 653 589 793 238 462 643 383 28                                                                                                Pi 的真实值
927877657573029! = 3.141 592 653 589 793 763154362168869644108970257452834353669890 * 10^13485028080082190            候选 d15
119027672349942! = 3.141 592 653 589 7 087155083108662556094834692416520748510                                                            Ickiverar 的 d13

Mathematica 代码验证：
10^FractionalPart[LogGamma[927877657573030`200]/Log[10]] = 3.141 592 653 589 793 763154362168869644108970257452834353669890485...

mathe

构造性方法可以得到任意长度的结果，当然不会是最小结果
比如选择
?lngamma(10^50)/log(10)
%26 = 4956570551809674817234887108108339491770560299419608.742478488800740938317582704430223027715000259746
? log(Pi)/log(10)
%28 = 0.4971498726941338543512682882908988736516783243804424461340534999249471120895526746555473864642912224
? 1+%28-.742478488800740938317582704430223027715000259746
%29 = 0.7546713838933929160336855838606758459366780646344424461340534999249471120895526746555473864642912224
? %*log(10)
%30 = 1.737695078662113277768684175517271964338182994044241372599195634164832766833857233716733633608547926
? x*(x+1)-2*%30*10^50
%36 = x^2 + x - 347539015732422655553736835103454392867636598808848.2745198391268329665533667714467433467267217095853
? polroots(%36)
%37 = [-18642398336384260823080966.41470932110761512429739385891105693533660359074166843711436414603488164818 + 0.E-115*I, 18642398336384260823080965.41470932110761512429739385891105693533660359074166843711436414603488164818 + 0.E-115*I]
? lngamma(10^50+18642398336384260823080966)/log(10)
%38 = 4956570551809674817234888040228256310983601453467909.497149872694133854351268207818544726381377479917
? 10^.497149872694133854351268207818544726381377479917
%39 = 3.141592653589793238462642801159841823296283273716797887197076777632859593592539171217958151957586201
? lngamma(10^50+18642398336384260823080967)/log(10)
%40 = 4956570551809674817234888040228256310983601453467959.497149872694133854351268288781451995721840878597
? 10^.497149872694133854351268288781451995721840878597
%41 = 3.141592653589793238462643386828058412090052015250231730023053334962070746313444107907787538049767930
? Pi
%42 = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068