Ickiverar 发表于 2019-12-18 08:36:49

13位小数的结果有了,119万亿左右。
d =0, n = 9
d =1, n = 62
d =2, n = 62
d =3, n = 10044
d =4, n = 50583
d =5, n = 1490717
d =6, n = 5573998
d =7, n = 65630447
d =8, n = 688395641
d =9, n = 5777940569
d = 10, n = 77773146302
d = 11, n = 1154318938997
d = 12, n = 1544607046599
d = 13, n = 119027672349942
done for n <      126438000000000, maxerr = 0.119045

liangbch 发表于 2019-12-19 09:23:46

我终于想到一个满意的算法,这个算法具有简单、高效、可靠的特点。它可能是终级算法。

liangbch 发表于 2019-12-20 14:44:13

完成了代码,速度和Ickiverar并不多,在我的电脑上,1分钟可搜索200亿左右。代码见https://www.bccn.net/paste/3202

mathe 发表于 2019-12-20 17:51:56

现在给出的算法都是每次n加一来验算,已经计算到$n=10^{15}$左右,继续用这样的算法计算下去效率会很低。
我们能否找到一个算法一次排除一个较大范围的n呢?

假设现在已经计算出某个$\lg(n!)$的小数部分$t_n$,而我们期望找到小数部分对数为$\lg(\pi)$的阶乘.
而我们接下去希望能够批量判断$\lg((n+1)!),\lg((n+2)!),\dots,\lg((n+d)!)$的小数部分是否同$\lg(\pi)$充分接近。
对于$n$远远大于d时,由于$\lg(n)\le \lg(n+h)\le\lg(n+d)$(对于一切$1\le h\le d$),于是我们就可以估计出
$ \lg(n!)+h\lg(n)\le \lg((n+h)!)\le \lg(n!)+h\lg(n)+h\lg(1+\frac{d}{n})\le\lg(n!)+h\lg(n)+\frac{hd}{n}\le\lg(n!)+h\lg(n)+\frac{d^2}{n}$

mathe 发表于 2019-12-20 18:22:20

现在如果我们能够找到$\lg(n)$充分好的分数逼近比如连分数逼近u/v,并且使得d倍误差加d^2/n不超过1/v,那么就可以保证每个值落在一段长度不超过1/v的区间。然后就可以使用中国剩余定理寻找合适的h有可能产生比较接近的值了

liangbch 发表于 2019-12-22 21:00:15

说说我的思路和算法,或者对114#的解读,如有不对,请指正。

对于一个实数x,为了方便起见,我们总是假定 $0<x<1$
如果x恰好能表示为$p/q$(p,q是整数,p与q互质)的形式,当n从1,2,3,
依次递增时,nx的小数部分呈现周期性,其周期为q的整数倍。q称为最小周期。
但对于不能表示为$p/q$形式的实数x,nx的小数部分不会呈现周期性。
但粗略的看,nx的分布仍然大体呈现周期性,$1/x$可称为最小周期。


由于周期不能为小数,我们定义 $p1=\lfloor 1/x \rfloor$, $p2= \ceil {1/x} $,
分别称之为最小欠周期和最小过周期。$n*p1$,可表示为一个整数加一个略小于0的数的和,
而$n*p2$可表示为一个整数加一个略大于0的数的和。如对于无理数$x=\lg(2)$,
$3x \approx 1+(-0.1)$
$4x \approx1+(0.2)$
进一步,使用连分数法,可求得x的的一系列渐进分数,如某个渐进分数$p/q>x$时,
q可称为x的欠周期,如某个渐进分数p/q<x时,q可称为x的过周期。例如$\lg(2)$的
渐进分数依次为$1/3$,$3/10$,$28/93$,$59/196$,$146/485$,$643/2136$.
$\lg(2)$乘以这些分数的分母,其值依次为
$\lg(2)*3   \approx 1 -0.0969100$,3为欠周期
$\lg(2)*10\approx 3 +0.0102999$, 10为过周期
$\lg(2)*93\approx 28 - 0.0042104$, 93为欠周期
$\lg(2)*196 \approx 59 + 0.0018791$, 196为过周期
$\lg(2)*485 \approx 146- 0.0004521$, 485为欠周期
$\lg(2)*2136 \approx 643+ 0.0000707$, 2136为过周期

指数形式
$2^3   = (1-0.2)*10^1$,    3为欠周期
$2^{10}   = (1+0.024)*10^3$,  10为过周期
$2^{93}\approx (1-0.01) * 10^{28}$, 93为欠周期
$2^{196} \approx (1+0.004) * 10^{59}$, 196为过周期
$2^{485} \approx (1-0.0011) * 10^{146}$,485为欠周期
$2^{2136}\approx (1+0.00016) * 10^{643}$,2136为过周期

类似的,若$\lg(n!)$的小数部分非常接近于$\lg(pi)$, 当n很大时, 
$\lg((n+d)!) \approx \lg(n!) + d*\lg(n)$,我们求得$\lg(n)$的斩近分数$p/q$,
q为$\ln(n)$的欠周期或者过周期,则,$\ln((n+q)!)$的小数部分会非常接近于$\lg(n!)$。
利用这一思路,若$n!$的10进制表示的前几位数字略小于$pi$,则我们不需要逐个检查
$(n+1)!$,$(n+2)!$,而只检查$(n+q)!$, q为$\lg(n)$的渐进分数的分母,大大降低了工作量。

上面假定,$\lg((n+d)!) \approx \lg(n!) * d*\lg(n)$,在实际计算过程中,我们需要做更精确的估计
 $\lg((n+d)!) = \lg(n!) + d*\lg(n)+ \lg(1+1/n) + \lg(1+2/n) + ... \lg(1+d/n)$
 $= \lg(n!) + d* \lg(n)+ ( \ln(1+1/n) + \ln(1+2/n) + ... \ln(1+d/n) )/ln(10) $
$= \lg(n!) + d* \lg(n)+ ( 1/n + 2/n + ... + d/n)/\ln(10) -0.5/\ln(10)( (1^2)/(n^2) +(2^2)/(n^2) + ...(d^2)/(n^2) ) + ... $
$= \lg(n!) + d*\lg(n)+ \frac{d(d+1)}{2n*\ln(10)} -0.5/\ln(10) ( {1^2}/{n^2} +{2^2}/{n^2} + ...{d^2}/{n^2} ) + ... $   
$=> \lg(n!) + d*\lg(n)+ \frac{d(d+1)}{2n*\ln(10)} > \lg(n+d)! $


问题描述
1. 已知 $\lg(n!)$的小数部分略小于$\lg(\pi)$
2. 求d, $\lg (n+d)!$的小数部分同样略小于$\lg(\pi)$, 并且对于 $n <i < n+d$, 不存在$\lg((n+d)!) < \lg(i) < \lg(\pi)$

步骤
1. 求得$b=\lg(n!)$和$x=\lg(n)$
2. 求得x的渐进分数序列,多个,具体计算多少个,根据需要做出限制。 
3. 对每一个渐进分数$p/q$, 求 $y= b + q*lg(n)+ \frac{q(q+1)}{2n*ln(10)}$,由上面的分析,
  y总是大于$\lg (n+q)! $, 若y的小数部分小于$\lg(\pi)$, 则$\lg (n+q)!$的小数部分一定小于$\lg(\pi)$
  对于每个q,求得一个y,在所有的y中,若q为d时,对应的y的小数部分取到最大,d即为答案。
4. 取n为n+d, 检查$\lg(n!)$的小数部分是否充分接近$\lg(\pi)$,并继续下一次迭代。

mathe 发表于 2019-12-23 11:07:56

分析得很好,关键是q得选择,不能太大,也不能太小。
对于给定的n,使用连分数,合适的q不一定存在,有可能我们需要寻找前后其它更好的n.所以这时n!不一定总能满足前面部分和$\pi$都一致的条件。

liangbch 发表于 2019-12-23 11:35:03

说中了关键,q选的太小,没有效率,q选的太大,可能会漏掉某些候选解。我其实已经意识到这个问题,随后补上一些修补方法。

liangbch 发表于 2019-12-23 12:05:56

使用连分数法,q增长的很快,但可能越过(漏掉)某些候选解,而使用法雷序列构造法,q则增长地很慢。适当地结合连分数法和法雷序列构造法,是一个考虑的方向。

Ickiverar 发表于 2019-12-23 18:30:43

我12月21日凌晨写的代码,就是用连分数展开近似loggamma的增长率,然后跳过不可能的n值。
实测效率极其低下。
原理是:
有公式$ ln(b!)-ln(a!) = (1/{2a}+lna) (b-a)+(2a-1)/(4a^2)(b-a)^2-...$,
如果${ln(a!)-ln(pi)}/ln10=F$已经非常接近整数,设其小数部分为$F-\text{round(F)}=f$接近0,那我们要连分数展开 $r=1/ln10(1/{2a}+lna)\approx p/q\approx {p'}/{q'}$,其中$p'q'$是比$pq$更不精确一级的连分数。连分数展开的误差是$epsilon=|p'-q' r|$,也就是说对任何比$q$小的整数$q_0$,$r q_0$都不可能在$epsilon$的距离内接近一个整数,也就是说

只要 $f+epsilon>0$ 且 $f-epsilon<0$,对所有 $q_0<q$,$f+r q_0$ 都不可能是整数。
只要 $f+epsilon>delta$ 且 $f-epsilon<-delta$,对所有 $q_0<q$,$f+r q_0$ 都不可能在$delta$的距离内接近整数。实际中,$delta$和想要阶乘跟pi前几位一样有关。

但由于上面那个对数阶乘的公式还有一个过估计的平方项,所以判据就要变成:
$f+epsilon>delta$ 且 $f-epsilon+q^2/(2a)<-delta$
分析一下这个公式:随着$q$的增大,$epsilon$是一直在减小的,$q^2/(2a)$是一直增大的,整个判据越来越不成立。这是自然的,所以$q$实际可取到的范围非常低。
实测中,从 d=12 的结果 a = 1544607046599 开算,$f=2.165\times 10^-14$,基本是0;根据连分数展开和判据可以取到$q=1091$,此时 $epsilon=4.13\times 10^-7$,$q^2/(2a)=3.85\times 10^-7$,这两个数都远大于$f$,可见限制 $q$ 的主要是那个平方偏离项,使得我们无法准确用连分数展开估计阶乘接近整数的程度。
而且一旦跳过了 $q$ 个数,下一个阶乘值的 $f$ 就会变得严重偏离0,达到 $0.000451$ ,此时 $f$ 的精度就严重限制了 $q$,现在 $q$ 只有37。进一步计算下去,之后 $q$ 几乎只能有个位数。

所以这个方法实际上非常不实用。
页: 2 3 4 5 6 7 8 9 10 11 [12] 13 14
查看完整版本: 阶乘和圆周率