[JZOJ5364]史莱姆

题目大意

给你一个棱长为2n的正方体，你可以进行m次攻击，求攻击完后正方体样子的方案数。
一次攻击定义为：你对一个棱长为2i的正方体进行分解，分解成8个2i−1次方个小正方体。
两种方案不同定义为：存在一对棱长为1的正方体单位x,y，在一种方案中他们属于同一个正方体而另一种方案他们属于不同正方体。
多组数据。
T≤100000,n≤40，m≤30000。

分析

一个小DP：
设f[i][j][k]表示考虑到大小为2n−i的大小的正方体，比他大的正方体一共用了j次攻击，当前大小的正方体有8*k个，的方案数。
很显然f[i][j][k]=∑f[i+1][j−k][l]∗c(k,l∗8)
这就是部分分了。
这东西没办法优化了，我们要换个思路。
看到统计方案的题，想想生成函数总是没问题的，首先要推出一个能够搞生成函数的dp。
我们反过来，合成一个2n的正方体。
设f[i][j]表示一个2i的正方体，用j次攻击的方案数，转移式是：
f[i][j]=∑a[1..8]|sum(a)=j−1(f[i−1][a[1]]∗⋅⋅⋅f[i−1][a[8]])，即枚举攻击一次后分解成的8个小正方体的分割情况。特别的，f[i][0]=1。
现在我们把f(i)看成一个多项式，f[i][j]表示f[i]的第j项的系数，就可以FFT了！
f(i)=f(i−1)8∗x+1，那么只要对f(i-1)求幂，然后整体右移一位，再令常数项+1，就是f(i)了。把答案记录下来直接输出即可。
这个NTT有大半年没打···低级错误犯了好多，回顾一下。
len为原多项式次数界之和,再调整到2^x。
递归过程中：对n次多项式带入w1n,w2n...wnn
我们DFT的过程中，求
t[i]=a[0],a[1]∗win,a[2]∗w2in...
t[i+n/2]=a[0],a[1]∗wi+n/2n,a[2]∗wi+nn...
i∈[0,len−1]
对t[i]分析，分组转化得
(a[0],a[2]∗win/2,a[4]∗w2in/2)
win∗(a[1],a[3]∗win/2,a[5]∗w2in/2)
括号里的东西分别储存在t′[i],t′[i+n/2]中。
t[i]=t′[i]+t′[i+n/2]∗win
t[i+n/2]=t′[i]+t′[i+n/2]∗wi+n/2n=t′[i]+t′[i+n/2]∗w−in
二进制表示的下标中，a为dft目标数组，t为过程数组，t[0101]=a[1010]
dft后，两个数组相乘，所有的点值都要乘，不然插值会出错。
插值就把wlen−i带入，最后出来的多项式系数要除以len。
预处理原根下的那些主次方根。

代码

#include<cstdio>#include<algorithm>#include<cmath>#include<cstring>#include<set>using namespace std;#define fo(i,j,k) for(i=j;i<=k;i++)#define fd(i,j,k) for(i=j;i>=k;i--)typedef long long ll;typedef double db;const int N=1e5+5,mo=998244353,rt=3;int Log,w[N],len,n,m,T,i,f[42][N],g[N],a[N],invl,j,k;int t[N];struct rec{    int n,m,id;}b[N];int ksm(int x,int y){    int ret=1;    while (y)    {        if (y&1) ret=1ll*ret*x%mo;        y>>=1;        x=1ll*x*x%mo;    }    return ret;}int mx2(int x){    int ret=1;Log=0;    while (ret<x) ret*=2,Log++;    return ret;}void predo(){    int i;    w[0]=1;    w[1]=ksm(rt,(mo-1)/len);    fo(i,2,len) w[i]=1ll*w[i-1]*w[1]%mo;}void DFT(int *a,int n,int sig){    int i,j,ws,half,siz,k,pos,u,v;    fo(i,0,n-1)    {        int tp=0;        fo(j,0,Log-1) tp=(tp<<1)+(((1<<j)&i)!=0);        t[tp]=a[i];    }    for(ws=1,siz=2;siz<=n;siz*=2,ws++)    {        half=siz>>1;        fo(i,0,half-1)        {            for(j=i;j<n;j+=siz)            {                k=j+half;                pos=len/siz*i*sig;                if (pos<0) pos+=len;                u=t[j];v=1ll*t[k]*w[pos]%mo;                t[j]=(u+v)%mo;                t[k]=(u-v+mo)%mo;            }        }    }    fo(i,0,n-1) a[i]=t[i];}int main(){    scanf("%d\n",&T);    fo(i,1,T)    {        scanf("%d %d",&b[i].n,&b[i].m);        n=max(n,b[i].n);        m=max(m,b[i].m);        b[i].id=i;    }    f[0][0]=1;    len=mx2(m*2+2);    invl=ksm(len,mo-2);    predo();    fo(i,1,n)    {        fo(j,0,m) g[j]=f[i-1][j];        fo(j,1,3)        {            fo(k,0,m) a[k]=g[k];            DFT(a,len,1);            fo(k,0,len-1) a[k]=1ll*a[k]*a[k]%mo;            DFT(a,len,-1);            fo(k,0,m) g[k]=1ll*a[k]*invl%mo;            fo(k,0,len-1) a[k]=0;        }        fo(j,1,m) f[i][j]=g[j-1];        f[i][0]=1;    }    fo(i,1,T) printf("%d\n",f[b[i].n][b[i].m]);}