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;
}
};
|