WNJXYK
Thanks to the cruel world.
WNJXYKのBlog
快速计算取特殊模数双阶乘
快速计算取特殊模数双阶乘

10^{18}的数量级已经没有办法使用暴力算法来求阶乘了,按照质因子统计也没法统计,束手无策。

分块

第一种方法的大致思路为使用分块的思路,将阶乘的所有数字按照每2^{16}个数字分成一个块,每个块内部快速计算,块与块之间的乘积也快速计算,然后再暴力将最后不能形成一整块的数字统计进入答案中。

考虑我们的问题,求:

1 \cdot 2 \cdot 3 \cdots (2m-1)​

然后,我们每2^{16}元素就分成一个块:

block = [(k\cdot 2^{17}+1)(k\cdot 2^{17}+3)\cdots (k\cdot 2^{17}+2^{17}-1)]

我们令进行换元q=2k+1,则:

block = [(q\cdot 2^{16}-(2^{16}-1))(q\cdot 2^{16}-(2^{16}-3))\cdots (q\cdot 2^{16}+(2^{16}-3))(q\cdot 2^{16}+(2^{16}-1))]

那么就可以首尾两两配对,使用平方差公式(a-x)(a+x)=a^2-x^2进行化简,最终可以得到:

block = [(q^2\cdot 2^{32}-1^2)(q^2\cdot 2^{32}-1^2)\cdots (q^2 \cdot 2^{32} – (2^{16}-1)^2)]

我们把块展开,那么其中包含的因子2^{64}的项在模意义下全部都等于了0,最终每个块的形式肯定为A+B\cdot q^2\cdot 2^{32}的形式,其中:

A = \prod_{i=1}^{2^{16}-1}(-i^2)

B = \sum_{i = 1}^{ 2^{16} – 1 }[ \prod_{j=1~and~j\neq i}^{2^{16}-1}(-i^2)]

经过计算得到:
A=7083242339890757633
B=556099689942450176

这样一来,一个块内部的乘积,我们就可以快速计算了,那么答案就是:

(2m-1)!! = block[k=0]\cdot block[k=1] \cdot \cdots \cdot block[k=\cdots] \cdot [incomplete~block]

最后一个不完整的块,我们可以使用暴力计算。块与块之间的乘积,我们可以使用一下数学公式进行计算:

1^2+3^2 +5^2 +\cdots +(2k-1)^2 = \frac{k(4k^2-1)}{3}

同样的道理,对于包含有因子2^{64}的项,在模运算下就可以忽略了。因此,假设我们的整块k=0\cdots p-1,最终答案应该为:
A^p + B\cdot A^{p-1}\cdot \frac{p(4p^2-1)}{3}\cdot 2^{32}

TIP : 除法部分可以直接将平方差公式分解成两项,其与中一定有一项是的倍数,直接除即可。

代码如下:

#define ull unsigned long long
const ull A=7083242339890757633;
const ull B=556099689942450176;
const ull S=1<<16;

inline ull quickPow(ull x,ull t){
    ull ret=1;
    for (;t;t/=2,x*=x) if (t&1) ret*=x;
    return ret;
}

inline ull solve(ull x){
    ull p=(x+1)/S/2,sumQ;
    if (p%3==0) sumQ=p/3*(4*p*p-1); else{
        if ((p*2-1)%3==0) sumQ=(p*2-1)/3*(p*2+1)*p; 
        else{
            sumQ=(p*2+1)/3*(p*2-1)*p;
        }
    }
    ull ret=quickPow(A,p)+B*quickPow(A,p-1)*sumQ*S*S;
    for (ull i=2ull*S*p+1;i<=x;i+=2) ret=ret*i;
    return ret;
}

构造多项式

我们构造多项式:P_k(x) = \prod_{i=0}^{k-1}(2x+2i+1)

那么(2k-1)!! = P_k(0),只要能计算出这个多项式的值就能求得结果。

定义中,每一个乘法部分里x一次项的系数均为2,所以在模2^{64}的意义下,x^i,(i\geq 64)的项在x\neq \pm 1的情况下是等于0的。所以我们的计算中,只需要时刻维护这个多项式的前64项即可。

根据多项式的定义 ,我们可以推知一下三条性质:

P_{2k}(x)=\prod_{i=0}^{2k}(2x+2i+1)=\prod_{i=0}^{k}(2x+2i+1) \cdot\prod_{i=0}^{k}(2(x+k)+2i+1) = P_k(x)\cdot P_k(x+k)
P_{m+n-1}(x)=\prod_{i=0}^{m+n-1}(2x+2i+1)=\prod_{i=0}^{m-1}(2x+2i+1) \cdot\prod_{i=0}^{n-1}(2(x+m)+2i+1) = P_m(x)\cdot P_n(x+m)
P_{k+1}(x)=P_k(x)\cdot (2x+2k+1)

通过这些性质,我们可以使用快速幂加速多项式计算的速度。

特殊性质

经过尝试我们可以发现,分段阶乘竟然存在同余关系。

\prod_{k=1}^{2^{21}}(2k-1) \equiv \prod_{k=2^{21}+1}^{2^{21}\cdot 2}(2k-1)\equiv \cdots (MOD~2^{64})

所以,我们只需要计算计算一下其中一块的值,然后快速幂+暴力就可以求得答案,非常简便。

赞赏
https://secure.gravatar.com/avatar/f83b57c055136369e9feba5d6671d6b5?s=256&r=g

WNJXYK

文章作者

一个蒟蒻

推荐文章

发表评论

textsms
account_circle
email

WNJXYKのBlog

快速计算取特殊模数双阶乘
$$10^{18}$$的数量级已经没有办法使用暴力算法来求阶乘了,按照质因子统计也没法统计,束手无策。 分块 第一种方法的大致思路为使用分块的思路,将阶乘的所有数字按照每$$2^{16}$$个数…
扫描二维码继续阅读
2018-07-18
<--! http2https -->