阶乘和圆周率

无心人

朴素的算法

Prelude> let facthi n m = if (n == m) then n else   if (n > m)  then 1 else   (facthi n (n + (div (m - n) 2))) * (facthi (n + (div (m - n) 2) + 1) m)
Prelude> let fact n = fact 1 n
Prelude> let isT a n = (==) a \$ take (length a) \$ show n
Prelude> let test a = head \$ filter ((isT a) . fact) [1..]

无心人

程序执行相当快，计算最后一项，只用几秒
修改start =  log10(3.14159265), end = log10(3.14159266)等参数可继续计算
===========================
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

void main(void)
{
   double start, end, lnf = 0.0, a = 1.0;
   int i = 0, find = 0;
   
   start = log10(3.14159);
   end = log10(3.14160);
   
   while (i < 100000000)
   {
   lnf += log10(a);
   lnf -= floor(lnf);
   if ((lnf >= start) && (lnf <= end))
   {
     find = 1;
     break;
   }
   i ++;
   a += 1.0;
   }
   if (find)
     printf("find: %8.0f  %0.15f", a, lnf);
   else
     printf("Not Find, Stop in ", a, lnf);
}

无心人

测试了Double
10^8次累加10^(-8)
最后结果是1.000000002289867

但是，使用斯特林公式，需要计算的对数值要多
不知道，使用高精度的对数累加和使用斯特林公式
在超过10^8下哪个比较好
另外，20#的结果需要人核验
本来算出个超过10^8的新结果
现在看来，是否正确需要检验呢

mathe

恩，反正都要计算所有的阶乘，直接对数累加应该是不错的方法。不过同样需要提高以下计算精度

无心人

目前看，32位十进制以上精度是必须的
但不知道mpfr的浮点速度如何？
mathe测试过么？

mathe

没试过。我现在机器重装以后，什么工具都没有了，也没有重装过

无心人

验证阶乘前8位数字的程序，适用于100000000!以下
======================
  import Time

  less32 a = until (<10^32) (\`div\` 10) a
  
  facthi :: Integer -> Integer -> Integer
  facthi n m = if (n == m) then (less32 n) else
               if (n > m)  then 1 else
                 less32 \$ (facthi n (n + (div (m - n) 2))) * (facthi (n + (div (m - n) 2) + 1) m)

  facth n = facthi 1 n

  main = do
           t <- getClockTime
           let f = facth 65630447
           print f
           print t
           t <- getClockTime
           print t

无心人

实际计算表明
选择最高位16个数字并不能保证最后结果前8位是正确的
所以，刚刚修改了上面的帖子
选择32位精度

测试程序在我下午5点多点在学校运行呢
结果，周一放上去，也同时确定了65630447!是否缺失以31415626开始
20#已经确认结果均正确
马上要用mpfr写个寻找31415265的

gxqcn

原帖由 无心人 于 2008-12-5 13:27 发表
开始记录结果
d = 3, n = 9
d = 31, n = 62
d = 314, n = 62
d = 3141, n = 10044
d =31415, n = 50583
d = 314159, n = 1490717
d = 3141592, n =   5573998  开始数字是3141592381
d = 31415926, n = 6563 ...

我记得 liangbch 曾写过指定精度的快速阶乘计算器，
可以指定保留任意位的有效精度，
用它来复测可能更好。

无心人

哦？

呵呵
下个候选是1.3亿附近的
再下个可能就要10多亿了
再下个可能就是计算的极限了，大概500亿内