TC12389

传送门
我想出了官方题解中第二种算法的一部分，苦于多项式太差，无法进行优化，最后写抄的是第一种方法。有机会还是学习一下多项式相关理论，学习一下第二种方法。
有几个要点要注意：
1.f[0]要视为常数项，而不能作为题目中所谓的n!0来看待，因为n!0并不直接参与递推，所以转移矩阵的第0列实际上是用于存放一些别的东西
2.根据题解里的定义，求出矩阵$A_0,A_1,A_2......$之后，一定要按编号从大到小乘，不可弄反。

#include<cstdio>
#include<cstring>
const int mo=1000000009,K=18,p[4]={2,3,5,7};
typedef long long ll;
inline int pow(int x,int y){
    int ans=1;
    for(;y;y>>=1,x=1ll*x*x%mo)if(y&1)ans=1ll*ans*x%mo;
    return ans;
}
struct matrix{
    int a[K][K];
    static int n,mo;
    inline void ini(){
        for(int i=0;i<=n;++i)memset(a[i],0,(n+1)<<2),a[i][i]=1;
    }
    inline void cle(){
        for(int i=0;i<=n;++i)memset(a[i],0,(n+1)<<2);       
    }
    inline matrix operator*(const matrix&rhs)const{
        static matrix ans;ans.cle();
        for(int i=0;i<=n;++i)
            for(int j=0;j<=n;++j)
                for(int k=0;k<=n;++k)
                    ans.a[i][k]=(ans.a[i][k]+1ll*a[i][j]*rhs.a[j][k])%mo;
        return ans;
    }
}a[100],ans;
int matrix::n=0,matrix::mo=0;
inline matrix pow(matrix x,int y){
    static matrix ans;ans.ini();
    for(;y;){
        if(y&1)ans=ans*x;
        if(y>>=1)x=x*x;
    }
    return ans;
}
inline int calc(int n,int k,int p,int mo){
    if(n<p)return 0;int b[100],i,j,l;a[0].cle();
    ll z=p;
    for(i=1;i<=k;++i)
        for(j=i;j<=k;++j)a[0].a[i][j]=1;a[0].a[0][0]=1;
    matrix::mo=mo;
    for(i=1;z<=n;z*=p,++i){
        a[i]=pow(a[i-1],p);
        for(j=1;j<=k;++j)a[i].a[0][j]=(a[i].a[0][j]+1)%mo;
    }
    for(l=--i,z/=p;i>=0;z/=p)b[i--]=n/z,n%=z;ans.ini();
    for(i=0;i<=l;++i)
        for(j=1;j<=b[i];++j)ans=a[i]*ans;
    return ans.a[0][k];
}
class MegaFactorial{
    public:inline int countTrailingZeros(int n,int k,int b){
        int i,j;
        for(i=3;i>=0;--i)if(b%p[i]==0)break;
        for(j=0;b%p[i]==0;b/=p[i],++j);matrix::n=k;
        if(j==1)return calc(n,k,p[i],mo);
            else return (calc(n,k,p[i],mo)+mo-calc(n,k,p[i],j))%mo*1ll*pow(j,mo-2)%mo;
    }
};