我在考虑用256位整数计算的可能，按道理，256位整数，连续乘2^64次，误差接近64位吧~~~~
本帖最后由 uk702 于 2021-2-3 08:34 编辑
我个人的看法 ，当 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 于 2021-2-3 08:30 编辑
但当 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。
在逐个 n 进行搜索时，每个 n 基本上只需要做一次乘法，就能基本判断出它是否接近 pi，每个 n 的计算量是极小的。
而通过上述方法，需要通过每间隔 10^5 进行一次校准，以避免误差的积累，校准的时候，必须有高精度的计算才能准确判断误差，同时，一旦我们发现某个候选的 n，也需要通过高精度的计算才能准确判断误差。
如果我们通过斯特林公式或者 loggamma 来计算高精度误差的话，那就需要进行1次的高精度开平方，1次还是2次的求对数，多次乘法和除法，多次加法，每次都要对几十位数进行计算，其计算量可想而知。
因此，如果搜索步长不能做到在很大范围内（我想至少需要 10^4）保持误差稳定可预测的话，这种方法肯定是赔本的买卖，并不能给我们带来收益。
使用 julia 进行初步测试，发现 2 个问题，1 是速度还是太慢不具备实用性，2 是计算过程漏解，参数敏感。
问题1 正在改成 c++，问题2 不好说，搞不好是方案性的错误，如果是这样，这种方案将演变成一种概率解。
本帖最后由 uk702 于 2021-2-3 10:00 编辑
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 的精度就足够了
lgpi=log(big(pi));lg10=log(big(10.));
h(n) = (n=big(n);log(sqrt(2.0*pi*n)) + n*log(n)) - n * (1-1/(0.4+12*n^2));
diff(n) = (n=big(n);delta = h(n) - lgpi; d10f = delta/lg10; d10i = round(BigInt, d10f); d10f-d10i);
M=10^5; N=M; ps=big(0.0002); nt = zeros(Int64, M*10);
function getnext(n, delta, dt, dd)
id = rem(floor(Int, delta*M), M, RoundDown);
# 偏差为 θ 的表已经建立，直接使用
if nt != 0
return nt
end
# 为 θ 建立候选者表
i = 0; j = 1; dta = big(0.0);
while true
delta += dt
dt += dd
i += 1
dta = delta - round(BigInt, delta)
#println("delta = $delta, dta = $dta, abs(dta) < ps = $(abs(dta) < ps)")
if abs(dta) < ps || i > N
nt = i
j = j+1
if j>10 || i > N
nt = i
break
end
end
end
return nt
end
let n=big(119027672349942); i = 0; dta = big(0.0); mindiff = big(100.0)
while true
delta = diff(n); absdt = abs(delta);
# println("test $(n), est. = $dta, real = $delta")
if absdt < 10^-8
println((n, absdt))
end
if absdt < mindiff
mindiff = absdt
end
#如果 delta\dt1\dt2 被约成 0 了，下面的计算就没有意义了
#当 n 充分大时，二次差分的累积应该相当小了，
#如果我们将步长限制在一个合理的范围，应该可以不考虑二次差分的影响
dt1=(h(n+1)-h(n))/lg10
dt2=(h(n+2)-h(n+1))/lg10
dd = dt2-dt1
#println((dt1, dt2, dd))
ds = getnext(n, delta, dt1, dd)
for k=1:10
if ds > 0
delta = diff(n+ds); absdt = abs(delta);
if absdt < 10^-8
println((n+ds, absdt))
end
if absdt < mindiff
mindiff = absdt
end
end
end
n += ds
end
end
println("use time = $(time()-t1)")
虽没找到 d14，总算找到超过 Ickiverar 的更优值：
n=119027672349942，3.141 592 653 589 7 087155083108662556094834692416520748510 Ickiverar 的 d13
n=523195175599421 3.141 592 653 589 7 2751919714105771491041459
本帖最后由 uk702 于 2021-2-5 11:09 编辑
重大成果！！！找到了候选的 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/Log] = 3.141 592 653 589 793 763154362168869644108970257452834353669890485...
构造性方法可以得到任意长度的结果，当然不会是最小结果
比如选择
?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