「Luogu P4245」【模板】任意模数NTT

发表于 |
字数统计: 962 | 阅读时长 ≈ 5

给定两个$n$次多项式$F(x),G(x)$的系数数列$a,b$和一个整数$p$，求$F(x)G(x)$在$\bmod p$意义下的值，不保证$p$可以分解成$a\cdot 2^k+1$的形式。

Data Range：$1\leq n\leq 10^5,0\leq a_i,b_i\leq 10^9,2\leq p\leq 10^9+9$

链接

Luogu P4245

题解

如果$p$是形如$a\cdot 2^k+1$的质数，可以直接刚$\texttt{NTT}$，但是布星啊，样例中的$p$都不是质数，所以说怎么办呢？

用$\texttt{NTT}$计算多项式乘法的最终结果是取模后的，所以考虑采用$\texttt{CRT}$合并答案。

但是，答案要在$\texttt{long long}$范围内唯一，所以取两个模数布星，要取三个，即所谓的三模数$\texttt{NTT}$。

这里我个人倾向于取$p_1=469762049$，$p_2=998244353$和$p_3=1004535809$，因为这三个数的原根是$3$。

于是可以通过$\texttt{NTT}$算出所给的多项式在模$p_1,p_2,p_3$意义下的乘积。

坑。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
#include<bits/stdc++.h>
using namespace std;
typedef long long int ll;
const ll MAXN=4e5+51,G=3;
const ll MOD1=469762049,MOD2=998244353,MOD3=1004535809;
const ll MOD=468937312667959297;
struct Polynomial{
    ll deg;
    ll coeff[MAXN];
};
Polynomial x,y,res1,res2,res3;
ll cnt,mod;
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
}
inline ll qadd(ll x,ll y,ll mod)
{
    return x+y>mod?x+y-mod:x+y;
}
inline ll qmin(ll x,ll y,ll mod)
{
    return x-y<0?x-y+mod:x-y;
}
inline ll qmul(ll x,ll y,ll mod)
{
    x%=mod,y%=mod;
    return ((x*y-(ll)((ll)((long double)x/mod*y+0.5)*mod))%mod+mod)%mod;
}
inline ll qpow(ll base,ll exponent,ll mod)
{
    if(!exponent)
    {
        return 1;
    }
    ll temp=qpow(base,exponent>>1,mod),res=temp*temp%mod;
    if(exponent&1)
    {
        res=res*base%mod;
    }
    return res;
}
inline void NTT(ll *cp,ll cnt,ll inv,ll mod)
{
    ll lim=0,cur=0,res=0,omg=0;
    while((1<<lim)<cnt)
    {
    	lim++;
    }
    for(register int i=0;i<cnt;i++)
    {
        cur=0;
        for(register int j=0;j<lim;j++)
        {
            if((i>>j)&1)
            {
                cur|=(1<<(lim-j-1));
            }
        }
        if(i<cur)
        {
            swap(cp[i],cp[cur]);
        }
    }
    for(register int i=2;i<=cnt;i<<=1)
    {
        cur=i>>1,res=qpow(G,(mod-1)/i,mod);
        for(register ll *p=cp;p!=cp+cnt;p+=i)
        {
            omg=1;
            for(register int j=0;j<cur;j++)
            {
                ll t=omg*p[j+cur]%mod;
                p[j+cur]=qmin(p[j],t,mod);
                p[j]=qadd(p[j],t,mod);
                omg=omg*res%mod;
            }
        }
    }
}
inline Polynomial mul(Polynomial x,Polynomial y,ll mod)
{
    Polynomial res;
    ll cnt=1,inv;
    static ll cpx[MAXN],cpy[MAXN];
    memset(cpx,0,sizeof(cpx)),memset(cpy,0,sizeof(cpy));
    for(register int i=0;i<=x.deg;i++)
    {
        cpx[i]=x.coeff[i];
    }
    for(register int i=0;i<=y.deg;i++)
    {
        cpy[i]=y.coeff[i];
    }
    while(cnt<=x.deg+y.deg)
    {
        cnt<<=1;
    }
    NTT(cpx,cnt,1,mod),NTT(cpy,cnt,1,mod);
    for(register int i=0;i<=cnt;i++)
    {
        cpx[i]=cpx[i]*cpy[i]%mod;
    }
    NTT(cpx,cnt,-1,mod);
    res.deg=x.deg+y.deg,inv=qpow(cnt,mod-2,mod);
    cpx[0]=cpx[0]*inv%mod;
    for(register int i=1;i<=cnt>>1;i++)
    {
        cpx[i]=cpx[i]*inv%mod;
        if(i!=cnt-i)
        {
            cpx[cnt-i]=cpx[cnt-i]*inv%mod;
        }
        swap(cpx[i],cpx[cnt-i]);
    }
    for(register int i=0;i<=res.deg;i++)
    {
        res.coeff[i]=qmin(cpx[i],0,mod);
    }
    return res;
}
inline ll CRT(ll r1,ll r2,ll r3,ll mod)
{
    ll inv1=208783132,inv2=395249030;
    ll r=qadd(qmul(qmul(qmin(r1,r2,MOD1),inv1,MOD1),MOD2,MOD),r2,MOD);
    ll k=qmul(qmin(r3,r,MOD3),inv2,MOD3);
    return ((k%mod)*(MOD%mod)%mod+r)%mod;
}
int main()
{
    x.deg=read(),y.deg=read(),mod=read();
    for(register int i=0;i<=x.deg;i++)
    {
        x.coeff[i]=read()%mod;
    }
    for(register int i=0;i<=y.deg;i++)
    {
        y.coeff[i]=read()%mod;
    }
    res1=mul(x,y,MOD1),res2=mul(x,y,MOD2),res3=mul(x,y,MOD3);
    for(register int i=0;i<=res1.deg;i++)
    {
        printf("%lld ",CRT(res1.coeff[i],res2.coeff[i],res3.coeff[i],mod));
    }
}