[2018集训队做业][UOJ450] 复读机 [DP+泰勒展开+单位根反演]

题面

传送门

思路

本文中全部$m$是原题目中的$k$

首先，这个一看就是$d=1,2,3$数据分治

d=1

不说了，很简单，$m^n$

d=2

先上个$dp$试试

设$dp[i][j]$表示前$i$个复读机用掉了$j$个机会，注意这个东西最后求出来的是分配方案，还要乘以一个$n!$

$dp[i][j]=\sum_{k=0}^j [d|k]\binom{n-j+k}{k}dp[i-1][j-k]$

$dp[i][j]=\sum_{k=0}^j [d|k]\frac{(n-j+k)!}{(n-j)!k!}dp[i-1][j-k]$

$(n-j)!dp[i][j]=\sum_{k=0}^j [d|k]\frac{1}{k!}(n-j+k)!dp[i-1][j-k]$

咱们令生成函数$A(x)=\sum_{i=0}^{\infty}[d|i]\frac{x^i}{i!}$，$B_i(x)=\sum{j=0}^{\infty}(n-j)!dp[i][j]$

那么能够发现$B_{i+1}(x)=B_i(x)\ast A(x)$

也就是答案等于$A^m(x)$的第$n$项系数

咱们看这个$A(x)$的形式，发现它下面有一堆阶乘，不禁得让咱们联想到泰勒展开

（我也不知道这个是怎么联想的不过就这样吧我会再写一篇博客解释的23333）

咱们发现$e^x=\sum_{i=0}^{\infty}\frac{x^i}{i!}$，同时$e^{-x}=\sum_{i=0}^{\infty} (-1)^i \frac{x^i}{i!}$

那么易得$A(x)=\frac{e^x+e^{-x}}{2}$

因此$A^m(x)=(\frac{e^x+e^{-x}}{2})^m=\frac{1}{2^m}\sum_{i=0}^m \binom{m}{i}e^{(2i-m)^x}$

而后咱们考虑$n$次项系数，发现最外面应该最后乘上去的$n!$和里面的$\frac{1}{n!}$抵消了，上面$e$的幂剩下的系数是$2i-m$

这样咱们能够获得答案的表达式$ANS=\sum_{i=0}^m \binom{m}{i} (2i-m)$

d=3

emmm

咱们亲爱的$e$好像用不了了

可是咱们这个时候有一个神秘的东西：单位根反演！

单位根反演的公式是：$[d|i]=\frac{1}{d}\sum_{j=0}^{d-1}\omega_d^{ij}$

其中的$\omega_d^j$表示$d$阶单位根的$j$次方

咱们代入上面的公式里面获得：

$A(x)=\sum_{i=0}^{\infty}\sum_{j=0}^{d-1}\frac{1}{d}\omega_d^{ij}\frac{x^i}{i!}$

$A(x)=\frac{1}{d}\sum_{i=0}{d-1}e^{\omega_d^ix}$

其实能够看到上面的$d=2$就是这个式子的特殊状况

那么$d=3$怎么搞呢？

咱们能够发现模数$19491001$是一个3的倍数+1的形式，那么必然存在一个三阶单位负数根根$g$（考虑费马小定理便可）

咱们把这个原根求出来，而后两次暴力展开二项式定理，最后能够获得：

$A^m(x)=\frac{1}{3^m}\sum_{i=0}^m \sum_{j=0}^{m-i} \binom{m}{i}\binom{m-i}{j}e^{(i+gj+g^2(m-i-j))x}$

而后就$O(m^2\log m)$作完了

Code

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cassert>
#define MOD 19491001
#define ll long long
using namespace std;
inline ll read(){
    ll re=0,flag=1;char ch=getchar();
    while(!isdigit(ch)){
        if(ch=='-') flag=-1;
        ch=getchar();
    }
    while(isdigit(ch)) re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
    return re*flag;
}
inline void add(ll &a,ll b){
    a+=b;
    if(a>=MOD) a-=MOD;
}
inline ll qpow(ll a,ll b){
    ll re=1;
    while(b){
        if(b&1) re=re*a%MOD;
        a=a*a%MOD;b>>=1;
    }
    return re;
}
ll n,m,d,f[1000010],finv[1000010],inv[1000010],g=7,inv3,w1,w2;
void init(){
    ll i,len=1000000;
    f[0]=f[1]=finv[0]=finv[1]=inv[1]=1;
    for(i=2;i<=len;i++) f[i]=f[i-1]*i%MOD;
    finv[len]=qpow(f[len],MOD-2);
    for(i=len;i>2;i--) finv[i-1]=finv[i]*i%MOD;
    for(i=2;i<=len;i++) inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
}
inline ll C(ll x,ll y){
    return f[x]*finv[y]%MOD*finv[x-y]%MOD;
}
int main(){
    n=read();m=read();d=read();
    ll i,j;ll ans=0,tmp;init();
    inv3=qpow(3,MOD-2);
    w1=qpow(g,(MOD-1)/3);
    w2=w1*w1%MOD;
    if(d==1) cout<<qpow(m,n)<<'\n';
    if(d==2){
        for(i=0;i<=m;i++){
            add(ans,C(m,i)*qpow((2*i-m+MOD)%MOD,n)%MOD);
        }
        cout<<ans*qpow(qpow(2,m),MOD-2)%MOD<<'\n';
    }
    if(d==3){
        for(i=0;i<=m;i++){
            for(j=0;j<=m-i;j++){
                tmp=(i+w1*j+w2*(m-i-j))%MOD;
                add(ans,C(m,i)*C(m-i,j)%MOD*qpow(tmp,n)%MOD);
            }
        }
        cout<<ans*qpow(qpow(3,m),MOD-2)%MOD<<'\n';
    }
}