固定小数点の bignumber ライブラリがあり、精度を落とさずに高速階乗を実装したいと考えています。
紙の上でいくつかの数学のトリックを行った後、次の式を取得しました。
(4N)!=((2N)!).((2N)!).{ (2N+1).(2N+3).(2N+5)...(4N-1) }.(2^N)/(N!)
これはすでにかなり高速であり、いくつかのプログラミングのトリックを使用すると、複雑さは~ O(log(n)).
明確にするために、私の現在の実装は次のとおりです。
//---------------------------------------------------------------------------
longnum fact(const DWORD &x,longnum &h) // h return (x>>1)! to speed up computation
{
if (x==0) { h=1; return 1; }
if (x==1) { h=1; return 1; }
if (x==2) { h=1; return 2; }
if (x==3) { h=1; return 6; }
if (x==4) { h=2; return 24; }
int N4,N2,N,i; longnum c,q;
N=(x>>2);
N2=N<<1;
N4=N<<2;
h=fact(N2,q); // get 2N! and N!
c=h*h; for (i=(N2+1)|1;i<=N4;i+=2) c*=i; c/=q; // c= ((2N!)^2)*T1 / N!
for (i=N4+1;i<=x;i++) c*=i; c.round(); c<<=N ; // convert 4N! -> x!, cut off precision losses
for (i=(N2+1)|1,N2=x>>1;i<=N2;i++) h*=i; h.round(); // convert 2N! -> (x/2)!, cut off precision losses
return c;
}
//---------------------------------------------------------------------------
longnum fact(const DWORD &x)
{
longnum tmp;
return fact(x,tmp);
}
//---------------------------------------------------------------------------
今私の質問：
この用語からすばやく取得する方法 はありますか?
N!T1 = { (2N+1).(2N+3).(2N+5)...(4N-1) }すでに回答済み。
明確にするために、この未知の用語を抽出する必要があります。
T2 = (4N)! / (((2N)!).((2N)!))
それで：
(4N)! = (((2N)!).((2N)!)).T2
.../(N!)階乗を計算する必要がないため、これは非常に役立ちます。
項は常に次のT1ように整数分解可能です。
T1 = T2 * N!
最後に、それは私を襲った :) 階乗の素数分解のための小さなプログラムを実行したところ、突然すべてがより明確になりました。
4! = 2!.2!.(2^1).(3^1) = 24
8! = 4!.4!.(2^1).(5^1).(7^1) = 40320
12! = 6!.6!.(2^2).(3^1).(7^1).(11^1) = 479001600
16! = 8!.8!.(2^1).(3^2).(5^1).(11^1).(13^1) = 20922789888000
20! = 10!.10!.(2^2).(11^1).(13^1).(17^1).(19^1) = 2432902008176640000
24! = 12!.12!.(2^2).(7^1).(13^1).(17^1).(19^1).(23^1) = 620448401733239439360000
28! = 14!.14!.(2^3).(3^3).(5^2).(17^1).(19^1).(23^1) = 304888344611713860501504000000
32! = 16!.16!.(2^1).(3^2).(5^1).(17^1).(19^1).(23^1).(29^1).(31^1) = 263130836933693530167218012160000000
36! = 18!.18!.(2^2).(3^1).(5^2).(7^1).(11^1).(19^1).(23^1).(29^1).(31^1) = 371993326789901217467999448150835200000000
40! = 20!.20!.(2^2).(3^2).(5^1).(7^1).(11^1).(13^1).(23^1).(29^1).(31^1).(37^1) = 815915283247897734345611269596115894272000000000
項の素数指数を分析したT2後 (半階乗の後の残り ^ 2)、それらの式を導き出します。
T2(4N) = multiplication(i=2,3,5,7,11,13,17,...) of ( i ^ sum(j=1,2,3,4,5,...) of (4N/(i^j))-(2N/(i^j)) )
- 掛け算はすべてを通して
primes <= 4N - 総和はどこまで
i^j <= 4N
問題は、除算4N/(i^j)と除算を整数演算2N/(i^j) で行う必要があるため、簡単に単純化できないことです。
だから私は別の質問があります：
どうすればこれを計算できますか:
exponent(i) = sum(j=1,2,3,4,5,...) of (N/(i^j))効果的に?iは任意の素数i<=Nです。簡単なはずです。
ここで、次のように項内eの素数の指数を計算します(ただし、これは私の好みには複雑すぎます)。iT2(N)
for (e=0,a=N/i,b=(N>>1)/i;(a)||(b);e+=a-b-b,a/=i,b/=i);
...実装T2しfact(x)て速度を比較してみます...