快速数论变换NTT

前言

如果你没有看过快速傅里叶变换FFT，那么建议先去了解一下再来看这篇文章。

首先恭喜你，背出了FFT模板学会了FFT。而我们发现，由于FFT所用的大量实数操作，有些时候无法满足精度需求。或者说，遇到模意义下的操作，FFT也无法满足需求。

所以就有了NTT。

原根

FFT能执行的原因就是巧妙地选择了单位复根。那么我们需要在模意义下找到单位复根的替代品。

我们首先看一看我们需要具有什么性质的数：

1、$\omega_n^n=1$。

2、$\omega_n^0,\omega_n^1,…,\omega_n^{n-1}$互不相同。

3、$\omega_{kn}^{ki}=\omega_n^i$。

4、$\omega_n^k=-\omega_n^{k+\frac{n}{2}}$。

5、$\sum_{i=0}^{n-1}\omega_n^{ti}={\Large\{}\begin{matrix}0&t\neq0\\n & t=0\end{matrix}$。

然而本蒟蒻太菜了，只能到处找资料。

下面的两项定义及原根的求法参考https://www.cnblogs.com/cytus/p/9296661.html

首先我们定义阶：

若$m>1$且$gcd(a,m)=1$，那么使得$a^r\equiv1(mod\,\,m)$的最小正整数$r$称为$a$对模$m$的阶，记做$\delta_m(a)$。

然后定义原根：

设$m$为正整数，$a$为整数，如果$\delta_m(a)=\phi(m)$，那么称$a$为模$m$的一个原根。

原根的求法：

设$\phi(m)=p_1^{e_1}p_2^{e_2}…p_n^{e_n}​$，若$g^{\frac{\phi(m)}{p_i}} \mod m \neq 1​$，那么$g​$是$m​$的一个原根。

常见的模数及其原根：http://blog.miskcoo.com/2014/07/fft-prime-table

原根性质的证明可以参考https://www.cnblogs.com/zhouzhendong/p/Fast-Fourier-Transform.html

参考代码

#include <cstdio>
#include <cstring>
#include <algorithm>
#define LL long long
using namespace std;

const LL Mod = 998244353;
const LL Maxn = 400010;
LL N, M, TotalLen, L;
LL A[ Maxn ], B[ Maxn ];
LL Index[ Maxn ], Omega[ Maxn ];
char Ch[ Maxn ];

void Exgcd( LL a, LL b, LL &x, LL &y ) {
    if( b == 0 ) {
        x = 1; y = 0; return;
    }
    Exgcd( b, a % b, y, x );
    y -= a / b * x;
    return;
}

LL Inv( LL a ) {
    LL x, y;
    Exgcd( a, Mod, x, y );
    if( x < 0 ) x += Mod;
    return x;
}

LL FastPow( LL x, LL y ) {
    if( y == 0 ) return 1;
    LL t = FastPow( x, y / 2 );
    t = t * t % Mod;
    if( y & 1 ) t = t * x % Mod;
    return t;
}

void Read() {
    scanf( "%s", Ch );
    N = strlen( Ch );
    for( LL i = 0; i < N; ++i ) A[ i ] = Ch[ N - 1 - i ] - '0';
    scanf( "%s", Ch );
    M = strlen( Ch );
    for( LL i = 0; i < M; ++i ) B[ i ] = Ch[ M - 1 - i ] - '0';
    return;
}

void Init() {
    TotalLen = N + M;
    for( L = 1; L < TotalLen; L <<= 1 );
    for( LL i = 0; i < L; ++i ) 
        Index[ i ] = ( Index[ i >> 1 ] >> 1 ) | ( ( i & 1 ) * ( L >> 1 ) );
    Omega[ 0 ] = 1; Omega[ 1 ] = FastPow( 3, ( Mod - 1 ) / L );
    for( LL i = 2; i < L; ++i ) 
        Omega[ i ] = Omega[ i - 1 ] * Omega[ 1 ] % Mod;
    return;
}

void NTT( LL *A ) {
    for( LL i = 0; i < L; ++i ) 
        if( Index[ i ] > i )
            swap( A[ i ], A[ Index[ i ] ] );
    for( LL HalfLen = 1; HalfLen < L; HalfLen <<= 1 ) 
        for( LL i = 0; i < L; i += HalfLen << 1 )
            for( LL j = 0; j < HalfLen; ++j ) {
                LL T = Omega[ L / 2 / HalfLen * j ] * A[ i + j + HalfLen ] % Mod;    
                LL t = A[ i + j ];
                A[ i + j ] = ( t + T ) % Mod;
                A[ i + j + HalfLen ] = ( t - T + Mod ) % Mod;
            }
    return;
}

void Transform() {
    for( LL i = 0; i < L; ++i ) A[ i ] = ( A[ i ] * B[ i ] ) % Mod;
    Omega[ 1 ] = Inv( Omega[ 1 ] );
    for( LL i = 2; i < L; ++i )
        Omega[ i ] = Omega[ i - 1 ] * Omega[ 1 ] % Mod;
    return;
}

void Print() {
    LL INV = Inv( L );
    for( LL i = 0; i < L; ++i ) 
        A[ i ] = A[ i ] * INV % Mod;
    for( LL i = 0; i < L; ++i ) {
        A[ i + 1 ] += A[ i ] / 10;
        A[ i ] %= 10;
    }
    while( L && !A[ L ] ) --L;
    for( LL i = L; i >= 0; --i ) printf( "%lld", A[ i ] );
    printf( "\n" );
    return;
}

int main() {
    Read();
    Init();
    NTT( A ); NTT( B );
    Transform();
    NTT( A );
    Print();
    return 0;
}