Topcoder SRM 569 Div1 - MegaFactorial (矩阵)
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 Bi+Bi,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)\)
当然更优的做法,咕咕咕
#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;
}
};