决战 dp+矩阵乘法+NTT

来源：互联网发布：网络综合布线设计图编辑：程序博客网时间：2024/04/29 19:39

题意

这里写图片描述
n<=2500

分析

有个比较显然的dp就是设f[i,j,s]表示dp到第i行，用了j个哲学家，以及第i行状态为s时的方案数。
如果直接暴力转移的话，复杂度就是O(64nm)，有点难以接受，考虑优化。
把状态设成f[i,s,j]，由于s的大小只有8，可以考虑用矩阵乘法来优化。我们可以把这状态看成是一个多项式fi,s(x)，那么第j项的系数就是用了j个哲学家的方案数。
转移的时候，由状态s1转移到s2，假设s2有c个1，那么一次转移就可以写成fi+1,s2(x)+=fi,s1(x)∗xc。
不难观察到，如果我们建出状态转移矩阵后A，Ai,j的值就为xcnt(j)。那么我们就可以通过矩阵快速幂来优化dp。
至于多项式乘法的话，我们可以先把多项式转化为点值表达式，然后两个多项式相乘就等于点值的直接相乘。

代码

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>using namespace std;typedef long long LL;const int N=2505;const int MOD=998244353;const int MAXM=8192;int a[4][4],hefa[8],trans[8][8],size[8]={0,1,1,2,1,2,2,3},t[8][MAXM],r[8][MAXM],bin[4],lg,M,rev[MAXM],n,m,e[N];struct Matrix{int a[8][8][MAXM];}g,f,tmp;void prework(){    for (int i=1;i<=3;i++)        for (int j=1;j<=3;j++)            scanf("%d",&a[i][j]);    a[1][2]=a[3][2]=a[1][2]|a[3][2];    a[2][1]=a[2][3]=a[2][1]|a[2][3];    a[1][1]=a[3][3]=a[1][1]|a[3][3];    a[1][3]=a[3][1]=a[1][3]|a[3][1];    for (int i=0;i<bin[3];i++)    {        hefa[i]=1;        for (int j=0;j<2;j++)            if ((i&bin[j])&&(i&bin[j+1])&&a[2][3]) hefa[i]=0;    }    for (int i=0;i<bin[3];i++)        for (int j=0;j<bin[3];j++)            if (hefa[j]&&hefa[i])            {                int flag=0;                for (int k=0;k<3;k++)                    for (int l=0;l<3;l++)                        if (abs(l-k)<=1&&(i&bin[k])&&(j&bin[l])&&a[3][2+l-k]) flag=1;                if (!flag) trans[i][j]=1;            }}int ksm(int x,int y){    int ans=1;    while (y)    {        if (y&1) ans=(LL)ans*x%MOD;        x=(LL)x*x%MOD;y>>=1;    }    return ans;}void NTT(int *a,int f){    for (int i=0;i<M;i++) if (i<rev[i]) swap(a[i],a[rev[i]]);    for (int i=1;i<M;i<<=1)    {        int wn=ksm(3,f==1?(MOD-1)/i/2:MOD-1-(MOD-1)/i/2);        for (int j=0;j<M;j+=(i<<1))        {            int w=1;            for (int k=0;k<i;k++)            {                int u=a[j+k],v=(LL)a[j+k+i]*w%MOD;                a[j+k]=(u+v)%MOD;a[j+k+i]=(u-v)%MOD;                w=(LL)w*wn%MOD;            }        }    }    int ny=ksm(M,MOD-2);    if (f==-1) for (int i=0;i<M;i++) a[i]=(LL)a[i]*ny%MOD;}void mul1(){    tmp=g;    for (int i=0;i<bin[3];i++)        for (int j=0;j<bin[3];j++)            for (int k=0;k<M;k++)                g.a[i][j][k]=0;    for (int i=0;i<bin[3];i++)        for (int k=0;k<bin[3];k++)            for (int j=0;j<bin[3];j++)                for (int l=0;l<M;l++)                    (g.a[i][j][l]+=(LL)tmp.a[i][k][l]*f.a[k][j][l]%MOD)%=MOD;}void mul2(){    tmp=f;    for (int i=0;i<bin[3];i++)        for (int j=0;j<bin[3];j++)            for (int k=0;k<M;k++)                f.a[i][j][k]=0;    for (int i=0;i<bin[3];i++)        for (int k=0;k<bin[3];k++)            for (int j=0;j<bin[3];j++)                for (int l=0;l<M;l++)                    (f.a[i][j][l]+=(LL)tmp.a[i][k][l]*tmp.a[k][j][l]%MOD)%=MOD;}void Matrix_ksm(int y){    g=f;y--;    while (y)    {        if (y&1) mul1();        mul2();y>>=1;    }}int main(){    bin[0]=1;    for (int i=1;i<=3;i++) bin[i]=bin[i-1]*2;    scanf("%d%d",&n,&m);    prework();    for (int i=0;i<bin[3];i++)        for (int j=0;j<bin[3];j++)            if (trans[i][j]) f.a[i][j][size[j]]=1;    for (M=1;M<=m||M<=8;M*=2,lg++);    for (int i=0;i<M;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));    for (int i=0;i<bin[3];i++)        for (int j=0;j<bin[3];j++)            NTT(f.a[i][j],1);    Matrix_ksm(n);    int ans=0;    r[0][0]=1;    NTT(r[0],1);    for (int i=0;i<bin[3];i++)    {        for (int j=0;j<bin[3];j++)            for (int k=0;k<M;k++)                (t[i][k]+=(LL)r[j][k]*g.a[j][i][k]%MOD)%=MOD;        NTT(t[i],-1);        (ans+=t[i][m])%=MOD;    }    printf("%d",(ans+MOD)%MOD);    return 0;}

阅读全文

0 0