阶乘和圆周率

wayne

当n比较小的时候,可以brute force. 当n很大的时候, log(n!) 的前面几位数字 是否可以借助LogGamma函数求得. 取了对数之后,至少能降低 问题的数量级,相信效率要高些

liangbch

使用斯特林公式的普通形式不比对数形式慢多少。计算x^n相当于log2(n)次 x×x的复杂度。

wayne

92# liangbch 明白. 还有一个好处就是不用牵涉 浮点运算. 我还是看完前面的帖子再想想辙吧. ===== 我也是才发现我前面都跟了帖.竟然浑然不记得.

liangbch

我的算法和无心人算法的差别，我的记作B，无心人的记作A。
    1. B使用斯特林公式的普通形式，而A使用斯特林公式的对数形式。
    2. B不计算n的阶乘的常用对数，而是直接计算n的阶乘的值，以10亿进制的表示形式。
    3. A的核心部分为计算64bit整数的的常用对数，使用C语言函数库，A受制于double型数的精度，无法提高精度。 B的和核心部分计算一个多精度数与64bit整数的乘积，可通过增加多精度数数的长度来提高精度。
    4. 性能的PK. 主要看 多精度数与64bit整数的乘积 VS log10(double),那个算得更快。...
liangbch 发表于 2012-12-12 21:34

关于斯特林级数（n的阶乘的普通形式，而非对数形式）的相关资料如下，这里记录下来。
$n! ~~ e^-n * n^(n+1/2) * sqrt(2*pi) * ( 1+ (1/12)/n+(1/288)/n^2+(-139/51840)/n^3+(-571/2488320)/n^4+... )$
请参阅 http://mathworld.wolfram.com/StirlingsSeries.html

$n^-k$ 的系数的分子为 http://oeis.org/A001163， 分母是 http://oeis.org/A001164

liangbch

求系数的PARI/GP 程序，根据 http://oeis.org/A001164，http://oeis.org/A001164改写

1. a_1163_1164_f(n)=
   
2. {
   
3.    local(A, m);
   
4.    if (n<1, n==0, A=vector(m=2*n+1, k, 1);
   
5.    for(k=2, m, A[k]=( A[k-1]-sum(i=2, k-1, i*A[i]*A[k+1-i]))/(k+1));
   
6.     print(numerator(A[m]*m!/2^n/n!));
   
7.     print(denominator(A[m]*m!/2^n/n!));
   
8.     print("__________"));
   
9. }
   
10. 
    
11. for (i=1,32,a_1163_1164_f(i));

复制代码

程序计算结果如下
1
12

1
288

-139
51840

-571
2488320

163879
209018880

5246819
75246796800

-534703531
902961561600

-4483131259
86684309913600

432261921612371
514904800886784000

6232523202521089
86504006548979712000

-25834629665134204969
13494625021640835072000

-1579029138854919086429
9716130015581401251840000

746590869962651602203151
116593560186976815022080000

1511513601028097903631961
2798245444487443560529920000

-8849272268392873147705987190261
299692087104605205332754432000000

-142801712490607530608130701097701
57540880724084199423888850944000000

2355444393109967510921431436000087153
13119320805091197468646658015232000000

2346608607351903737647919577082115121863
155857531164483425927522297220956160000000

-2603072187220373277150999431416562396331667
1870290373973801111130267566651473920000000

-73239727426811935976967471475430268695630993
628417565655197173339769902394895237120000000

34856851734234401648335623107688675640839679447003
2601648721812516297626647395914866281676800000000

909773124599542506852275229422593983242880452145053
811714401205505084859513987525438279883161600000000

-1527335577854677023023224272800947125313629267269390501
9740572814466061018314167850305259358597939200000000

-183856455668177802003316143799518064719008299958634826921
14026424852831127866372401704439573476381032448000000000

2583312098861137963745902036370496943872138148651712093816393
1178219687637814740775281743172924172016006725632000000000

5180134290822682443757710427952467581918233549140896702364013
28277272503307553778606761836150180128384161415168000000000

-527550309097873396592733540579928993424142983691519876840948418433873
14613128884259277641708402381685690086746366936130519040000000000

-2114866241537081164613223324215572812504648703648482437460602956015127
701430186444445326802003314320913124163825612934264913920000000000

180394412915538782140015777241228025103785450235726235175126981743099027459
260932029357333661570345232927379682188943128011546547978240000000000

3226140192053936286912811949056082647586604417173687729452086326364208020303641
55891640688340870308367948893044727924871618020073270576939008000000000000

-10218654456520534088469164280902985100842191028132480093114328858063973003580356809
670699688260090443700415386716536735098459416240879246923268096000000000000

-11294270192060551526555825418377569122875449049360883617198277415758886356961061261
8880988975581887254515845120660348492338221235741297614432239616000000000000

郭先抢

95# liangbch


牛人

倪举鹏

In[1]:= 5385973!

Out[1]= 314159123460773280930725870820910993495089529408514282818639392019470668675908814965545330194114\

mathe

无心人 发表于 2008-12-5 13:27
开始记录结果
d = 3, n = 9
d = 31, n = 62

这个应该提交结果到OEIS
https://oeis.org/A328955

无心人

mathe 发表于 2019-11-1 16:54
这个应该提交结果到OEIS
https://oeis.org/draft/A328955

哈，多久的帖子了，提交到OEIS也是算了结了。
另外现在的计算资源是当时的几十倍了~

无心人

尝试了一天，搜索了94亿数字