[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';
    }
}