Topcoder SRM 569 Div1

Topcoder SRM 569 Div1 - MegaFactorial (矩阵)

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

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

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

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

$egin{aligned} lfloor frac{sum_{i=1}sum_{j=1}F(jcdot p^i)}{k} floor end{aligned} $

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

则(lfloor frac{(ak+b)mod km }{k} floor =a+lfloor frac{(bmod km)}{k} floor =a)

由于(kleq 3)，扩大模数后可以用unsigned int 存

下面考虑用矩阵求解上式

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

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

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

递推的转移式子是(f(n,k) ightarrow 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(jcdot B^i))的每个(i)求解，一共有(frac{n}{B^i})个(j)，每个(j)出现的递推层数为等差数列

即$nmod B^i,nmod B^i+Bi,nmod B^i +2 cdot B^icdots $

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

然后是依次累和，把矩阵的转移中(0 ightarrow 0)的转移赋为1即可做到

tips:首项是(nmod B^i)

一共有(log _Bn)种不同的(i)，因此复杂度为(O(log_Bnlog_2 ncdot 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;
		}
};