MTT:任意模数NTT

MTT:任意模数NTT

概述

有时我们用FFT处理的数据很大,而模数可以分解为\(a\cdot 2^k+1\)的形式。次数用FFT精度不够,用NTT又找不到足够大的模数,于是MTT就应运而生了。

MTT没有模数的限制,比NTT更加自由,应用广泛,可以用于任意模数或很大的数。

MTT

MTT是基于NTT的,其思想很简单,就是做多次NTT,每次使用不同的素数,然后使用CRT合并解,在合并的过程中模最终模数,或是对于无模数的情况使用高精度。

做NTT的次数取决于最大可能答案的大小,所用的所有素数之积必须大于答案

实现

此处以取三个素数为例
我们可以做三次NTT,相邻次之间改变素数,但这样常数太大,于是我们常常选择封装(适合于模数不太多的情况)。
我们定义一个结构体node,有三个成员a,b,c,分别代表三个模数下的值,同时,我们定义模数的结构体与之一一对应。

struct node{
    LL a,b,c;
    node(){
        a=b=c=0;
    }
    node(LL x){
        a=b=c=x;
    }
    node(LL x,LL y,LL z){
        a=x;
        b=y;
        c=z;
    }
}MOD=node(167772161,469762049,998244353),BASE=node(3),INV=node(116878283,426037461,929031873);

我们还要定义关于此结构体的运算,其中成员之间互不影响,只和操作对象里对应的成员产生运算

inline node operator+(node x,node y){
    return node(x.a+y.a,x.b+y.b,x.c+y.c);
}
inline node operator-(node x,node y){
    return node(x.a-y.a,x.b-y.b,x.c-y.c);
}
inline node operator*(node x,node y){
    return node(x.a*y.a%MOD.a,x.b*y.b%MOD.b,x.c*y.c%MOD.c);
}
inline node operator%(node x,node y){
    return node(x.a%y.a,x.b%y.b,x.c%y.c);
}
inline node operator/(node x,node y){
    return node(x.a/y.a,x.b/y.b,x.c/y.c);
}
inline node operator-(node x,LL y){
    return node(x.a-y,x.b-y,x.c-y);
}
inline node operator*(node x,LL y){
    return node(x.a*y,x.b*y,x.c*y);
}
inline node operator/(node x,LL y){
    return node(x.a/y,x.b/y,x.c/y);
}
inline node operator%(node x,LL y){
    return node(x.a%y,x.b%y,x.c%y);
}

然后套用NTT的板子,最后用CRT合并。

假设这一位的答案是\(x\),三个模数分别为\(A,B,C\),那么:
\[\begin{aligned}x\equiv x_1\pmod{A} \\ x\equiv x_2\pmod{B} \\ x\equiv x_3\pmod{C}\end{aligned}\]
先把前两个合并:
\[\begin{aligned}x_1+k_1A=x_2+k_2B\\x_1+k_1A\equiv x_2\pmod{B}\\k_1\equiv \frac{x_2-x_1}A\pmod{B}\end{aligned}\]
于是求出了\(k_1\),也就求出了\(x\equiv x_1+k_1A\pmod{AB}\),记\(x_4=x_1+k_1A\)

\[\begin{aligned}x_4+k_4AB=x_3+k_3C\\x_4+k_4AB\equiv x_3\pmod{C}\\k_4\equiv \dfrac{x_3-x_4}{AB}\pmod{C}\end{aligned}\]
因为\(x<ABC\),所以
\[x=x_4+k_4AB\]

LL CRT(node x){
    LL mod1=MOD.a,mod2=MOD.b,mod3=MOD.c,mod_1_2=mod1*mod2;
    LL inv_1=inv(mod1,mod2),inv_2=inv(mod1*mod2%mod3,mod3);
    LL A=x.a,B=x.b,C=x.c;
    LL x4=(B-A+mod2)%mod2*inv_1%mod2*mod1+A;
    return ((C-x4%mod3+mod3)%mod3*inv_2%mod3*(mod_1_2%Last_Mod)+Last_Mod+x4)%Last_Mod;
}

于是我们就能写出完整代码了。

// luogu-judger-enable-o2
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const int INF=1e9+7,MAXN=3e6+10/*Min:2^20+10*/;
void exgcd(LL a,LL b,LL &x,LL &y){
    if(!b){
        x=1;
        y=0;
        return;
    }
    exgcd(b,a%b,y,x);
    y-=a/b*x;
}
inline LL inv(LL x,LL p){
    LL a,b;
    exgcd(x,p,a,b);
    return (a%p+p)%p;
}
struct node{
    LL a,b,c;
    node(){
        a=b=c=0;
    }
    node(LL x){
        a=b=c=x;
    }
    node(LL x,LL y,LL z){
        a=x;
        b=y;
        c=z;
    }
}MOD=node(167772161,469762049,998244353),BASE=node(3),INV=node(116878283,426037461,929031873);
inline node operator+(node x,node y){
    return node(x.a+y.a,x.b+y.b,x.c+y.c);
}
inline node operator-(node x,node y){
    return node(x.a-y.a,x.b-y.b,x.c-y.c);
}
inline node operator*(node x,node y){
    return node(x.a*y.a%MOD.a,x.b*y.b%MOD.b,x.c*y.c%MOD.c);
}
inline node operator%(node x,node y){
    return node(x.a%y.a,x.b%y.b,x.c%y.c);
}
inline node operator/(node x,node y){
    return node(x.a/y.a,x.b/y.b,x.c/y.c);
}
inline node operator-(node x,LL y){
    return node(x.a-y,x.b-y,x.c-y);
}
inline node operator*(node x,LL y){
    return node(x.a*y,x.b*y,x.c*y);
}
inline node operator/(node x,LL y){
    return node(x.a/y,x.b/y,x.c/y);
}
inline node operator%(node x,LL y){
    return node(x.a%y,x.b%y,x.c%y);
}
LL Last_Mod;
LL CRT(node x){
    LL mod1=MOD.a,mod2=MOD.b,mod3=MOD.c,mod_1_2=mod1*mod2;
    LL inv_1=inv(mod1,mod2),inv_2=inv(mod1*mod2%mod3,mod3);
    LL A=x.a,B=x.b,C=x.c;
    LL x4=(B-A+mod2)%mod2*inv_1%mod2*mod1+A;
    return ((C-x4%mod3+mod3)%mod3*inv_2%mod3*(mod_1_2%Last_Mod)+Last_Mod+x4)%Last_Mod;
}
inline LL fpm_(LL base,LL p,LL mod){
    LL ret=1;
    while(p){
        if(p&1)
            ret=ret*base%mod;
        base=base*base%mod;
        p>>=1;
    }
    return ret%mod;
}
inline node fpm(LL base,node p){
    return node(fpm_(base,p.a,MOD.a),fpm_(base,p.b,MOD.b),fpm_(base,p.c,MOD.c));
}
int N,M,lim=1,lg,rev[MAXN];
node Wn[MAXN];
inline void NTT(node *a,int type){
    for(int i=0;i<lim;i++)
        if(i<rev[i])
            swap(a[i],a[rev[i]]);
    for(int mid=1;mid<lim;mid<<=1){
        int len=mid<<1/*n*/;
        node Wn=fpm(3,(MOD-1)/(LL)len);
        for(int j=0;j<lim;j+=len){
            node w=node(1);
            for(int k=0;k<mid;k++){
                node x=a[j+k],y=w*a[j+k+mid]%MOD;
                a[j+k]=(x+y)%MOD;
                a[j+k+mid]=(x-y+MOD)%MOD;
                w=w*Wn%MOD;
            }
        }
    }
    if(type==-1){
        reverse(a+1,a+lim);
        node lim_inv=node(inv(lim,MOD.a),inv(lim,MOD.b),inv(lim,MOD.c));
        for(int i=0;i<lim;i++)
            a[i]=a[i]*lim_inv;
    }
}
node a[MAXN],b[MAXN];
int main(){
    scanf("%d%d%lld",&N,&M,&Last_Mod);
    for(int i=0;i<=N;i++){
        LL ii;
        scanf("%lld",&ii);
        a[i]=node(ii%Last_Mod)%MOD;
    }
    for(int i=0;i<=M;i++){
        LL ii;
        scanf("%lld",&ii);
        b[i]=node(ii%Last_Mod)%MOD;
    }
    while(lim<=N+M){
        lim<<=1;
        lg++;
    }
    for(int i=0;i<lim;i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(lg-1));
    Wn[0]=node(1);
    for(int i=1;i<lim;i++)
        Wn[i]=Wn[i-1]*INV;
    NTT(a,1);
    NTT(b,1);
    for(int i=0;i<lim;i++)
        a[i]=a[i]*b[i]%MOD;
    NTT(a,-1);
    for(int i=0;i<=N+M;i++)
        printf("%lld ",CRT(a[i]));
    return 0;
}

例题

模板题:P4245?【模板】任意模数NTT
模板题:【模板】多项式求逆(加强版)

原文地址:https://www.cnblogs.com/guoshaoyang/p/11296140.html

时间: 2024-10-01 07:45:50

MTT:任意模数NTT的相关文章

洛谷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?

任意模数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

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

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

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

三模数NTT模板

求两个多项式的卷积对任意数p取模 两个好记的FNT模数: 5*2^25+1 7*2^26+1 原根都为3 1 //Achen 2 #include<algorithm> 3 #include<iostream> 4 #include<cstring> 5 #include<cstdlib> 6 #include<vector> 7 #include<cstdio> 8 #include<queue> 9 #include&

任意模数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&q

Noip前的大抱佛脚----数论

数论 知识点 Exgcd \(O(logn)\)求解\(Ax+By=C\)的问题 1.若\(C\%gcd(A,B)!=0\)则无解 2.\(Gcd=gcd(A,B);A/=Gcd,B/=Gcd,C/=Gcd\) 3.代入下面代码求\(Ax+By=1\) 4.\(x*C\),得到一组特解 5.通解为\(\begin{cases}x=x_0+k*B \\y=y_0+k*A\end{cases}\) void Exgcd(ll a,ll b,ll &x,ll &y) { if(!b){x=1;y

多项式总结(持续更新)

多项式的一堆乱七八糟的操作学了一部分了--(多点求值和快速插值还没有) 打算写下来整理一下.不过因为还有一些没学的以及没完全理解的--只好先持续更新了. 不扯淡了,直接开始. 1.NTT FFT咱就不说了,有兴趣可以看兔哥博客. NTT和FFT很相似.但是因为FFT涉及到复数运算所以会有一些精度误差,然后有的时候也会遇到需要取模的情况--于是快速数论变换NTT应运而生.因为单位根和原根有相似的性质,所以NTT使用原根取代了单位根进行运算.模数998244353的原根是3,每次取原根的\(\fra