BZOJ3992(NTT+DP+快速幂）

最新推荐文章于 2019-03-14 10:02:00 发布

Fat_tu

最新推荐文章于 2019-03-14 10:02:00 发布

阅读量429

点赞数

CC 4.0 BY-SA版权

分类专栏：快速数论变换倍增 dp 文章标签： ntt dp 快速幂

本文链接：https://blog.youkuaiyun.com/q582116859/article/details/77287562

dp 同时被 3 个专栏收录

24 篇文章

订阅专栏

快速数论变换

15 篇文章

订阅专栏

倍增

8 篇文章

订阅专栏

该博客探讨了一道数学与算法结合的题目，涉及质数M、集合元素构成的序列以及模运算。通过设立DP状态，将乘积问题转换为加法的卷积形式，利用NTT（快速傅里叶变换）进行状态转移优化，并结合快速幂处理大数计算，降低复杂度至Mlog(M)log(N)。此外，还提及了寻找质数原根的方法。

题面
题意是给你一个质数M,和一个元素都小于M 的集合，大小为S，用集合中的元素构成长度为N的序列（可以用重复的元素），使其乘积模M为x，问方案数,答案模479*2^21+1。
M<=8000,N<=1e9

这题一看就很套路，模数为费马素数，大概就和NTT有关吧。
先考虑简单的DP，设f[i][j]为长度为i的序列，乘积模M等于j的方案数
枚举k，有

f [i + 1] [j * k 模 M] + = f [i] [j]

$f[i+1][j*k模M]+=f[i][j]$

看到M为质数，就必然存在原根，且集合中的元素都小于M，均可化成M的某次幂，且不重复，x同理。就把原问题转化成了一个长为N的序列，和为x的方案数，把乘法转换成了加法。有

f [i] [j] = \sum k = 0 j f [i - 1] [k] * f [i - 1] [j - k] + \sum k = 0 j + M - 2 f [i - 1] [k] * f [i - 1] [j + M - 2 - k]

$f[i][j]=\sum_{k=0}^{j}f[i-1][k]*f[i-1][j-k]+ \sum_{k=0}^{j+M-2}f[i-1][k]*f[i-1][j+M-2-k]$

是一个卷积形式，就可以用NTT优化状态转移。
看到N<=1e9，每次转移过程都类似，可以用快速幂优化。相当于枚举前半段的和与后半段的和，NTT进行状态转移。

总复杂度Mlog(M)log(N)。

对于求质数P原根，我也不太懂，大概这样吧
听说原根一般不大，从2开始枚举原根x，枚举P-1的约数i，i!=p-1。若x^i%P=1，则x不是原根。

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const LL mo=1004535809,N=100100;

int s,g,len,n,m,x,y;
int rev[N],gg=3;
LL ans[N],b[N],bb[N],little_busters[N];

LL cheng(LL a,LL b,LL p)
{
    LL ans=1ll;
    for(;b;b>>=1,a=a*a%p)
    if(b&1)
    ans=ans*a%p;
    return ans;
}

bool ok(int x,int p)
{
    for(int i=2;i*i<=p;i++)
    if((p-1)%i==0&&cheng(x, (p-1)/i,p)==1)
    return 0;
    return 1;
}

int find(int p)
{
    if(p==2)
    return 1;
    int res=2;
    for(;!ok(res,p);)
    res++;
    return res;
}

void init(int lim)
{
    n=1;
    int k=-1;
    while(n<lim)
    n<<=1,k++;
    for(int i=0;i<n;i++)
    rev[i]=(rev[i>>1] >> 1) | ((i&1)<<k);
}

void ntt(LL *a,bool ops)
{
    for(int i=0;i<n;i++)
    if(i<rev[i])
    swap(a[i],a[rev[i]]);

    for(int l=2;l<=n;l<<=1)
    {
        int m=l>>1;
        int wn;
        if(ops)
        wn=cheng(gg,(mo-1)/l,mo);
        else
        wn=cheng(gg,mo-1-(mo-1)/l,mo);
        for(int i=0;i<n;i+=l)
        {
            LL w=1;
            for(int k=0;k<m;k++)
            {
                LL t=w*a[i+k+m]%mo;
                a[i+k+m]=(a[i+k]-t+mo)%mo;
                a[i+k]=(a[i+k]+t)%mo;
                w=w*wn%mo;
            }
        }
    }
    if(!ops)
    {
        LL Inv=cheng(n,mo-2,mo);
        for(int i=0;i<n;i++)
        a[i]=a[i]*Inv%mo;
    }
}

void kiss(LL *a,LL *c)
{
    ntt(a,1);
    if(a!=c)
    {
        for(int i=0;i<n;i++)
        little_busters[i]=c[i];
        ntt(c,1);
    }
    for(int i=0;i<n;i++)
    a[i]=a[i]*c[i]%mo;
    ntt(a,0);
    for(int i=m-1;i<n;i++)
    {
        a[i-m+1]=(a[i-m+1]+a[i])%mo;
        a[i]=0;
    }

    if(a!=c)
    for(int i=0;i<n;i++)
    c[i]=little_busters[i];
}

int main()
{
    cin>>len>>m>>x>>s;

    g=find(m);

    for(int i=1;i<=s;i++)
    {
        int u;
        scanf("%d",&u);
        bb[u]=1;
    }

    y=0;
    for(;cheng(g,y,m)!=x;)
    y++;

    for(int i=0;i<m-1;i++)
    if(bb[cheng(g,i,m)])
    b[i]=1;

    init(m+m);

    ans[0]=1;

    for(;len;len>>=1)
    {
        if(len&1)
        kiss(ans,b);
        kiss(b,b);
    }
    cout<<ans[y]<<endl;

    return 0;
}