[模板] 快速数论变换 (NTT)

原根

原根的计算

$$\forall_{i \in [1, p-2]} g^i \not \equiv 1$$

$$g^{d^x} \equiv g^{dx} \equiv g^{p-1} \equiv 1 (mod \ p)$$

$$\forall g^q \not \equiv 1 (mod \ p)$$

$$p – 1 = \prod_{i=1}^r p_{i}^{k_i}$$

$$\forall_{i \in [1, r]} g^{\frac{p-1}{p_i}} \not \equiv 1 (mod \ p)$$

inline int root(int x){
for(register int i=2; i<=x; i++){
int tmp = x - 1;

bool flag = true;

for(register int k=2; k * k <= (x - 1); k++){
if(tmp % k == 0){
if(powx(i, (x - 1) / k, x) == 1){
flag = false;
break;
}

while(tmp % k == 0)
tmp /= k;
}
}

if(flag && (tmp == 1 || powx(i, (x - 1) / tmp, x) != 1)){
return i;
}
}
}


原根的性质

性质三

$$\omega_{n}^{n} \equiv g^{p-1} \equiv 1 (mod \ p)$$

$$\omega_{n}^{\frac{n}{2}} \equiv \pm 1 (mod \ p)$$

$$\omega_{n}^{\frac{n}{2}} \not \equiv \omega_{n}^{0} (mod \ p)$$

$$\omega_{n}^{\frac{n}{2}} \equiv – 1 (mod \ p)$$

$$\omega_{n}^{k + \frac{n}{2}} \equiv – \omega_{n}^{k} (mod \ p)$$

性质四

\begin{aligned} S(\omega_n^k) & = 1 + \omega_n^k + (\omega_n^k)^2 + \cdots + (\omega_n^k)^{n-1} \\ & = \frac{1-(\omega_n^k)^n}{1-\omega_n^k} \\ & = \frac{(\omega_n^k)^n – 1}{\omega_n^k – 1} \\ \end{aligned}

$${(\omega_{n}^{k})}^{n} – 1 \equiv 0 (mod \ p)$$

快速数论变换

$NTT$ 常用的模数为 $998244353$ 和 $1004535809$ ，它们的原根为 $3$

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;
const int MOD = 998244353;
const int MAXN = 4194304 + 10;

inline long long powx(long long a, long long b, long long mod){
long long ans = 1;

for(; b; b >>= 1){
if(b & 1) ans *= a, ans %= mod;
a *= a, a %= mod;
}

return ans;
}

inline long long exgcd(long long a, long long b, long long &x, long long &y){
if(!b){
x = 1, y = 0;
return a;
}

long long d = exgcd(b, a%b, y, x);
y -= x * (a / b);

return d;
}

inline long long inv(long long a, long long mod){
long long x, y;
exgcd(a, mod, x, y);
return (x + mod) % mod;
}

inline int root(int x){
for(register int i=2; i<=x; i++){
int tmp = x - 1;

bool flag = true;

for(register int k=2; k * k <= (x - 1); k++){
if(tmp % k == 0){
if(powx(i, (x - 1) / k, x) == 1){
flag = false;
break;
}

while(tmp % k == 0)
tmp /= k;
}
}

if(flag && (tmp == 1 || powx(i, (x - 1) / tmp, x) != 1)){
return i;
}
}
throw;
}

namespace NTT{
int place[MAXN];
long long omega[MAXN];
long long omegaInverse[MAXN];

inline void init(int n){
int k = 0;

while((1 << k) < n)
k++;

for(register int i=0; i<n; i++){
for(register int j=0; j<k; j++){
if(i & (1 << j)){
place[i] |= 1 << (k - j - 1);
}
}
}

long long g = root(MOD);
long long tmp = powx(g, (MOD - 1) / n, MOD);

for(register int i=0; i<n; i++){
omega[i] = (i == 0) ? 1 : omega[i - 1] * tmp % MOD;
omegaInverse[i] = inv(omega[i], MOD);
}
}

inline void transform(long long *a, int n, long long *omega){
for(register int i=0; i<n; i++){
if(i < place[i]){
swap(a[i], a[place[i]]);
}
}

for(register int range=2; range <= n; range <<= 1){
int mid = range >> 1;

for(register long long *p = a; p != a + n; p += range){
int k = n / range;

for(register int i=0; i<mid; i++){
int tmp = i + mid;

long long t = omega[k * i] * p[i + mid] % MOD;

p[tmp] = (p[i] - t + MOD) % MOD;
p[i] = (p[i] + t) % MOD;
}
}
}
}

inline void dft(long long *a, int n){
transform(a, n, omega);
}

inline void idft(long long *a, int n){
transform(a, n, omegaInverse);

long long tmp = inv(n, MOD);

for(register int i=0; i<n; i++){
a[i] = a[i] * tmp % MOD;
}
}
}

inline int read(){
int x = 0;
int p = 1;
char ch = getchar();

while(ch < '0' || ch > '9'){
if(ch == '-')
p = 0;
ch = getchar();
}

while('0' <= ch && ch <= '9'){
x = x*10 + ch - '0';
ch = getchar();
}

return p ? x : (-x);
}

int n, m, tot;

long long a[MAXN], b[MAXN];

int main(){
n = read();
m = read();

n++;
m++;

for(register int i=0; i<n; i++)
a[i] = read();

for(register int i=0; i<m; i++){
b[i] = read();
}

tot = 1;
while(tot < n + m){
tot <<= 1;
}

NTT::init(tot);

NTT::dft(a, tot);
NTT::dft(b, tot);

for(register int i=0; i<tot; i++){
a[i] = a[i] * b[i] % MOD;
}

NTT::idft(a, tot);

for(register int i=0; i<n+m-1; i++){
printf("%lld ", a[i]);
}
return 0;
}