BZOJ3992序列统计(动态规划,NTT)

Posted 小蒟蒻yyb的博客

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了BZOJ3992序列统计(动态规划,NTT)相关的知识,希望对你有一定的参考价值。

【BZOJ3992】序列统计(动态规划,NTT)

题面

BZOJ

题解

最裸的暴力
\(f[i][j]\)表示前\(i\)个数,积在膜意义下是\(j\)的方案数
转移的话,每次枚举一个数,直接丢进去就好
复杂度\(O(nm|S|)\),10pts

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define RG register
#define MOD 1004535809
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
int n,m,X,T;
int f[2][10000];
int a[10000];
int main()
{
    n=read();m=read();X=read();T=read();
    for(int i=1;i<=T;++i)a[i]=read()%m;
    for(int i=1;i<=T;++i)f[1][a[i]]++;
    for(int i=2;i<=n;++i)
    {
        for(int j=0;j<m;++j)f[i&1][j]=0;
        for(int j=1;j<=T;++j)
            for(int k=0;k<m;++k)
                (f[i&1][k*a[j]%m]+=f[(i+1)&1][k])%=MOD;
    }
    printf("%d\n",f[n&1][X]);
    return 0;
}

发现每一步的转移是相同的,
因此可以矩阵快速幂
时间复杂度\(O(lognm^3)\),30pts
我懒得写了


我们都发现了转移是相同的
那么不一定只能用矩阵快速幂呀
我们的转移也是满足结合律的
所以可以把转移跑快速幂
复杂度\(O(lognm^2)\),60pts

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define RG register
#define MOD 1004535809
#define MAX 10000
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
int n,m,X,T;
int f[MAX],s[MAX];
int a[MAX],ret[MAX];
int* zy(int *a,int *b)
{
    memset(ret,0,sizeof(ret));
    for(int i=0;i<m;++i)
        for(int j=0;j<m;++j)
            (ret[i*j%m]+=1ll*a[i]*b[j]%MOD)%=MOD;
    return ret;
}
int main()
{
    n=read();m=read();X=read();T=read();
    for(int i=1;i<=T;++i)a[i]=read()%m;
    for(int i=1;i<=T;++i)f[a[i]]++;
    bool fl=false;
    int b=n;
    while(b)
    {
        if(b&1)
        {
            if(fl)
            {
                memset(ret,0,sizeof(ret));
                for(int i=0;i<m;++i)
                    if(s[i])
                        for(int j=0;j<m;++j)
                            (ret[i*j%m]+=(1ll*s[i]*f[j])%MOD)%=MOD;
                for(int i=0;i<m;++i)s[i]=ret[i];
            }
            else
            {
                for(int i=0;i<m;++i)s[i]=f[i];
                fl=true;
            }
        }
        memset(ret,0,sizeof(ret));
        for(int i=0;i<m;++i)
            if(f[i])
                for(int j=0;j<m;++j)
                    (ret[i*j%m]+=(1ll*f[i]*f[j])%MOD)%=MOD;
        for(int i=0;i<m;++i)f[i]=ret[i];
        b>>=1;
    }
    printf("%d\n",s[X]);
    return 0;
}

现在就是最大的问题了
\(n\)已经优化到了\(logn\)
转移现在才是最大的问题
我们发现转移是这样的:
\(f[i]*f[j]\to f[i*j]\)
如果它长成这个样子:
\(f[i]*f[j]\to f[i+j]\)
这样子的话就会做啦
这样就可以跑一遍多项式的卷积

怎么转换呢?
题目给定的条件\(m\)是质数
我们知道\(x^{\varphi(m)}\% m=1\)
而如果\(x\)\(m\)的原根
那么,对于\(0~\varphi(m)\)
每一个原根的若干次幂恰好对应一个数
那么,这样的话,乘法可以转换成幂的加法

于是,直接跑多项式的卷积就好了
因为要取膜,只能跑\(NTT\)

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<set>
#include<map>
#include<vector>
#include<queue>
using namespace std;
#define ll long long
#define RG register
#define MOD 1004535809
#define MAX 100000
inline int read()
{
    RG int x=0,t=1;RG char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}
const int pr=3;
const int phi=MOD-1;
int n,m,X,T;
int f[MAX],s[MAX];
int a[MAX],mp[MAX],b[MAX];
int N,M,l,r[MAX],ret[MAX];
int fpow(int a,int b,int P)
{
    int s=1;
    while(b){if(b&1)s=1ll*s*a%P;a=1ll*a*a%P;b>>=1;}
    return s;
}
int ys[MAX],yst;
int getroot(int n)
{
    int tmp=n-1;
    for(int i=2;i*i<=tmp;++i)
        if(tmp%i==0)
        {
            ys[++yst]=i;
            while(tmp%i==0)tmp/=i;
        }
    if(tmp>1)ys[++yst]=tmp;
    for(int g=2;g<=n-1;++g)
    {
        bool fl=true;
        for(int i=1;i<=yst;++i)
            if(fpow(g,(n-1)/ys[i],n)==1){fl=false;break;}
        if(fl)return g;
    }
    return -1;
}
void getmap()
{
    int prm=getroot(m);
    for(int i=0;i<m-1;++i)
        mp[fpow(prm,i,m)]=i;
}
void preNTT()
{
    M=2*(m-2);
    for(N=1;N<=M;N<<=1)++l;
    for(int i=0;i<N;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}
void NTT(int *P,int opt)
{
    for(int i=0;i<N;++i)if(i<r[i])swap(P[i],P[r[i]]);
    for(int i=1;i<N;i<<=1)
    {
        int W=fpow(pr,phi/(i<<1),MOD);
        for(int p=i<<1,j=0;j<N;j+=p)
        {
            int w=1;
            for(int k=0;k<i;++k,w=1ll*w*W%MOD)
            {
                int X=P[j+k],Y=1ll*P[i+j+k]*w%MOD;
                P[j+k]=(X+Y)%MOD;P[i+j+k]=(X-Y+MOD)%MOD;
            }
        }
    }
    if(opt==-1)
    {
        reverse(&P[1],&P[N]);
        int inv=fpow(N,MOD-2,MOD);
        for(int i=0;i<N;++i)P[i]=1ll*P[i]*inv%MOD;
    }
}
void zy(int *a1,int *a2,int *c)
{
    memset(a,0,sizeof(a));memset(b,0,sizeof(b));
    for(int i=0;i<m-1;++i)a[i]=a1[i],b[i]=a2[i];
    NTT(a,1);NTT(b,1);
    for(int i=0;i<N;++i)a[i]=1ll*a[i]*b[i]%MOD;
    NTT(a,-1);
    memset(ret,0,sizeof(ret));
    for(int i=0;i<m-1;++i)ret[i]=(a[i]+a[i+m-1])%MOD;
    for(int i=0;i<m-1;++i)c[i]=ret[i];
}
int main()
{
    n=read();m=read();X=read();T=read();
    getmap();preNTT();
    for(int i=1;i<=T;++i)
    {
        int x=read()%m;
        if(x)f[mp[x]]++;
    }
    s[mp[1]]=1;
    while(n)
    {
        if(n&1)zy(s,f,s);
        zy(f,f,f);
        n>>=1;
    }
    printf("%d\n",s[mp[X]]);
    return 0;
}

以上是关于BZOJ3992序列统计(动态规划,NTT)的主要内容,如果未能解决你的问题,请参考以下文章

BZOJ 3992: [SDOI2015]序列统计 NTT+快速幂

[BZOJ3992][SDOI2015]序列统计(DP+原根+NTT)

[bzoj3992][SDOI2015]序列统计——离散对数+NTT

BZOJ3992[SDOI2015]序列统计 NTT+多项式快速幂

BZOJ[3992][SDOI2015]序列统计 生成函数+NTT

BZOJ 3992: [SDOI2015]序列统计 快速幂+NTT(离散对数下)