算法一
定义 (f[i][j]) 表示前 (i) 个珍珠，个数为奇数的种类有 (j) 种，这样状态是 (n*D) 的，转移是 (O(1)) 的，直接讨论当前选的的颜色是在奇数的 (j) 种中还是偶数的 ((D-j)) 种中
答案即为 (sum_{i=0}^{n-m*2}f[n][i])
复杂度 (O(n*D)) ,期望得分 (52)
算法二
(D leq 300，n,m leq 10^9)
按照算法一的定义，矩阵快速幂+极致卡常
这一部分复杂度是 (O(D^3*logn)) ，结合算法一期望得分 (60)
算法三
特判 (m=0) 的，(ans=D^n) ， (m=1) 的， (ans=D^n-A_D^n)
这一部分复杂度 (O(1)) ，结合算法一，二，三，期望得分 (64)
算法四
二项式反演+生成函数
(n-m*2) 即为合法情况下，最多有几种宝石的个数为奇数
首先特判掉 (n-m*2<0，ans=0) 和 (n-m*2geq D，ans=D^n)
定义 (f_i) 表示 (D) 种宝石，恰好 (i) 种个数为奇数的方案数
(g_i) 表示 (D) 种宝石，钦定 (i) 种个数为奇数，其余的任意选的方案数（至少 (i) 种为奇数）
对于一种颜色， (cnt) 为奇数的方案数对应的 (EGF) 为 (frac{e^x-e^{-x}}{2})， (cnt) 任意的方案数对应的 (EGF) 为 (e^x)
则
二项式定理展开 ((e^x-e^{-x})^i)
提出 (-1) 然后合并 (e) 的幂
将 ((e^x)^{D-2(i-j)}) （系数为 (frac{[D-2(i-j)]^n}{n!}) ）和组合数展开
换个枚举量，枚举 (i-j)，
(sum) 后面就是卷积的形式，直接 (NTT) 算就好了
得到了 (g)，然后需要二项式反演求出 (f)，但是不能直接 (2) 层 (for) 循环去求，那样是 (O(d^2)) 的
注意反演的式子
后面又是一个卷积，但是是差一定，所以枚举 (j) 从 (0) 到 (D-i) ， (NTT) 求出 (Ans_{D-i}=sum_{j=0}^{D-i}A_jB_{D-i-j})
其中 (A_i=frac{(-1)^i}{i!})， (B_{D-i}=i!g_i)
那么 (f_i=frac{1}{i!}Ans_{D-i})
复杂度 (O(DlogD)) ，期望得分 (100)
代码
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define ll long long
using namespace std;
const int maxn=1<<18|10,mod=998244353,Gi=332748118,G=3;
int n,m,D;
int rev[maxn],a[maxn],b[maxn],f[maxn],g[maxn],fac[maxn],fy[maxn];
int fp(int a,int x,int ans=1){
for(;x;x>>=1,a=1ll*a*a%mod) if(x&1) ans=1ll*ans*a%mod;
return ans;
}
void NTT(int *a,int lim,int opt){
for(int i=0;i<lim;++i) if(i<rev[i]) swap(a[i],a[rev[i]]);
for(int d=1;d<lim;d<<=1){
int st=d<<1,w1=fp(opt==-1?Gi:G,(mod-1)/st);
for(int i=0;i<lim;i+=st){
int w=1;
for(int j=0;j<d;++j,w=1ll*w*w1%mod){
int x=a[i+j],y=1ll*w*a[i+j+d]%mod;
a[i+j]=(x+y)%mod;
a[i+j+d]=(x+mod-y)%mod;
}
}
}
if(opt!=-1) return;
for(int i=0,inv=fp(lim,mod-2);i<lim;++i) a[i]=1ll*a[i]*inv%mod;
}
void Init(int n){
fac[0]=fy[0]=1;
for(int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%mod;
fy[n]=fp(fac[n],mod-2);
for(int i=n;i>1;--i) fy[i-1]=1ll*fy[i]*i%mod;
}
int main(){
scanf("%d%d%d",&D,&n,&m);
Init(D);
if(n-2*m<0) return puts("0"),0;
if(n-2*m>=D) return printf("%d
",fp(D,n)),0;
for(int i=0;i<=D;++i){
if(i&1) a[i]=1ll*(mod-fp((D-2*i+mod)%mod,n))*fy[i]%mod;
else a[i]=1ll*fp((D-2*i+mod)%mod,n)*fy[i]%mod;
b[i]=fy[i];
}
int lim=1,bit=0;
while(lim<=D*2) lim<<=1,++bit; --bit;
for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<bit);
NTT(a,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;++i) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,lim,-1);
for(int i=0;i<=D;++i) g[i]=1ll*fac[D]*fy[D-i]%mod*fp(fp(2,i),mod-2)%mod*a[i]%mod;
for(int i=0;i<=D;++i){
if(i&1) a[i]=mod-fy[i];
else a[i]=fy[i];
b[D-i]=1ll*fac[i]*g[i]%mod;
}
for(int i=D+1;i<lim;++i) a[i]=b[i]=0;
NTT(a,lim,1),NTT(b,lim,1);
for(int i=0;i<lim;++i) a[i]=1ll*a[i]*b[i]%mod;
NTT(a,lim,-1);
for(int i=0;i<=D;++i) f[i]=1ll*a[D-i]*fy[i]%mod;
int ans=0;
for(int i=0;i<=n-2*m;++i) (ans+=f[i])%=mod;
printf("%d
",(ans+mod)%mod);
return 0;
}