阶乘和圆周率

Ickiverar

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

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

liangbch

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

mathe

现在给出的算法都是每次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

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

liangbch

说说我的思路和算法，或者对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

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

liangbch

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

liangbch

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

Ickiverar

我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$ 几乎只能有个位数。

所以这个方法实际上非常不实用。