任意模数fft

1.正常拆系数fft,8次dft

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize(4)
//#pragma GCC optimize("unroll-loops")
//#pragma comment(linker, "/stack:200000000")
//#pragma GCC optimize("Ofast,no-stack-protector")
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")
#include<bits/stdc++.h>
//#include <bits/extc++.h>
#define fi first
#define se second
#define db double
#define mp make_pair
#define pb push_back
#define mt make_tuple
//#define pi acos(-1.0)
#define ll long long
#define vi vector<int>
#define mod 1000000007
#define ld long double
//#define C 0.5772156649
#define ls l,m,rt<<1
#define rs m+1,r,rt<<1|1
#define sqr(x) ((x)*(x))
#define pll pair<ll,ll>
#define pil pair<int,ll>
#define pli pair<ll,int>
#define pii pair<int,int>
#define ull unsigned long long
#define bpc __builtin_popcount
#define base 1000000000000000000ll
#define fin freopen("a.txt","r",stdin)
#define fout freopen("a.txt","w",stdout)
#define fio ios::sync_with_stdio(false);cin.tie(0)
#define mr mt19937 rng(chrono::steady_clock::now().time_since_epoch().count())
inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline void sub(ll &a,ll b){a-=b;if(a<0)a+=mod;}
inline void add(ll &a,ll b){a+=b;if(a>=mod)a-=mod;}
template<typename T>inline T const& MAX(T const &a,T const &b){return a>b?a:b;}
template<typename T>inline T const& MIN(T const &a,T const &b){return a<b?a:b;}
inline ll qp(ll a,ll b){ll ans=1;while(b){if(b&1)ans=ans*a%mod;a=a*a%mod,b>>=1;}return ans;}
inline ll qp(ll a,ll b,ll c){ll ans=1;while(b){if(b&1)ans=ans*a%c;a=a*a%c,b>>=1;}return ans;}

using namespace std;
//using namespace __gnu_pbds;

const ld pi=acos(-1);
const ull ba=233;
const db eps=1e-5;
const ll INF=0x3f3f3f3f3f3f3f3f;
const int N=400000+10,maxn=2000000+10,inf=0x3f3f3f3f;

struct cd{
    ld x,y;
    cd(ld _x=0.0,ld _y=0.0):x(_x),y(_y){}
    cd operator +(const cd &b)const{
        return cd(x+b.x,y+b.y);
    }
    cd operator -(const cd &b)const{
        return cd(x-b.x,y-b.y);
    }
    cd operator *(const cd &b)const{
        return cd(x*b.x - y*b.y,x*b.y + y*b.x);
    }
    cd operator /(const db &b)const{
        return cd(x/b,y/b);
    }
}a[N],b[N],c[N],d[N];
int rev[N],A[N],B[N],C[N];
void getrev(int bit)
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1) | ((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int step=1;step<n;step<<=1)
    {
        cd wn(cos(dft*pi/step),sin(dft*pi/step));
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);
            for(int k=j;k<j+step;k++)
            {
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;a[k+step]=x-y;
                wnk=wnk*wn;
            }
        }
    }
    if(dft==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;
}
void mtt(int n,int m,int p) {
    int sz=0;
    while((1<<sz)<=n+m)sz++;getrev(sz);
    int len=1<<sz;
    for(int i=0;i<len;i++)
    {
        int t1=A[i]%p,t2=B[i]%p;
        a[i]=cd(t1>>15,0),b[i]=cd(t1&0x7fff,0);
        c[i]=cd(t2>>15,0),d[i]=cd(t2&0x7fff,0);
    }
    fft(a,len,1);fft(b,len,1);fft(c,len,1);fft(d,len,1);
    for(int i=0;i<len;i++)
    {
        cd aa=a[i],bb=b[i],cc=c[i],dd=d[i];
        a[i]=aa*cc;b[i]=bb*cc;c[i]=aa*dd;d[i]=bb*dd;
    }
    fft(a,len,-1);fft(b,len,-1);fft(c,len,-1);fft(d,len,-1);
    for(int i=0;i<len;i++)
    {
        ll aa=(ll)(a[i].x+0.5),bb=(ll)(b[i].x+0.5);
        ll cc=(ll)(c[i].x+0.5),dd=(ll)(d[i].x+0.5);
        aa=(aa%p+p)%p;bb=(bb%p+p)%p;cc=(cc%p+p)%p;dd=(dd%p+p)%p;
        C[i]=((((aa<<15)%p)<<15)%p+(bb<<15)%p+(cc<<15)%p+dd)%p;
    }
}
int main()
{
    int n,m,p;scanf("%d%d%d",&n,&m,&p);
    for(int i=0;i<=n;i++)scanf("%d",&A[i]);
    for(int i=0;i<=m;i++)scanf("%d",&B[i]);
    mtt(n,m,p);
    for(int i=0;i<=n+m;i++)printf("%d ",C[i]);
    return 0;
}
/********************

********************/

2.黑科技优化拆系数fft,4次dft

//#pragma GCC optimize(2)
//#pragma GCC optimize(3)
//#pragma GCC optimize(4)
//#pragma GCC optimize("unroll-loops")
//#pragma comment(linker, "/stack:200000000")
//#pragma GCC optimize("Ofast,no-stack-protector")
//#pragma GCC target("sse,sse2,sse3,ssse3,sse4,popcnt,abm,mmx,avx,tune=native")
#include<bits/stdc++.h>
//#include <bits/extc++.h>
#define fi first
#define se second
#define db double
#define mp make_pair
#define pb push_back
#define mt make_tuple
//#define pi acos(-1.0)
#define ll long long
#define vi vector<int>
#define mod 1000000007
#define ld long double
//#define C 0.5772156649
#define ls l,m,rt<<1
#define rs m+1,r,rt<<1|1
#define sqr(x) ((x)*(x))
#define pll pair<ll,ll>
#define pil pair<int,ll>
#define pli pair<ll,int>
#define pii pair<int,int>
#define ull unsigned long long
#define bpc __builtin_popcount
#define base 1000000000000000000ll
#define fin freopen("a.txt","r",stdin)
#define fout freopen("a.txt","w",stdout)
#define fio ios::sync_with_stdio(false);cin.tie(0)
#define mr mt19937 rng(chrono::steady_clock::now().time_since_epoch().count())
inline ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
inline void sub(ll &a,ll b){a-=b;if(a<0)a+=mod;}
inline void add(ll &a,ll b){a+=b;if(a>=mod)a-=mod;}
template<typename T>inline T const& MAX(T const &a,T const &b){return a>b?a:b;}
template<typename T>inline T const& MIN(T const &a,T const &b){return a<b?a:b;}
inline ll qp(ll a,ll b){ll ans=1;while(b){if(b&1)ans=ans*a%mod;a=a*a%mod,b>>=1;}return ans;}
inline ll qp(ll a,ll b,ll c){ll ans=1;while(b){if(b&1)ans=ans*a%c;a=a*a%c,b>>=1;}return ans;}

using namespace std;
//using namespace __gnu_pbds;

const ull ba=233;
const db eps=1e-5;
const ld pi=acos(-1);
const ll INF=0x3f3f3f3f3f3f3f3f;
const int N=270000+10,maxn=2000000+10,inf=0x3f3f3f3f;

struct cd{
    ld x,y;
    cd(ld _x=0.0,ld _y=0.0):x(_x),y(_y){}
    cd operator +(const cd &b)const{
        return cd(x+b.x,y+b.y);
    }
    cd operator -(const cd &b)const{
        return cd(x-b.x,y-b.y);
    }
    cd operator *(const cd &b)const{
        return cd(x*b.x - y*b.y,x*b.y + y*b.x);
    }
    cd operator /(const db &b)const{
        return cd(x/b,y/b);
    }
}a[N],b[N],dfta[N],dftb[N],dftc[N],dftd[N];
cd conj(cd a){return cd(a.x,-a.y);}
int rev[N],A[N],B[N],C[N];
void getrev(int bit)
{
    for(int i=0;i<(1<<bit);i++)
        rev[i]=(rev[i>>1]>>1) | ((i&1)<<(bit-1));
}
void fft(cd *a,int n,int dft)
{
    for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
    for(int step=1;step<n;step<<=1)
    {
        cd wn(cos(dft*pi/step),sin(dft*pi/step));
        for(int j=0;j<n;j+=step<<1)
        {
            cd wnk(1,0);
            for(int k=j;k<j+step;k++)
            {
                cd x=a[k];
                cd y=wnk*a[k+step];
                a[k]=x+y;a[k+step]=x-y;
                wnk=wnk*wn;
            }
        }
    }
    if(dft==-1)for(int i=0;i<n;i++)a[i]=a[i]/n;
}
void mtt(int n,int m,int p) {
    int sz=0;
    while((1<<sz)<=n+m)sz++;getrev(sz);
    int len=1<<sz;
    for(int i=0;i<len;i++)
    {
        int x=(i>n?0:A[i]%p),y=(i>m?0:B[i]%p);
        a[i]=cd(x&0x7fff,x>>15);
        b[i]=cd(y&0x7fff,y>>15);
    }
    fft(a,len,1);fft(b,len,1);
    for(int i=0;i<len;i++)
    {
        int j=(len-i)&(len-1);
        cd aa,bb,cc,dd;
        aa = (a[i] + conj(a[j])) * cd(0.5, 0);
        bb = (a[i] - conj(a[j])) * cd(0, -0.5);
        cc = (b[i] + conj(b[j])) * cd(0.5, 0);
        dd = (b[i] - conj(b[j])) * cd(0, -0.5);
        dfta[j] = aa * cc;dftb[j] = aa * dd;
        dftc[j] = bb * cc;dftd[j] = bb * dd;
    }
    for(int i=0;i<len;i++)
    {
        a[i] = dfta[i] + dftb[i] * cd(0, 1);
        b[i] = dftc[i] + dftd[i] * cd(0, 1);
    }
    fft(a,len,1);fft(b,len,1);
    for(int i=0;i<len;i++)
    {
        int da = (ll)(a[i].x / len + 0.5) % p;
        int bb = (ll)(a[i].y / len + 0.5) % p;
        int dc = (ll)(b[i].x / len + 0.5) % p;
        int dd = (ll)(b[i].y / len + 0.5) % p;
        C[i] = (da + ((ll)(bb + dc) << 15) + ((ll)dd << 30)) % p;
        C[i] = (C[i]+p)%p;
    }
}
int main()
{
    int n,m,p;scanf("%d%d%d",&n,&m,&p);
    for(int i=0;i<=n;i++)scanf("%d",&A[i]);
    for(int i=0;i<=m;i++)scanf("%d",&B[i]);
    mtt(n,m,p);
    for(int i=0;i<=n+m;i++)printf("%d ",C[i]);
    return 0;
}
/********************

********************/

原文地址:https://www.cnblogs.com/acjiumeng/p/11109702.html

时间: 2024-07-31 13:04:05

任意模数fft的相关文章

MTT:任意模数NTT

MTT:任意模数NTT 概述 有时我们用FFT处理的数据很大,而模数可以分解为\(a\cdot 2^k+1\)的形式.次数用FFT精度不够,用NTT又找不到足够大的模数,于是MTT就应运而生了. MTT没有模数的限制,比NTT更加自由,应用广泛,可以用于任意模数或很大的数. MTT MTT是基于NTT的,其思想很简单,就是做多次NTT,每次使用不同的素数,然后使用CRT合并解,在合并的过程中模最终模数,或是对于无模数的情况使用高精度. 做NTT的次数取决于最大可能答案的大小,所用的所有素数之积必

任意模数NTT(MTT)模板

记住代码里3个模数,它们的原根都是3.考虑通过3个模数下的答案用中国剩余定理乱搞,得出答案.常数较大.有个什么拆系数FFT不会. P4245 [模板]任意模数NTT #include<bits/stdc++.h> using namespace std; typedef long long LL; #define rint register int const int mod[3]={469762049,998244353,1004535809}; const int N=400010; in

【51nod】1123 X^A Mod B (任意模数的K次剩余)

题解 K次剩余终极版!orz 写一下,WA一年,bug不花一分钱 在很久以前,我还认为,数论是一个重在思维,代码很短的东西 后来...我学了BSGS,学了EXBSGS,学了模质数的K次剩余--代码一个比一个长-- 直到今天,我写了240行的数论代码,我才发现数论这个东西= =太可怕了 好吧那么我们来说一下任意模数的K次剩余怎么搞 首先,如果模数是奇数,我们可以拆成很多个质数的指数幂,再用解同余方程的方法一个个合起来,细节之后探讨 但是如果,模数有偶数呢 因为要输出所有解,解的个数不多,我们可以倍

51nod 1172 Partial Sums V2 任意模FFT

Partial Sums V2 题目连接: https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1172 Description 给出一个数组A,经过一次处理,生成一个 数组S,数组S中的每个值相当于数组A的累加,比如:A = {1 3 5 6} => S = {1 4 9 15}.如果对生成的数组S再进行一次累加操作,{1 4 9 15} => {1 5 14 29},现在给出数组A,问进行K次操作后的结果.(输出结果

洛谷P4245 【模板】MTT(任意模数NTT)

题目背景 模板题,无背景 题目描述 给定 22 个多项式 F(x), G(x)F(x),G(x) ,请求出 F(x) * G(x)F(x)∗G(x) . 系数对 pp 取模,且不保证 pp 可以分解成 p = a \cdot 2^k + 1p=a⋅2k+1 之形式. 输入输出格式 输入格式: 输入共 33 行.第一行 33 个整数 n, m, pn,m,p ,分别表示 F(x), G(x)F(x),G(x) 的次数以及模数 pp .第二行为 n+1n+1 个整数, 第 ii 个整数 a_iai?

luogu P4245 【模板】任意模数NTT MTT

Code: #include<bits/stdc++.h> #define setIO(s) freopen(s".in","r",stdin) #define maxn 1000000 #define M 32768 #define double long double #define ll long long using namespace std; namespace poly{ const double pi=acos(-1); int rev[

[题解] Luogu P4245 [模板]任意模数NTT

三模NTT 不会... 都0202年了,还有人写三模NTT啊... 讲一个好写点的做法吧: 首先取一个阀值\(w\),然后把多项式的每个系数写成\(aw + c(c < w)\)的形式,换句话说把多项式\(f(x)\)写成两个多项式相加的形式: \[ f(x) = wf_0(x) + f_1(x) \] 这样在这道题中取\(W = 2^{15}\)就可以避免爆long long了. 乘起来的话就是 \[ f \cdot g = (w f_0 + f_1)(wg_0 + g_1) = (f_0 g

【codeforces 623E】dp+FFT+快速幂

题目大意:用$[1,2^k-1]$之间的证书构造一个长度为$n$的序列$a_i$,令$b_i=a_1\ or\ a_2\ or\ ...\ or a_i$,问使得b序列严格递增的方案数,答案对$10^9+7$取模. 数据范围,$n≤1^{18}$,$k≤30000$. 考虑用dp来解决这一题,我们用$f[i][j]$来表示前$i$个数中,使用了$j$个二进制位(注意!并不是前$j$个),那么答案显然为$\sum_{i=0}^{k} \binom{n}{i} \times f[n][i]$. 考虑

[拉格朗日反演][FFT][NTT][多项式大全]详解

1.多项式的两种表示法 1.系数表示法 我们最常用的多项式表示法就是系数表示法,一个次数界为\(n\)的多项式\(S(x)\)可以用一个向量\(s=(s_0,s_1,s_2,\cdots,s_n-1)\)系数表示如下:\[S(x)=\sum_{k=0}^{n-1}s_kx^k\] 系数表示法很适合做加法,可以在\(O(n)\)的时间复杂度内完成,表达式为:\[S(x)=A(x)+B(x)=\sum_{k=0}^{n-1}(a_k+b_k)x^k\] 当中\[s_k=a_k+b_k\] 但是,系数