icontofig | 发布于 2019-09-26 19:59:05 | 阅读量 556 | 生成函数 NTT
发布于 2019-09-26 19:59:05 | 生成函数 NTT

题解

如果你曾经经过大量的FFT/NTT的生成函数训练,应该一眼能够看出这道题一定是个多项式卷积用NTT优化……

如果是相加的话,我们可以直接设多项式S = Σaixi,(生成函数在这)其中ai表示数字i出现过多少次。

如果我们是选2个数的话,对于每一个和的方案数,很明显就是将S卷上S,然后取对应多项式项的系数即为方案数。

因为这个题是乘积,而且对于乘积还有对质数m取模操作,所以我们不能简单的按照上面的方法进行生成函数构造。

如何解决乘积问题呢?

我们知道多项式FFT/NTT都是建立在指数(下标)相加的基础上(也就是普通的多项式乘法的法则)进行的,所以我们想能不能把乘积变成加法,这样就可以按照经典套路进行生成函数构造了。

想了一下,取对数貌似可以把乘法变成加法来算了……

但是因为此处还是对质数m取模之后结果的统计,而且NTT也不支持浮点数,所以我们需要考虑在模意义下的对数,也就是离散模对数。

如何处理离散模对数呢……这里如果不了解的话就需要补一下基础的数论知识了。

这里我们需要求m的原根,求法也非常简单,直接暴力枚举的每个数,看对于每个质数pi | (m-1),是否g的(m-1)/pi次方mod m等于1,如果等于则说明不是原根,否则说明是原根,原根通常非常小,所以不用担心会在处理原根的时候TLE。

运用原根,我们可以把每个数在mod m意义下的离散模对数求出来,然后就变成正常的多项式卷积啦!

于是这个时候我们就可以用NTT来做啦。

但是由于要选N次,而N的数据范围到了1e9次,所以我们不能直接暴力NTT,要用快速幂的思想对多项式进行快速幂形式的卷积操作。

这里要注意每一次卷积完成的时候(IDNT之后),我们只是正常的进行了一次多项式卷积,统计答案的时候还是有问题的,因为指数并不能取模,所以每一卷积,都会有答案落在[m-1,2*m-2]上,我们要把这一区域上的答案对应地加到前面[0,m-2]上(对应了题目中的乘积取模操作)。

总时间复杂度为O(nlog2n)。

代码

#include <bits/stdc++.h>
using namespace std;
const int maxm = 16005;
typedef long long ll;
const ll gf = 3;
int n,m,x,s;
ll tmp;
int sq;
const ll mod = 1004535809;
ll fpow(ll a,ll b,ll mo){
    ll ret = 1;
    while(b){
        if(b & 1)ret = ret * a % mo;
        a = a * a % mo;
        b >>= 1;
    }
    return ret;
}
int rev[maxm<<2];
int len,lg;
ll idx[maxm<<2],b[maxm<<2],f[maxm<<2];
void init(int q,int k){
    for(len = 1;len < q + k;len <<= 1,++lg);
    for(int i = 0;i < len;++i)
        rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (lg - 1));
    return;
}
void NTT(ll *a,int len,int f){
    for(int i = 0;i < len;++i)
        if(i < rev[i])
            swap(a[i],a[rev[i]]);
    for(int i = 2;i <= len;i <<= 1){
        int now = i >> 1;
        ll wn = fpow(gf,(mod-1)/i,mod);
        if(f == -1)wn = fpow(wn,mod-2,mod);
        for(int j = 0;j < len;j += i){
            ll wk = 1,x,y;
            for(int k = j;k < j + now;k++,wk = wk * wn % mod){
                x = a[k];y = a[k + now] * wk % mod;
                a[k] = (x + y) % mod;
                a[k + now] = (x - y + mod) % mod;
            }
        }
    }
    if(f == -1){
        ll inv = fpow(len,mod-2,mod);
        for(int i = 0;i < len;++i)
            a[i] = a[i] * inv % mod;
    }
    return;
}
vector<ll>qk;
ll find_root(ll x){
    ll tmp = x - 1;
    for(ll i = 2;i * i <= x;++i){
        if(tmp % i == 0){
            qk.push_back(i);
            while(tmp % i == 0)tmp /= i;
        }
    }
    if(tmp != 1)qk.push_back(tmp);
    for(ll g = 2;g < x;++g){
        bool flag = true;
        for(int j = 0;j < qk.size();++j){
            if(fpow(g,(x-1)/qk[j],x) == 1){
                flag = false;
                break;
            }
        }
        if(flag)return g;
    }
}
void work(const ll *a1,const ll *a2,int len,ll *ans){
    static ll b1[maxm<<2],b2[maxm<<2];
    fill(b1,b1+len,0);
    fill(b2,b2+len,0);
    copy(a1,a1+len,b1);
    copy(a2,a2+len,b2);
    NTT(b1,len,1);
    NTT(b2,len,1);
    for(int i = 0;i < len;++i)
        b1[i] = b1[i] * b2[i] % mod;
    NTT(b1,len,-1);
    for(int i = m-1;i < 2*m-1;++i)
        b1[i-(m-1)] = (b1[i-(m-1)] + b1[i]) % mod;
    for(int i = 0;i < m-1;++i)
        ans[i] = b1[i];
    return;
}
int main(){
    ios_base::sync_with_stdio(0);
    cin.tie(0);
    cin >> n >> m >> x >> s;
    ll g = 0;
    g = find_root(m);
    tmp = 1;
    init(m-1,m-1);
    for(int i = 0;i < m-1;++i){
        idx[tmp] = i;
        tmp = tmp * g % m;
    }
    fill(f,f+len,0);
    for(int i = 1;i <= s;++i){
        cin >> sq;
        if(sq == 0)continue;
        f[idx[sq]] += 1;
    }
    fill(b,b+len,0);
    b[0] = 1;
    while(n){
        if(n & 1)
            work(b,f,len,b);
        work(f,f,len,f);
        n >>= 1;
    }
    cout << b[idx[x]] % mod << endl;
    return 0;
}



内容更新于: 2019-09-26 19:59:05
链接地址: http://blog.leanote.com/post/icontofig/%5BSDOI2015%5D%E5%BA%8F%E5%88%97%E7%BB%9F%E8%AE%A1-NTT-%E7%94%9F%E6%88%90%E5%87%BD%E6%95%B0

上一篇: 牛客国庆欢乐赛5 2017四川省赛 K-2017Revenge bitset加速dp、原根

下一篇: BZOJ 3160 万径人踪灭 FFT + 回文自动机 + 生成函数

556 人读过
立即登录, 发表评论.
没有帐号? 立即注册
0 条评论
文档导航