FFT 瞎讲

FFT

蒟蒻学了 FFT 赶紧做个总结,要不然就忘了...

多项式的表示

首先学习两种多项式的表达法

  1. 系数表达法
  2. 点值表达法

系数表达法

比如一个 \(n\) 次多项式 \(A(x)\) 他有 \(n + 1\) 项 于是 设每一项的系数为 \(a_i\) 则有 \(A(x) = \sum _{i = 0}^ n a_i x^ i\)

利用这种方法计算多项式乘法复杂度为 \(O(n^2)\) 即利用(第一个多项式中每个系数都需要与第二个多项式的每个系数相乘)

用这种方法,计算 \(A(x_0)\) 的值是非常好用的,我们使用霍纳法则即可。

点值表达法

我们知道如果给出 \(n + 1\) 个互不相同的点,他们确定一个 \(n\) 次多项式,这个结论是显然的

我们在下文中对于次数的规定稍微放开,规定 A 的次数界 (即最大次数 < 次数界)为 \(K_a\) 而同理,有 \(B_k\) 。通常的,把 \(A(x)*B(x)\) 的次数界规定为 \(K_a + K_b\)

对于许多多项式的操作,点值表达式十分方便,例如有一个点值表达式 A

\[\{\{x_0, y_0\}, \{x_1, y_1\}\cdots\{x_{n - 1}, y_{n-1}\}\} \]

还有一个点值表达式 B

\[\{\{x_0, y'_0\}, \{x_1, y_1'\}\cdots\{x_{n - 1}, y'_{n-1}\}\} \]

那么 \(A + B\) 的表达就是

\[\{\{x_0, y_0 + y'_0\}, \{x_1, y_1 + y'_1\}\cdots\{x_{n - 1}, y_{n-1} + y'_{n-1}\}\} \]

表达法END

复数浅讲

上过高中的请跳过,我是初中生,会讲锅的。

定义一个复数 \(u = a + bi\) 其中 \(i = \sqrt{-1}\) 。你其实可以按向量那么理解。

然后我们复数乘法显得巧妙

比如说一个

复数 \(b\) * 复数 \(c\) 可以理解为就是旋转,比如这个图

看到了吗,旋转,模长相乘。

我们还有一类特殊的向量——\(n\) 次本源单位根 \(\omega\)

-.-。我谔谔,表示什么,就是方程 \(x^n = i\) 的解

这就很显然了

有几个关于本源单位根的定理

  • 消去引理: \(\omega ^{dk}_{dn} = \omega ^k_n\)
  • 对于任意偶数 \(n > 0\)\(\omega ^{\frac {n}{2}}_n = \omega^1_2 = -1\)
  • 折半引理:如果 \(n\) 为偶数,则有 \((\omega^k_n)^2 = (\omega ^{k + n/2}_n)^2\) 从复数相乘的角度上看这是显然的
  • 求和引理:对于任意整数 \(n\ge1\) 且 不能被 \(n\) 整除的 \(k\),有 \(\sum _{j=0}^{n-1} (w^k_n )^j = 0\) 为啥,由于我们可以用等比数列求和公式来计算为 \(\sum ^{n-1}_0 (\omega_n^k)^j = \frac {(w^k_n)^n - 1}{w_n^k -1} = 0\)

DFT

离散小波变换 离散傅里叶变换

计算次数界为 \(n\) 的 的多项式 \(A(x) = \sum _0^{n-1} a_j x^j\) 我们不妨带入每个单位根得到 \(y_k = A(\omega_n^k) = \sum _{j = 0} ^ {n-1}a_j \omega _n ^{kj}\)

FFT

快速傅里叶变换、利用单位复数根的性质,我们可以在 \(O(n\log n)\) 的时间内计算 出 \(DFT_n\)

不妨假设都为 \(n\) 以次数界 , 且 \(n\) 为 2 的幂次。

不妨计算 \(A(\omega_n^k)\) 然后按下标的奇偶分类为 \(A_1A_2\)

\[A(\omega_n^k) = A_1(w_n^{2k}) + \omega_n^k A_2(\omega_n^{2k}) \]

我们通过观察 有

\[A(\omega_n^{k + n/2}) = A_1(w_n^{2k}) - \omega_n^k A_2(\omega_n^{2k}) \]

所以,我们在计算上一个时,可以 \(O(1)\) 得到下一个 于是就很明显了

于是明天再写 逆FFT变换

#include <bits/stdc++.h>

template <typename T>
inline T read()
{
	T x = 0;
	char ch = getchar();
	bool f = 0;
	while(ch < '0' || ch > '9')
	{
		f = (ch == '-');
		ch = getchar();
	}
	while(ch <= '9' && ch >= '0')
	{
		x = (x << 1) + (x << 3) + ch - '0';
		ch = getchar();
	}
	return f ? -x: x;
}

template <typename T>
void put(T x)
{
	if(x < 0) x = -x, putchar('-');
	if(x < 10)
	{
		putchar(x + 48);
		return ;
	}
	put(x / 10);
	putchar(x % 10 + 48);
	return;
}

#define rd read <int>
#define pt(x) put <int> (x), putchar(' ')

const int Maxn = 1e6 + 111;

int n, m, limit;

typedef double db;

namespace myspace
{

struct complex
{
	db x, y;
	complex(db xx, db yy){ x = xx, y = yy; }
	complex(){ x = 0, y = 0; }
	complex operator * (const complex &a) { return complex(x * a.x - y * a.y, x * a.y + y * a.x); }
	complex operator + (const complex &a) { return complex(x + a.x , y + a.y); }
	complex operator - (const complex &a) { return complex(x - a.x , y - a.y); }
	db abs() { return sqrt(x * x + y * y); }
};

};

myspace::complex a[Maxn], b[Maxn];

const db Pi = 3.1415926535897;

void fft(int limit, myspace::complex *a, int type)
{
	if(limit == 1) return;
	myspace::complex a1[(limit >> 1)], a2[(limit >> 1)];
  // 写的很巧妙,注意一定是 < limit
	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);
	myspace::complex w(1, 0), wn(cos(2.0 * Pi / limit), type * sin(2.0 * Pi / limit));
  // 记录本源单位根,同时计算a
	for(int i = 0; i < limit >> 1; ++i, w = w * wn)
    a[i] = a1[i] + w * a2[i], a[i + (limit >> 1)] = a1[i] - w * a2[i];
	return;
}

int main()
{
#ifdef _DEBUG
	freopen("in.txt", "r", stdin);
	// freopen("out.txt", "w", stdout);
#endif
	n = rd();
	m = rd();
	for(int i = 0; i <= n; ++i) a[i].x = rd();
	for(int i = 0; i <= m; ++i) b[i].x = rd();
	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) pt((int) (a[i].x / limit + 0.5));
	putchar('\n');
	return 0;
}

猜你喜欢

转载自www.cnblogs.com/zhltao/p/12633158.html
FFT