FFT/ NTT

1.FFT

首先，\(FFT\)是用来求出\(h(x)=\sum_{i=1}^nf(i)*g(n-i)\)这样形式的卷积

\(h(x)\)中的每一项，用非常朴素的做法。两个多项式：

\[f(x)=a_0+a_1x+a_2x^2...a_nx^n\]

\[g(x)=b_0+b_1x+b_2x^2...b_mx^m\]

对于这两个多项式相乘，求出来每\(k\)次项前的系数就是算出来的系数，这个复杂

度可想而知\(O(n*m)\)

于是我们可以考虑这是两个分别具有\(n,m\)个系数的多项式，我们于是考虑代入求

值将自变量相同的\(f,g\)相乘，得到项数为\(n+m\)的点值式。每次选值代入，复杂度

\(O(n*n)\)

于是引入一个复数得到概念\(a+bi\)其中\(i=\sqrt{-1}\).把实部看作\(x\)方向，虚部看

作\(y\)方向的。复数就是和在一起的一个向量。

1.复数相加，直接实部加实部，虚部加虚部。

2.复数相乘，幅角相加，模长相乘。

\(c=r_1(cos\alpha+i*sin\alpha)\)

\(d=r_2(cos\beta+i*sin\beta)\)

\(c*d=r_1*r_2*(cos\alpha*cos\beta-sin\alpha*sin\beta+i*sin\alpha*cos\beta+i*sin\beta*cos\alpha)\)

于是：

\(c*d=r_1*r_2*(cos(\alpha+\beta)+i*sin(\alpha+\beta))\)

其中\(W_n=sin(2*\pi/n)+i*cos(2*\pi/n)\)

\(W_n^k\)就是\(W_n\)的\(k\)次方（好像不太妥）

1.\(W_n^0=W_n^n = 1\)

2.\(W_{n/2}^{k/2}=W_n^k\)

3.\(W_{n}^{k+n/2}=-W_n^k\)

于是对于一个多项式(\(n\)为偶数)

\[A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1}+a_nx^n\]

将其奇偶分类：

\[A_0(x)=a_0+a_2*x+a_4*x^2...+a_n*x^{n/2}\]

\[A_1(x)=a_1+a_3*x+a_5*x^2...+a_{n-1}*x^{n/2-1}\]

\[A(x)=A_1(x^2)*x+A_0(x^2)\]

将\(W_n^k,k\in[0,n)\),

1.\(k<n/2\)

\[A(W_n^k)=A_1(W_n^{2k})*W_n^k+A_0(W_n^{2k})\]

于是

\[A(W_n^k)=A_1(W_{n/2}^{k})*W_{n}^k+A_0(W_{n/2}^{k})\]

2.\(k>=n/2\)(此\(k\)非比下面的\(k\))

\[A(W_n^{k+n/2}) =-A_1(W_n^{2k})*W_n^k+A_0(W_n^{2*k})\]

于是

\[A(W_n^{k+n/2})=-A_1(W_{n/2}^{k})*W_n^k+A_0(W_{n/2}^k)\]

似乎前面两项只有一个符号的区别，再此观察奇偶分类以后的两个多项式，他们

可以继续递归分解，而得出来的系数恰好可以用来计算原来的多项式点值。

然后求出来的点值最后相乘的出来的新的点值在带回去求出现在的系数。

最后要记得除\(limit\)

至今只会递归写法：

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std ;
int n , m ;
#define N 3000000
inline int read()
{
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}
const double Pi = acos(-1.0) ;
struct node{
double x , y ;
node (double xx = 0 , double yy = 0 ){ x = xx , y = yy ; }
}a[N] , b[N] ;
node operator + (node a , node b){ return node(a.x + b.x , a.y + b.y) ; }
node operator - (node a , node b){ return node(a.x - b.x , a.y - b.y) ; }
node operator * (node a , node b){ return node(a.x * b.x - a.y*b.y , a.x *b.y + a.y * b.x ) ; }
void FFT(int limit , node *a , int type ){
    if( limit == 1 ) return ;
    node a1[(limit>>1)+30] , a2[(limit>>1) + 30] ;
    for(int i = 0 ; i <= limit ; i += 2 )
        a1[i>>1] = a[i] , a2[i>>1] = a[i + 1] ;
    FFT( limit >> 1 , a1 , type ) ;
    FFT( limit >> 1 , a2 , type ) ;
    node W = node( cos( 2.0 * Pi / limit ) , type * sin( 2.0 * Pi / limit ) ) ;
    node beg = node( 1 , 0 ) ;
    for(int i = 0 ; i < ( limit >> 1 ) ; i++ , beg = beg * W )
        a[i] = a1[i] + beg * a2[i] ,
        a[i + (limit>>1)] = a1[i] - beg * a2[i] ;
    return ;
}
int main()
{
    n = read() , m = read() ;
    for(int i = 0 ; i <= n ; i++ ) a[i].x = read() ;
    for(int i = 0 ; i <= m ; i++ ) b[i].x = read() ;
    int limit = 1 ;
    while( limit <= n + m ) limit <<= 1 ;
    FFT( limit , a , 1 ) ;
    FFT( limit , b , 1 ) ;
    for(int i = 0 ; i <= limit ; i++ ) a[i] = a[i] * b[i] ;
    FFT( limit , a , -1 ) ;
    for(int i = 0 ; i <= n + m ; i++ ) printf("%d " , (int)(a[i].x / limit + 0.5) ) ;
    return 0 ;
}

1.P3338 [ZJOI2014]力

\[E_i=\frac{F_i}{q_i}=\sum_{j<i}\frac{q_j}{(i-j)^2}+\sum_{j>i}\frac{q_j}{(i-j)^2}\]

令

\[p_i=q_{n-i+1},t_i=\frac{1}{i*i}\]

\(So\)

\[E_i=\sum_{j=1}^{i-1}t_j*p_{i-j}-\sum_{j=1}^{n-i}t_j*p_{n-j-i+1}\]

上面就是一个简单的卷积形式

1. NTT

由于\(FFT\)最大的毛病就是存在精度误差。于是我们能不能不用\(W_n\)这种存在精

度误差，而却依然满足某些性质的东东。

有个叫原根 的好东西

原根，是一个数学符号。设m是正整数，a是整数，若a模m的阶等于\(\phi(m)\)，则

称a为模m的一个原根。

食用\(NTT\)的话模数\(P\)一般取\(k*2^n+1\)，且\(P\)为质数

假如我们现在找到了一个在模\(p\)意义下的原根\(g\)。其中\(P=2^b*k+1,n=2^b\)

令\(g_n\equiv g^k\pmod{P}\)

于是

\[g_n^n\equiv g^{nk}\equiv g^{P-1}\equiv1\pmod{P}\]

由于原根本就具有的性质，且\(g_n^{nk}\)中\(nk<P\)

所以

\[g_n,g_n^1,g_n^2,g_n^3...\]

互不相等。

\[g_n^{2*k}\equiv g_{n/2}^k \pmod P\]

我们知道：

\[g^{P-1} \equiv 1 \pmod{P}\]

So

\[(g^{(P-1)/2}+1)*(g^{(P-1)/2}-1)\equiv 0 \pmod{P}\]

但是，由于

\[g^{(P-1)/2}\not\equiv1 \pmod{P} \]

所以：

\[g^{(P-1)/2}\equiv -1 \pmod P\]

end:

\[g_{n}^{k+n/2}\equiv g_n^k*g_n^{n/2}\equiv -g_n^k \pmod P\]

好的当我们要把\(NTT\)当\(FFT\)用时，可以取一个较大的模数，且其中\(b\)(上面那

个隐藏的卧底)合适

贴代码，风格突变，不管。。。

依然只会递归版

#include<bits/stdc++.h>
using namespace std ;

inline int read()
{
    char c=getchar();int x=0,f=1;
    while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
    while(c>='0'&&c<='9'){x=x*10+c-'0';c=getchar();}
    return x*f;
}

#define Mod 998244353
#define G 3
#define invG 332748118
#define N 10000000
#define ll long long

ll a[N] , b[N] ;
int n , m ;

ll KSM( ll k , ll x ){
    ll ans = 1 ;

    while( x ) {
        if( x & 1 ) ans = ans * k % Mod ;
        k = k * k % Mod ;
        x = x >> 1 ;

    }

    return ans % Mod ;
}
void NTT(int limit , ll *a , int type ){

    if( limit == 1 ) return ;

    ll a1[(limit>>1) + 30] , a2[(limit>>1) + 30] ;
    for(int i = 0 ; i <= limit ; i += 2 ){
        a1[i >> 1] = a[i] ;
        a2[i >> 1] = a[i + 1] ;

    }
    NTT( limit >> 1 , a1 , type ) ;
    NTT( limit >> 1 , a2 , type ) ;

    ll Wn = KSM( type == 1 ? G : invG , ( Mod - 1 ) / limit ) ;
    ll Gn = 1 ;

    for(int i = 0 ; i < ( limit >> 1 ) ; i++ , Gn = Gn * Wn % Mod ){
    ll x = a1[i] % Mod , y = a2[i] * Gn % Mod ;
    a[ i ] = ( x + y + Mod ) % Mod ;
    a[ i + ( limit >> 1 ) ] = ( x - y + Mod ) % Mod ;

    }
    return ;

}

int main()
{

    n = read() , m = read() ;
    for(int i = 0 ; i <= n ; i++ ) a[i] = read() ;
    for(int i = 0 ; i <= m ; i++ ) b[i] = read() ;

    int limit = 1 ;
    while( limit <= n + m ) limit = limit << 1 ;

    NTT( limit , a , 1 ) ;
    NTT( limit , b , 1 ) ;

    for(int i = 0 ; i <= limit ; i++ ) a[i] = ( a[i] * b[i] ) % Mod ;
    NTT( limit , a , -1 ) ;

    ll inv = KSM( limit , Mod - 2 ) ;
    for(int i = 0 ; i <= n + m ; i++ ) printf("%d " , ( a[i] * inv ) % Mod ) ;

    return 0 ;
}

转载于:https://www.cnblogs.com/powerYao/p/11445073.html