LOJ6696 复读机 加强版

Link
考虑单个元素的EGF，显然为\(f(x)=\sum\limits_{i\ge0}[d|i]\frac{x^i}{i!}=\frac1d\sum\limits_{i=0}^{d-1}\exp(\omega_d^ix)\)。
那么答案就是\([x^n]f(x)^k\)。
考虑分圆多项式\(\Phi_d(x)=\prod\limits_{k\in[0,d)\wedge\gcd(d,k)=1}(x-\omega_d^k)=\prod\limits_{k|d}(x^{\frac dk}-1)^{\mu(k)}\)。
那么\(\omega_d^k=(x^k\bmod\Phi_d(x))\mid_{x=\omega}\)。
很显然\(\deg(\Phi_d(x))=\varphi(d)\)，且\(\Phi_d(x)\)的系数都是整数，因此\(\omega_d^k\)可以由\(\omega_d^0,\cdots,\omega_d^{\varphi(d)-1}\)以整数线性表示。
实际上可以证明，对于域\(\mathbb N\)上的线性空间\(V=\operatorname{span}(\omega_d^0,\cdots,\omega_d^{d-1})\)一定有\(\dim V=\varphi(d)\)。
这样我们就可以把\(f(x)\)写成\(F(\exp x,\exp(\omega x),\cdots,\exp(\omega_d^{d-1}x))\)的形式了。
设\(G=F^k\)，那么\(F\frac{\partial G}{\partial x}=kG\frac{\partial F}{\partial x}\)。
由这个式子列出递推式然后解得系数是trivial的。
时间复杂度为\(O(k^{\varphi(d)}\log n)\)，\(\log n\)是快速幂的。
当\(d=2\)时可以通过线性筛幂函数做到\(O(k)\)。

#include<cstdio>
using i64=long long;
const int N=2007,P=19491001;
int n,k,d;i64 w,ans,fac[N],ifac[N],inv[N];
int abs(int x){return x<0? -x:x;}
void inc(i64&a,i64 b){a+=b-P,a+=a>>63&P;}
i64 pow(i64 a,i64 b){i64 c=1;for(a%=P;b;b>>=1,a=a*a%P)if(b&1)c=c*a%P;return c;}
i64 binom(int n,int m){return fac[n]*ifac[m]%P*ifac[n-m]%P;}
void init(int n=2000)
{
    fac[0]=ifac[0]=inv[0]=fac[1]=ifac[1]=inv[1]=1,w=pow(7,(P-1)/d);
    for(int i=2;i<=n;++i) fac[i]=i*fac[i-1]%P,inv[i]=(P-P/i)*inv[P%i]%P,ifac[i]=inv[i]*ifac[i-1]%P;
}
namespace Sub1
{
    void solve()
    {
	for(int a=0;a<=k;++a) for(int b=2-((k^a)&1);a+b<=k;b+=2) inc(ans,pow(a*w+b,n)*binom(k,(k+a+b)/2)%P*binom(k,(k+abs(a-b))/2)%P);
	printf("%lld",4*pow(d,P-k-1)*ans%P);
    }
}
namespace Sub2
{
    i64 a[N][N],b[N][N];
    void solve()
    {
	for(int i=0;i<=k;++i) a[0][i]=a[i][0]=binom(k,i),b[i][0]=i*a[i][0]%P;
	for(int i=1;i<=k+k;++i)
	    for(int j=1;j<=i&&j<=k;++j)
	    {
		i64 x=k*a[i-1][j]-b[i-1][j]-b[i][j-1];
		if(j>=2) x+=k*a[i-1][j-2]-b[i-1][j-2];
		if(i>=2) x+=2*k*a[i-2][j-1]-b[i-2][j-1];
		if(i>=2&&j>=2) x+=2*k*a[i-2][j-2]-b[i-2][j-2];
		x=(x+4*P)%P,b[i][j]=x,a[i][j]=inv[i]*x%P;
		if(i<=k) a[j][i]=a[i][j],b[j][i]=j*a[j][i]%P;
	    }
	for(int i=0;i<=k+k;++i) for(int j=0;j<=k;++j) inc(ans,(j==k? 1:2)*pow((P+i-k)*w+k-j,n)*a[i][j]%P);
	printf("%lld",ans*pow(d,P-k-1)%P);
    }
}
int main()
{
    scanf("%d%d%d",&n,&k,&d),init();
    if(n%d) puts("0"); else if(d==4) Sub1::solve(); else if(d==6) Sub2::solve();
}