Topcoder SRM 569 Div1 - MegaFactorial (矩阵)

首先是对于末尾0个数的处理，设最后得到的数中包含$i$的指数为$F(i)$

对于$B=2,3,5,7$的情况，可以直接计算答案$\sum_{i=1}\sum_{j=1}F(j\cdot B^i)$

对于$B$为质因子组合的情况，即$B=6(2\times 3),10(2\times 5)$，因为$F(i)$实际有单调性，可以直接取较大的质因子

对于$B$为质因子次方的情况，即$B=2^2,2^3,3^2$的情况，设$B=p^k$则答案可以表示为

$\begin{aligned} \lfloor \frac{\sum_{i=1}\sum_{j=1}F(j\cdot p^i)} {k}\rfloor \end{aligned} $

由于要取模，实际上要做一点魔改，设模数为$m$，答案可以表示为$c=ak+b$的形式，这个式子求出的是$a$

则$\lfloor \frac{(ak+b)\mod km } {k}\rfloor =a+\lfloor \frac{(b\mod km)} {k}\rfloor =a$

由于$k\leq 3$，扩大模数后可以用unsigned int 存

下面考虑用矩阵求解上式

$nk!$(下面用$f(n,k)$表示)这个东西可以看作从$n$向下的一个递推式

因此考虑以$k$为矩阵元素，求出每个$f(n,k)$被调用的次数

注意这样递推就是反向的了

递推的转移式子是$f(n,k)\rightarrow f(n,k-1),f(n-1,k)$，其中$f(n,k-1)$的转移需要在层内完成

据此构造矩阵即可，注意$f(n,0)$不能向$f(n-1,0)$转移

考虑对于$\sum_{i=1}\sum_{j=1}F(j\cdot B^i)$的每个$i$求解，一共有$\frac{n} {B^i}$个$j$，每个$j$出现的递推层数为等差数列

即$n\mod B^i,n\mod B^i+B^i,n\mod B^i +2 \cdot B^i\cdots $

我们要求的其实是每一层的$f(i,0)$，所以考虑求出每次$B^i$层的转移矩阵

然后是依次累和，把矩阵的转移中$0\rightarrow 0$的转移赋为1即可做到

tips:首项是$n\mod B^i$

一共有$\log _Bn$种不同的$i$，因此复杂度为$O(\log_Bn\log_2 n\cdot k^3)$

当然更优的做法，咕咕咕

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include<bits/stdc++.h>
using namespace std;

#define reg register
typedef long long ll;
typedef unsigned int U;
#define rep(i,a,b) for(reg int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(reg int i=a,i##end=b;i>=i##end;--i)

#define pb push_back
template <class T> inline void cmin(T &a,T b){ ((a>b)&&(a=b)); }
template <class T> inline void cmax(T &a,T b){ ((a<b)&&(a=b)); }

char IO;
int rd(){
	int s=0;
	int f=0;
	while(!isdigit(IO=getchar())) if(IO=='-') f=1;
	do s=(s<<1)+(s<<3)+(IO^'0');
	while(isdigit(IO=getchar()));
	return f?-s:s;
}

const int N=20;


int n,d;
U P=1e9+9;
struct Mat{
	U a[N][N];
	Mat(){ memset(a,0,sizeof a);} 
	void One(){ rep(i,0,d) a[i][i]=1; }
	U* operator [] (int x){ return a[x]; }
	Mat operator * (const Mat &x) const {
		Mat res;
		rep(i,0,d) rep(j,0,d) rep(k,0,d) res.a[i][k]=(res.a[i][k]+1ll*a[i][j]*x.a[j][k])%P;
		return res;
	}
	void Show(){
		rep(i,0,d) { rep(j,0,d) printf("%d ",a[i][j]); puts(""); }
	}
} A,B,C;

Mat qpow(Mat x,int k){
	Mat res; res.One();
	for(;k;k>>=1,x=x*x) if(k&1) res=res*x;
	return res;
}
int Factor(int &x) {
	int p=-1,c=0;
	rep(i,2,x) if(x%i==0) {
		while(x%i==0) c++,x/=i;
		p=i;
		break;
	}
	return x=p,c;
}

class MegaFactorial {
	public:
		int countTrailingZeros(int N, int K, int b) {
			A=Mat(),n=N,d=K;

			if(b==10) b=5;
			if(b==6) b=3;
			int t=Factor(b);
			P*=t;

			drep(i,d,0) {
				A[i][i]=1;
				rep(j,0,d) A[j][i]+=A[j][i+1];
			}
			A[0][0]=0;
			ll ans=0;
			for(ll i=b;i<=n;i*=b) {
				B=qpow(A,i); B[0][0]=1;
				Mat res=qpow(A,n%i)*qpow(B,n/i-1);
				rep(i,0,d) (ans+=res[i][0])%=P;
			}
			P/=t,ans/=t;

			return ans;
		}
};