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

第一次自己切NTT,感觉NTT就是FFT的取模版本,具体可以看 http://blog.csdn.net/acdreamers/article/details/39026505 讲得很清楚

然后系数本来是可以 预处理出来,时间可以缩很多,但我为了模板的简洁还是在想求系数了。

说了那么多废话,下面我们进入正题:题目最简单的dp可以设做dp[i][j]:前i位积mod m 为j的方案数。由于乘法再取模不能直接转移,所以我们找到m的原根,把每个数换成m的原根的x次幂 这样就能直接转移了。我们发现这个转移的复杂度是n*m^2的,由于m^2的部分可以看做是个多项式的乘法,可以用NTT简化成mlogm 对于n 我们可以发现n次转移都是相同的 所以可以用快速幂求出。然后此题就顺利解决了。

#include <iostream>
#include <cstring>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
#define maxn 200000
#define mod 1004535809
#define ll long long
int quickpow(int a,int b,int c)
{
    int res=1;
    for(;b;b>>=1)
    {
        if(b&1) res=(1LL*res*a)%c;
        a=(1LL*a*a)%c;
    }
    return res;
}
int n,invn,nn,R[maxn],p,root,x,m;
void NTT(int *a,int f)
{
    int id=0;
    for(int i=1;i<nn;i++) if(i<R[i]) swap(a[i],a[R[i]]);
    for(int i=1;i<nn;i<<=1)
    {
        id++;
        int wn=quickpow(3,f==1 ? (mod-1)/(1<<id) : mod-1-(mod-1)/(1<<id),mod);
        for(int j=0;j<nn;j+=(i<<1))
        {
            int w=1;
            for(int k=0;k<i;k++)
            {
                int x=a[j+k],y=(1LL*w*a[j+k+i])%mod;
                a[j+k]=(1LL*(x+y))%mod;
                a[j+k+i]=((1LL*(x-y))%mod+mod)%mod;
                w=(1LL*w*wn)%mod;
            }
        }
    }
    if(f==-1)
    {
        for(int i=0;i<nn;i++) a[i]=(1LL*a[i]*invn)%mod;
    }
}
int T[maxn],mm;
int aa[maxn],bb[maxn],cc[maxn],dd[maxn],pos[maxn];
inline bool judge(int x, int p)
{
    for (int i = 2; i * i <= p; i ++)
        if ((p - 1)%i==0&&quickpow(x,(p-1)/i,p)==1) return 0;
    return 1;
}
inline int Find_Root(int p)
{
    if (p == 2) return 1;
    int res = 2;
    for (; !judge(res, p); res ++) ;
    return res;
}
int num[maxn];
void init()
{
    scanf("%d%d%d%d",&n,&p,&x,&m);
    for(int i=1;i<=m;i++)
        scanf("%d",&T[i]);
    mm=p*2-2;
    int l=0;
    for(nn=1;nn<=mm;nn<<=1)l++;
    for(int i=1;i<nn;i++)
        R[i]=(R[i>>1]>>1)|((i&1)<<(l-1));
    invn=quickpow(nn,mod-2,mod);
    root = Find_Root(p);
    for (int i=0;i<p-1;i++)
    {
        num[i]=!i ? 1 : num[i-1]*root%p;
        pos[num[i]] = i;
    }
}
void mul(int *ret,int *num1,int *num2)
{
    for (int i = 0; i < nn; i ++)
        cc[i] = num1[i],dd[i]=num2[i];
    NTT(cc,1);NTT(dd,1);
    for(int i=0;i<nn;i++)
    {
        ret[i]=(1LL*cc[i]*dd[i])%mod;
    }
    NTT(ret,-1);
    for (int i=nn-1;i>=p-1;i--)
    {
        ret[i-p+1]=1LL*(ret[i-p+1]+ret[i])%mod;
        ret[i] = 0;
    }
}
void quickpow2(int *a,int b)
{
    aa[0]=1;
    for(;b;b>>=1)
    {
        if(b&1)
        {
            mul(aa,a,aa);
        }
        mul(a,a,a);
    }
}
void solve()
{
    for (int i=1;i<=m;i++)
    {
        if(T[i] == 0) continue ;
        bb[pos[T[i]]]++;
    }
    quickpow2(bb,n);
    int ans=aa[pos[x]];
    printf("%d\n",ans);
}
int main()
{
    init();
    solve();
    return 0;
}


  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值