无心人 发表于 2021-2-2 15:55:55

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

uk702 发表于 2021-2-3 07:58:52

本帖最后由 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:29:28

本帖最后由 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。

uk702 发表于 2021-2-3 08:48:59

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

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

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

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

uk702 发表于 2021-2-3 09:01:40

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

uk702 发表于 2021-2-3 09:12:25

本帖最后由 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)")

uk702 发表于 2021-2-4 16:42:33

虽没找到 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:02:35

本帖最后由 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...

mathe 发表于 2021-4-17 15:05:36

构造性方法可以得到任意长度的结果,当然不会是最小结果
比如选择
?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
页: 4 5 6 7 8 9 10 11 12 13 [14]
查看完整版本: 阶乘和圆周率