快速傅里叶变换

最新推荐文章于 2024-08-12 16:55:46 发布

Mongo_w

最新推荐文章于 2024-08-12 16:55:46 发布

阅读量239

点赞数

分类专栏：笔记文章标签： c++ 算法机器学习

本文链接：https://blog.csdn.net/weixin_45916935/article/details/112578903

版权

笔记专栏收录该内容

15 篇文章 0 订阅

订阅专栏

FFT

最后的最后还是回到了这个玩意上…

前置

单位根

记为 $w_n$ ，表示将圆 $n$ 等分后的第一块对应的复数.

于是我们可以写出 $w_n=\cos \frac{2\pi}{n}+i\sin \frac{2\pi}{n}$ .

同样可以定义 $k$ 次方，根据复数的三角运算可知：

$w_n^k=\cos k\times \frac{2\pi}{n}+i\sin k\times \frac{2\pi}{n}$ .

$w_{2n}^{2k}=w_n^{k}$ .

注意到 $w_{2n}^{2k}=\cos 2k\times \frac{2\pi}{2n}+i\sin 2k\times \frac{2\pi}{2n}=\cos k\times \frac{2\pi}{n}+i\sin k\times \frac{2\pi}{n}$ .

$w_n^{k+\frac{n}{2}}=w_{n}^{k}$ .

显然 $w_{n}^{\frac{n}{2}}=\cos \frac{n}{2}\times \frac{2\pi}{n}+i\sin \frac{n}{2}\times \frac{2\pi}{n}=-1+0=-1$ ，乘起来即可.

快速傅里叶变换(FFT)

这一步需要我们将系数表达式变成点值表达式.

设 $F(x)=\sum\limits_{k=0}^n a_kx^k$ 为我们要处理的 $n$ 次多项式，我们的操作是将他们的系数奇偶分类：

$F(x)=(a_0+a_2x^2+\ldots)+x(a_1+a_3x^2+\ldots)$ .

注意到后面我们设成新的两个多项式 $A (x)$ 和 $B (x)$ ，则整体就可以写成:

$F(x)=A(x^2)+xB(x^2)$ .

把 $w_n^k$ 扔进去得到：

$F(w^k_n)=A(w_{\frac{n}{2}}^{k})+w_{n}^{k}B(w_{\frac{n}{2}}^k)$ .

然后还有一个：

$F(w_{n}^{k+\frac{n}{2}})=A(w_{\frac{n}{2}}^k)-w_{n}^{k}B(w_\frac{n}{2}^k)$ .

所以我们在处理 $k\in [0,\frac{n}{2}]$ 的时候可以顺便把 $k\in (\frac{n}{2},n]$ 的给弄出来.

递归处理一边即可，时间复杂度 $O(n\log n)$ .

至于怎么写，一下给出参考代码（不保证能过编译.

struct complex{
    double x , y;
    complex (double xx = 0 , double yy = 0){x = xx , y = yy;}
} a[MAXN] , b[MAXN] , tmp[MAXN];
/*
math operation
*/
void fft(complex *A , int l , int r)
{
    if(l == r)
    {
        
	}
    int cur = 0 , h = l , t = (l + r >> 1) + 1;
    for(int i = l ; i <= r ; i++)
    {
        if((i - l) & 1) tmp[t++] = A[i];
        else tmp[l++] = A[i];
	}
    fft(A , l , l + r >> 1);
    fft(A , (l + r >> 1) + 1 , r);
    for(int i = l ; i <= r ; i++)
    {
        /*
        fft operation
        */
	}
}

现在转成了点值表达式，但我需要的是系数表达式啊.

所以还有快速傅里叶逆变换.

快速傅里叶逆变换(IFFT)

我们现在得到了点值表达 $(y_0,y_1,\ldots,y_{n-1})$ ，构造另一组向量 $(c_0,c_1,\ldots,c_{n-1})$ 使其满足：

$c_k=\sum\limits_{i=0}^{n-1}y_i(w_n^{-k})^i$ ，即将 $c_k$ 看成以 $(y_0,y_1\ldots,y_{n-1})$ 为系数表达的多项式 $B (x)$ ，其取 $x=w_n^{-k}$ 时的值.

至于这是干什么的，先展开化简：

$c_k=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(w_{n}^{i})^j)(w_{n}^{-k})^i$

$=\sum\limits_{i=0}^{n-1}(\sum\limits_{j=0}^{n-1}a_j(w_{n}^j)^i(w_{n}^{-k})^i)$

$=\sum\limits_{i=0}^{n-1}\sum\limits_{j=0}^{n-1}a_j(w_{n}^{j-k})^i$

$=\sum\limits_{j=0}^{n-1}a_j(\sum\limits_{i=0}^{n-1}(w_n^{j-k})^i)$ .

记 $S(x)=\sum\limits_{i=0}^{n-1}x^i$ .

那上面也可以写成：

$c_k=\sum\limits_{j=0}^{n-1}a_jS(w_n^{j-k})$ .

而 $S(w_n^k)=1+w_n^k+(w_n^{k})^2+\cdots+(w_n^{k})^{n-1}$ .

发挥人类智慧，当 $k\not= 0$ 时，同时乘一个 $w_n^k$ 调项得到：

$1-w_{n}^k)S(w_n^k)=1-(w_n^k)^n$ .

$S(w_n^k)=\frac{1-1}{1-w_n^k}=0$ .

而当 $k = 0$ 时直接得到：

$S (1) = n$ .

所以后面那一坨可以扔掉，原式子写成：

$c_k=a_kn$ .

即 $a_k=\frac{c_k}{n}$ ，至此我们实现了点值到系数的转化.

实现把 $w_n^k$ 换成 $w_n^{-k}$ 即可，其余和 $\text{FFT}$ 是一样的.

实现常数太大T了怎么办？

我们又充分发挥人类智慧，观察 $\text{FFT}$ 的过程，操作前相当于将系数按一个次序排列，操作的时候则相当于直接覆盖（专业术语蝴蝶变换.

于是我们找到这个规律即可，由观察法，记 $p (x)$ 表示系数下标 $x$ 的最终下标，则

$p (x)$ 等价于将 $x$ 二进制翻转后的值.

预处理这货即可：

int l = 0 , lmt = 1 , p[N];
int main()
{
    while(lmt <= N + M) lmt <<= 1 , l++;
    for(int i = 0 ; i <= lmt ; i++)
        p[i] = (p[i >> 1] >> 1) | ((i & 1) << (l - 1));
}

那 $\text{FFT}$ 就这么写：

/*
structrue definition and complex math operation
*/
const double Pi = acos(-1.0);
void fft(complex *A , int type)
{
    for(int i = 0 ; i <= lmt ; i++)
        if(i < p[i]) swap(A[i] , A[p[i]]);
   	for(int mid = 1 ; mid < lmt ; mid <<= 1)//这是待合并的区间的一半
    {
        complex Wn(cos(pi / mid) , type * sin(pi / (mid << 1)));
        for(int R = mid << 1 , j = 0 ; j <= lmt ; j += R)//R是合并区间的总长度，j用来跳段
        {
            complex w(1 , 0);
            for(int k = 0 ; k < mid ; k++ , w = w * Wn)//枚举左半部分
            {
                complex x = A[j + k] , y = w * A[j + mid + k];//人类智慧
                A[j + k] = x + y; A[j + mid + k] = x - y;
            }
        }
    }
}

于是 $\text{FFT}$ 就介绍完了.

然而，信竞里边设计卷积计算方案数的基本都要带模数.

而涉及模数，就意味着会乘爆.

所以普通 $\text{FFT}$ 是不够的，需要用上数论快速傅里叶变换(NTT)，最终目的即摆脱单位根.

于是我们用上了原根…

前置

阶

对于 $(a, p) = 1$ ，定义 $a^n\equiv 1(\mod p)$ 的最小正整数 $n$ 为 $a$ 在模 $p$ 意义下的阶，记作 $\delta_p(a)$ .

例如 $\delta_7(3) = 6$ .

$3\mod 7=3$ .

$3^2\mod 7=2$ .

$3^3\mod 7=6$ .

$3^4\mod 7=4$ .

$3^5\mod 7=5$ .

$3^6\mod 7=1$ .

原根

Definition

设 $p$ 为正整数， $a$ 为整数，若 $\delta_p(a)=\phi(p)$ ，则称 $a$ 为模 $p$ 意义下的一个原根.

由上方知 $\delta_7(3)=\phi(7)=6$ ，所以 $3$ 是模 $7$ 意义下的一个原根.

其有一个性质：

若 $P$ 为一个素数，假设一个数 $g$ 是 $P$ 的原根，那么 $g^i\mod P$ 的结果两两不同.

完全剩余系的基本定义.

再者，我们发现 $g^k$ 的阶恰好等于 $\frac{p-1}{\gcd(k,p-1)}$ ，并且其为 $p - 1$ 的约数.

于是对于形如 $n=2^k$ 的， $p$ 取 $998244353=2^{23}\times 7\times 17 + 1$ 就能基本实现自律.

此时得到 $k=\frac{p-1}{n}$ ，对应原根 $g = 3$ .

至于怎么摆脱原根…只能祈祷出题人不会这么阴间了.

#define ll long long
const int p = 998244353 , G = 3;
inline int fpow(int x , int y = p - 2)
{
    int ans = 1;
    for(; y ; y >>= 1 , x = 1ll * x * x % p) if(y & 1) ans = 1ll * ans * x % p;
    return ans;
}
const int invG = fpow(G);
void NTT(ll *A , int type)
{
    for(int i = 0 ; i < lmt ; i++)
        if(i < r[i]) swap(A[i] , A[r[i]]);
    for(int mid = 1 ; mid < lmt ; mid <<= 1)
    {
        ll Wn = fpow(type == 1 ? G : invG , (p - 1) / (mid << 1));
        for(int R = (1 << mid) , j = 0 ; j <= lmt ; j += R)
        {
            ll w = 1;
            for(int k = 0 ; k < mid ; k++ , w = 1ll * w * Wn % mod)
            {
                ll x = A[j + k] , y = w * A[j + mid + k] % p;
                A[j + k] = (x + y) % p; A[j + mid + k] = (x - y) % p;
            }
        }
	}
}

以及提供一个封装库：

namespace Poly
{
	#define ll long long
	const int p = 998244353 , G = 3;
	inline ll fpow(ll x , ll y = p - 2)
	{
		ll ans = 1;
		for(; y ; y >>= 1 , x = x * x % p) if(y & 1) ans = ans * x % p;
		return ans;
	}
	const ll invG = fpow(G);
	int r[N] , l , lmt = 1;
	inline void prerev(int n , int m)
	{
		while(lmt <= n + m) l++ , lmt <<= 1;
		for(int i = 0 ; i <= lmt ; i++)
			r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
	}
	inline void NTT(ll *A , int type)
	{
		for(int i = 0 ; i <= lmt ; i++)
			if(i < r[i]) swap(A[i] , A[r[i]]);
		for(int mid = 1 ; mid < lmt ; mid <<= 1)
		{
			ll Wn = fpow(type == 1 ? G : invG , (p - 1ll) / (mid << 1));
			for(int R = mid << 1 , j = 0 ; j < lmt ; j += R)
			{
				ll w = 1;
				for(int k = 0 ; k < mid ; k++ , w = w * Wn % p)
				{
					ll x = A[j + k] , y = w * A[j + mid + k] % p;
					A[j + k] = (x + y) % p; A[j + mid + k] = (x - y + p) % p;
				}
			}
		}
	}
	inline void solve(ll *a , ll *b , int n , int m)
	{
		prerev(n , m);
		NTT(a , 1); NTT(b , 1);
		for(int i = 0 ; i <= lmt ; i++) a[i] = a[i] * b[i] % p;
		NTT(a , -1);
		ll inv = fpow(lmt);
		for(int i = 0 ; i <= lmt ; i++) a[i] = (a[i] * inv) % p;
	}
}
int main()
{
    /*
    */
    Poly::NTT(a , b , n , m);
    /*
    */
}

Mongo_w

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
快速傅里叶变换

这里有一个问题：我们有两个多项式，分别计作Ai(x),Bj(x)A_{i}(x),B_{j}(x)Ai(x),Bj(x)，这里的下标表示对应多项式最高阶数，如Ai(x)A_{i}(x)Ai(x)的最高阶数为iii即：Ai(x)=a0+a1x1+a2x2+……+aixiA_{i}(x)=a_{0}+a_{1}x^{1}+a_{2}x^{2}+……+a_{i}x^{i}Ai(x)=a0+a1x1+a2x2+……+aixiBj(x)=b0+b1x1+b2x2+……+bjxjB_{j}(x)
复制链接

扫一扫

专栏目录