MegaFactorial

作者：谢兴宇

关键词：矩阵乘法 多项式 倍增 伯努利数 数论

题目简述

我们作出如下递归定义：

$f(n,k)=\left \{ \begin{aligned} f(n,k-1)f(n-1,k) & ,n>0,k>0 \\ 1 & ,n=0 \\ n & ,k=0 \end{aligned} \right.$

给出N,K,B，求f(N,K)在B进制下末尾0的个数，答案对$10^9+9$取模。 $1\le N \le 10^9,1\le K \le 16,2\le B \le 10$

初步分析

求B进制末尾的0的数量就是求最大的c，使得$B^c|f(n,k)$。

如果B是质数的话（2,3,5,7），答案就是将f(n,k)质因式分解后B的指数，我们不妨假定我们已经可以处理这种情况了。

如果B不是质数，因为$B\le 10$，所以只有两种情况：

一种是$B=p_1p_2(2 * 3,2 * 5)$，即b为两个不同的质数的乘积，这样的话显然（如果不显然的话，可以通过算法二中的分析知道）答案是f(n,k)中较大的质数的指数，又可以归约为b是质数的情况；

或者$B=p^x（2^2,2^3,3^3）$，这时答案就是f(n,k)中p的指数除以x下取整。考虑模意义下的下取整就比较麻烦了。

对于这个可以有两种处理方式：设f(n,k)中p的指数为c，$m=10^9+9$，那么$\lfloor{c \over x}\rfloor \ mod \ m =(c-(c \ mod \ x)) \ mod \ m =((c \ mod m)-(c \ mod \ x)) \ mod \ m$，所以对于模m和模x的情况计算两遍即可；或者我们也可以直接计算$\lfloor{c\ mod \ xm \over x}\rfloor$，这是因为：

设$c=ax+b,0\le b<x$,

$\lfloor {(ax+b) \ mod \ xm \over x} \rfloor \\ =\lfloor {x(a \ mod \ m)\over x}+{b \ mod \ xm\over x}\rfloor \\ = \lfloor (a \mod m) +{b \over x} \rfloor \\ =(a \ mod \ m)+ \lfloor {b \over x } \rfloor \\ = \lfloor {c \over x} \rfloor \ mod \ m$

因为x最多只有3，所有$(xm)^2$刚好不会爆long long。。

那么我们现在的问题就是求f(n,k)中p的指数$c \ mod \ m$，m中可能有一个因子是2或有一个因子是3。

我们不妨直接改写f，令f(n,k)为原来的f(n,k)中p的指数。那么递归式就可以改写为：

$f(n,k)=\left \{ \begin{aligned} f(n,k-1)+f(n-1,k) & ,n>0,k>0 \\ 0 & ,n=0 \\ g(n) & ,k=0 \end{aligned} \right.$

g(n)为n中p的指数。

注意到递归式的第一行，我们可以这样展开：$f(n,k)=\sum_{i=1}^k f(n-1,i)+f(n,0)$；也可以这样：$f(n,k)=\sum_{i=1}^n f(i,k-1)$。即可以转化成某一维是前缀和的形式。那么根据展开的维的不同，就可以有不同的做法。

算法一

考虑$f(n,k)=\sum_{i=1}^k f(n-1,i)+f(n,0)$，如果我们把$f(n,0)$看成常数项，那么就意味着$f(n,0),f(n,1),...,f(n,k)$都是$f(n-1,0),f(n-1,1),...,f(n-1,k)$的线性表示，而n又这么大，所以我们不妨尝试用矩阵乘来解决这个问题。

设由$F(0)$转移到$F(p^t)$的矩阵是$A_{t}$，那么$A_{t+1}$就是先算出$A_t^p$，然后再给所有线性表示中的常数项+1。

最终的答案矩阵就是将n p进制分解后把矩阵从大到小依次乘起来即可。

时间复杂度$O(k^3 \log_pn \log_2p)=O(k^3 \log_2n)$。

算法二

由伯努利数的相关知识可以知道，一个关于n的k次多项式的1~n的前缀和是关于n的k+1次多项式。所以看到$f(n,k)=\sum_{i=1}^nf(i,k-1)$，我们不妨尝试用多项式做这个题。

我们再仔细看这个递推式$f(n,k)=f(n,k-1)+f(n-1,k)$，考虑$f(i,0)$对$f(n,k)$的贡献的话（即令$f(n,0)=[n=i]$），这实际上就是一个杨辉三角，或者说等于从$(i,1)$出发，只能向上或向右走，到达$f(n,k)$的路径数。所以$f(i,0)$对$f(n,k)$的贡献是$C_{n-i+k-1}^{k-1}$，即i会对$f(n,k)$贡献$g(i)C_{n-i+k-1}^{k-1}$。

所以$f(n,k)=\sum_{i=1}^ng(i)C_{n-i+k-1}^{k-1}=\sum_{i=1}^{\lfloor log_pn \rfloor}\sum_{j=1}^{\lfloor {n \over p^i} \rfloor }C_{n-jp^i+k-1}^{k-1}$。（通过这个式子就可以看出为什么当p增大，$f(n,k)$不会更大）

$C_{n-jp^i+k-1}^{k-1}={\prod_{i=1}^{k-1}(n-jp^i+i) \over (k-1)!}$，也就是说这可以看做是一个k-1次多项式。这里我们先忽略分母中的$(k-1)!$，因为这涉及到逆元的问题，我们先认为分母等于1...

那么我们可以倍增来处理这个式子。我们起初只有$F(x)=\prod_{i=1}^{k-1}(x+i)$，最终要求$\sum_{j=1}^{\lfloor {n\over p^i} \rfloor }F(x+jp^i)$。如果我们有m，那么$F(x,2m)=F(x,m)+F(x+mp^i,m)$。

这样的话时间复杂度是$O(k^2\log_2^2n)$。

但是其实我们没有必要这样做，我们回到这个式子$f(n,k)=\sum_{i=1}^ng(i)C_{n-i+k-1}^{k-1}$，使用类似算法一的思路，会有$F(x,k)=\sum_{i=0}^{p-1}F(x+ip^{k-1},k-1)+F(x,1),F(x,1)=\prod_{i=1}^{k-1}(x+i)$。

这样时间复杂度就是$O(k^2log_2n)$。

最后我们再来处理逆元的问题，实际上有这样的式子：

${ax \over x } \ mod \ m ={x\over x}(a \ mod \ m) = { { ax \ mod \ mx} \over x}$

所以我们可以把$(k-1)!$中的x都乘到模数里，需要的时候再除下来。$(k-1)!$中的x并不会很大，当x=2时只有$2^{10}$，当x=3时只有$3^6$。

注意这里如果使用只求一遍（模$2*(10^9+9)/3*(10^9+9)$）然后直接下取整的做法的话，因为模数又要乘$2^{10}/3^6$，计算过程中可能会爆long long，就需要用到$O(1)$的快速乘。

算法三

我们再回到这个式子。。。

$f(n,k)=\sum_{i=1}^{\lfloor log_pn \rfloor}\sum_{j=1}^{\lfloor {n \over p^i} \rfloor }{\prod_{t=1}^{k-1} (n-jp^i+t)\over (k-1)!}$

注意到我们带进多项式中的是一个等差数列（$\prod_{t=1}^{k-1}(x+t)$），或者更进一步，是首项为0的等差数列（$\prod_{t=1}^{k-1}(x+(n \ mod \ p^i)+t)$）。那么我们不妨把公差提出来：$\sum_{t=0}^{k-1}a_tp^t\sum_{j=0}^{\lfloor {n\over p^i}\rfloor-1}j^t$（$a_t$是多项式的第$t$项）。那么后面就是一个伯努利数！伯努利数不妨用经典的$O(k^2)$递推求出。

时间复杂度$O(k^2log_2n)$。

不过伯努利数在递推的过程中要除一个$k!$，导致模数上要乘的东西更大了。。不过还好，$k!(k-1)!$中的2也只有26个，3只有12个。还是可以接受的。