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 B^i+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;
		}
};