阶乘和圆周率

无心人

朴素的算法 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 #include #include 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亿内