复数

YZY Lv3

一、引言:超越一维的数字

我们对数字的认知,始于在数轴上左右移动的整数与实数。这是一个一维的世界。然而,在广阔的算法问题域中,我们频繁面对二维乃至更高维度的结构:平面上的点、向量的旋转、周期性的变换。将这些二维问题强行压缩到一维的实数轴上处理,往往会使问题变得复杂、表达式冗长、几何直觉尽失。

复数,正是为了优雅地描述二维世界而生的数学工具。一个复数 内含了两个独立的信息(实部 和虚部 ),天然地对应了平面上的一个点或一个向量 。这种对应关系并非简单的“记号游戏”,它深刻地改变了我们对运算的理解。复数的加法完美对应了向量的平移,而其乘法,则蕴含了二维空间中至关重要的操作——旋转与伸缩

这使得复数成为连接代数与几何的桥梁。许多看似棘手的几何问题,在复数域内可以转化为简洁的代数运算。更令人惊叹的是,复数域中的一个特殊概念——单位根,其完美的对称性与周期性,构成了计算机科学中最强大的算法之一“快速傅里叶变换”(FFT)的基石,彻底改变了我们处理多项式与卷积问题的方式。

本文将带领你系统地探索复数在算法设计中的应用。我们将从复数的基础表示法出发,深入其几何内涵,剖析其如何简化二维向量运算,并最终为你揭示单位根的奥秘,为理解FFT铺平道路。

二、复数的表示与基本运算

掌握任何工具的第一步,是熟悉它的形态与基本操作。复数有三种相互关联的表示方式:代数表示、几何表示和极坐标表示。

2.1 代数表示

这是复数最基本的定义。一个复数 被写为:

其中 是实数,分别称为 实部(Real Part)和虚部(Imaginary Part),记作 虚数单位,其定义为

所有复数的集合构成了复数域,用 表示。实数域 可以看作是复数域中虚部为零的子集。

复数的加减法遵循实部与虚部分别相加减的原则,这与向量的加减法完全一致:

复数的乘法按多项式乘法展开,并利用 进行化简:

复数的除法通过一个技巧——乘以分母的共轭(Conjugate)——来实现分母的实数化。复数 的共轭为

2.2 几何表示:复平面

代数表示揭示了复数的“构成”,而几何表示则赋予其“生命”。我们可以用一个二维平面来表示复数,这个平面被称为复平面阿尔冈图(Argand Diagram)。平面上,水平的 轴代表实部,称为实轴;垂直的 轴代表虚部,称为虚轴

复数 就唯一地对应了复平面上的一个点 ,或者说,一个从原点 指向点 的向量。

在这种视角下:

  • (Modulus):复数 对应的向量的长度,记作 。根据勾股定理, 。模表示了点到原点的距离。
  • 幅角(Argument):向量与实轴正方向的夹角,记作 。它定义在 区间内。幅角表示了点的方向。通过反正切函数计算,atan2(y, x) 是一个在编程中常用的函数,它比 atan(y/x) 更健壮,能正确处理 的情况并返回覆盖 的完整角度。

2.3 极坐标表示

有了模和幅角,我们可以用极坐标 来描述一个复数。根据三角函数定义,我们有 。代入代数形式 ,得到:

这就是复数的极坐标表示。它将一个复数分解为大小 () 和方向 () 两个独立的部分,为后续理解乘法几何意义打下基础。

2.4 欧拉公式:最美的数学公式

瑞士数学家欧拉发现了一个深刻而优美的联系,它将指数函数、三角函数和虚数单位 完美地结合在一起:

这个公式被称为欧拉公式。它断言,一个纯虚指数的指数函数,实际上是单位圆上的一个点,其实部是 ,虚部是

欧拉公式的推导超出了本文范畴(通常使用泰勒级数展开证明),但我们必须理解它的核心意义:它将旋转运动(由三角函数描述)转化为了指数运算

将欧拉公式代入极坐标表示,我们得到复数最紧凑、最强大的形式:

这个形式清晰地揭示了复数的本质: 是大小(伸缩因子), 是方向(旋转因子)。

2.5 C++中的<complex>

为了让使用者能够将精力集中在算法逻辑而非复数的底层实现上,C++ 标准库在 <complex> 头文件中提供了一个强大且易用的模板类 std::complex。熟练掌握它,是高效实现复数相关算法的前提。我们通常使用 std::complex<double>,并为其定义一个短别名 cd

1
2
#include <complex>
using cd = std::complex<double>;

下面将系统地介绍其常用功能。

2.5.1 构造与初始化

std::complex 提供了多种灵活的构造方式。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
// 1. 默认构造函数,值为 (0, 0)
cd z1; // z1 = 0 + 0i

// 2. 使用实部和虚部构造
cd z2(3.0, 4.0); // z2 = 3 + 4i

// 3. 仅使用实部构造 (虚部默认为0)
cd z3(5.0); // z3 = 5 + 0i

// 4. 列表初始化 (C++11)
cd z4 = {2.0, -1.0}; // z4 = 2 - i

// 5. 从另一个复数拷贝构造
cd z5 = z2; // z5 = 3 + 4i

// 6. 从实数隐式转换
cd z6 = 7.0; // z6 = 7 + 0i

2.5.2 成员函数与访问器

可以直接访问或修改复数的实部和虚部。

  • T real() const; / void real(T);
  • T imag() const; / void imag(T);
1
2
3
4
5
6
7
8
9
10
11
cd z(3, 4);

// 访问 (const版本)
double r = z.real(); // r = 3.0
double i = z.imag(); // i = 4.0

// 修改
z.real(5); // z 变为 5 + 4i
z.imag(-2); // z 变为 5 - 2i

cout << z; // C++重载了 << 操作符, 输出 (5,-2)

2.5.3 运算符重载

std::complex 重载了所有常见的算术运算符,使得复数运算如同原生类型一样自然。

赋值运算符: operator=, operator+=, operator-=, operator*=, operator/=
一元运算符: operator+ (正号), operator- (负号)
二元运算符: operator+, operator-, operator*, operator/

这些运算符支持 complexcomplexcomplexT (实数)、Tcomplex 之间的运算。

1
2
3
4
5
cd z1(1, 1), z2(3, 4);
z1 += z2; // z1 变为 (4, 5)
cd sum = z1 + z2; // (4,5) + (3,4) = (7, 9)
cd prod = z1 * 2.0; // (4,5) * 2.0 = (8, 10)
cd quot = z1 / z2; // (4,5) / (3,4)

2.5.4 非成员函数 (核心功能)

除了基本的运算符,<complex> 还提供了一系列以自由函数形式存在的强大工具。

2.5.4.1 几何相关函数
  • T abs(const complex<T>& z): 计算复数 z (modulus)

    1
    2
    cd z(3, 4);
    double m = abs(z); // m = 5.0
  • T arg(const complex<T>& z): 计算复数 z幅角 (argument) ,单位是弧度,范围

    1
    2
    cd z(1, 1);
    double a = arg(z); // a = PI/4 (约 0.785)
  • T norm(const complex<T>& z): 计算复数 z模的平方 。在比较模的大小时,使用 normabs 更高效,因为它避免了开方运算。

    1
    2
    cd z(3, 4);
    double n = norm(z); // n = 25.0
  • complex<T> conj(const complex<T>& z): 计算复数 z共轭 (conjugate)

    1
    2
    cd z(3, 4);
    cd c = conj(z); // c = 3 - 4i
  • complex<T> polar(const T& r, const T& theta = 0): 从极坐标 构造复数 。这是创建旋转因子或单位根时极其有用的函数。

    1
    2
    const double PI = acos(-1.0);
    cd z = polar(5.0, PI / 2); // 模为5, 幅角90度, z 约等于 0 + 5i
  • complex<T> proj(const complex<T>& z): 计算复数 z黎曼球面上的投影。在大多数竞赛场景下不常用。

2.5.4.2 指数、对数、幂、开方函数
  • complex<T> exp(const complex<T>& z): 计算
  • complex<T> log(const complex<T>& z): 计算自然对数(主值)
  • complex<T> pow(const complex<T>& base, ...): 计算幂
  • complex<T> sqrt(const complex<T>& z): 计算平方根(主值)。
1
2
3
4
5
6
7
8
const double PI = acos(-1.0);
cd i(0, 1);
// 验证欧拉公式 e^(i*pi) = -1
cout << exp(cd(0, PI)) << endl; // 输出约 (-1,0)
// 计算 i^2
cout << pow(i, 2.0) << endl; // 输出 (-1,0)
// 计算 sqrt(-1)
cout << sqrt(cd(-1, 0)) << endl; // 输出 (0,1)
2.5.4.3 三角函数与双曲函数

标准库也为复数定义了全套的三角函数(sin, cos, tan)、反三角函数(asin, acos, atan)和双曲函数。它们在一些高等数学或物理相关的题目中可能出现,但在主流的FFT体系中不常用。

2.5.5 使用建议

  • 别名: 始终使用 using cd = complex<double>; 来简化代码。
  • 几何意义: 牢记 polar, abs, arg 的几何意义,它们是连接代数与几何的桥梁。在处理旋转、单位根等问题时,polar 是你的首选工具。
  • 效率: 在比较模长时,优先使用 norm 避免不必要的 sqrt 计算。
  • 统一性: std::complex 的设计与原生数值类型高度统一,你可以放心地将它用于模板元编程或泛型算法中。

通过熟练运用 std::complex,我们可以将复杂的复数运算逻辑委托给标准库,从而专注于算法本身的设计与优化。

三、复数乘法:旋转与伸缩的艺术

复数加减法与向量加减法别无二致,相对直观。但复数乘法的几何意义,才是其在算法中大放异彩的关键。

3.1 代数视角

我们回顾一下复数乘法的代数定义:

这个公式是正确的,但它看起来杂乱无章,缺乏任何直观的几何解释。我们很难从这个式子中看出输入 与输出 之间的几何关系。

3.2 几何视角:模长相乘,幅角相加

现在,让我们切换到极坐标表示 ,看看乘法会发生什么。
。它们的乘积是:

这个结果优雅得令人屏息。它清晰地揭示了复数乘法的几何本质:

  1. 模相乘:新复数的模是原来两个复数模的乘积,即
  2. 幅角相加:新复数的幅角是原来两个复数幅角的和,即

结论:将一个复数 乘以另一个复数 ,在几何上等价于对 所对应的向量进行两个独立的操作:

  1. 伸缩:将其长度(模)乘以
  2. 旋转:将其绕原点逆时针旋转 的角度。

这个深刻的洞见是复数在计算几何中应用的核心。原本需要通过旋转矩阵、三角函数变换才能完成的旋转操作,现在仅仅是一个复数乘法。

3.3 应用:向量旋转

这是复数乘法最直接、最经典的应用。

题意概括:
给定平面上的一个点 和一个旋转角度 (以弧度为单位)。求 绕原点逆时针旋转 角度后得到的新点 的坐标。

解题思路:

  1. 将点 视为复数
  2. 构造一个用于旋转的复数“操作符” 。我们希望这个操作符只旋转而不伸缩,所以它的模必须为1。根据欧拉公式,这个操作符就是
  3. 执行复数乘法
  4. 新复数 的实部和虚部就是新点 的坐标。

C++ 代码示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using cd = complex<double>;

// PI 的高精度定义
const double PI = acos(-1.0);

// p: 待旋转的点 (x, y)
// a: 旋转角度 (弧度)
void solve() {
double x, y, a;
cin >> x >> y >> a;

// 1. 将点 (x, y) 表示为复数 p
cd p(x, y);

// 2. 创建旋转操作符 o
// 模为 1, 幅角为 a
// polar(r, theta) 函数可以直接从极坐标创建复数
cd o = polar(1.0, a);

// 3. 执行乘法完成旋转与伸缩
cd res = p * o;

// 4. 输出新坐标
cout << fixed << setprecision(10) << res.real() << " " << res.imag() << '\n';
}

int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);

solve();

return 0;
}

代码解释:

  • cd p(x, y); 将输入的笛卡尔坐标点 直接构造成复数
  • cd o = polar(1.0, a); 是本解法的核心。polar(r, theta) 函数根据模 和幅角 theta 生成复数 。我们想要纯旋转,所以模为1.0,幅角为给定的旋转角 a。这比手动计算 cos(a)sin(a) 来创建 cd(cos(a), sin(a)) 更简洁、意图更明确。
  • cd res = p * o; 一行代码完成了过去需要矩阵乘法或繁琐三角函数推导才能完成的旋转操作。这就是复数的威力。
  • res.real()res.imag() 分别提取结果的 坐标。

复杂度分析:

  • 时间复杂度: 。所有操作都是常数时间的浮点数运算。
  • 空间复杂度: 。只使用了几个变量来存储复数。

相比之下,传统的旋转公式是:


如果我们展开复数乘法 ,会得到 ,其结果与旋转公式完全吻合。这再次证明了复数乘法内在地执行了旋转操作。

四、单位根:周期性与对称性的化身

在复数域中,有一类特殊的数,它们是代数与几何完美结合的典范,并且是FFT算法的灵魂。它们就是单位根

4.1 定义

在复数域中,方程 的所有解被称为 次单位根(n-th roots of unity)。其中 是一个正整数。

换言之, 次单位根是那些自身相乘 次后恰好等于 1 的复数。

4.2 求解与几何意义

如何找到这些解?让我们再次请出强大的极坐标表示
方程 可以写作:


我们将右边的 1 也写成极坐标形式。1 的模是 1,幅角是 。但由于幅角的周期性(转一圈 回到原位),它的幅角也可以是 。所以,,其中 是任意整数。

于是方程变为:

要使等式成立,两边的模和幅角必须分别相等。

  • 模相等: 。因为 是非负实数,所以 。这说明所有单位根都在复平面的单位圆上
  • 幅角相等: 。因此

我们只需取 就能得到 个不同的解。因为当 时,,这与 时的 对应同一个复数。

所以, 次单位根共有 个,它们是:

我们通常用 (Omega) 来表示第 次单位根。 被称为 次单位根

几何意义:
个单位根在复平面上构成了什么图形?

  • 它们的模都为 1,所以它们均匀地分布在以原点为中心的单位圆上。
  • 它们的幅角分别为 。相邻两个单位根的幅角差是恒定的

这意味着,** 次单位根在复平面上构成了一个内接于单位圆的正 边形**,其中一个顶点在 处(即 )。

下图展示了 时的情况,8个8次单位根构成了一个正八边形。

4.3 单位根的性质

单位根的优美远不止于其几何形态,它们拥有的代数性质是FFT算法能够成立的理论基石。假设 是一个偶数。

  1. 基本性质: 。这由指数运算律 直接得到。

  2. 周期性 (Periodicity):
    这说明单位根的指数对模 具有周期性。

  3. 对称性 (Symmetry):
    证明:

    根据欧拉公式,
    所以,
    几何解释: 是将 在单位圆上再旋转半圈( 弧度),恰好到达其关于原点的对称点。这个性质是FFT中“蝶形运算”的直接来源。

  4. 折半引理 (Halving Lemma):
    证明:

    意义: 这个引理是分治思想的体现。它告诉我们,一个规模为 的问题中的单位根的平方,恰好是规模为 的子问题中使用的单位根。这是FFT能够递归分解的基础。

  5. 求和引理 (Summation Lemma):

    证明:

    • 时,。求和式变为
    • 时,。求和式是一个等比数列。根据求和公式

      意义: 这个性质在傅里叶逆变换(IFFT)中扮演了关键角色,用于从变换后的结果还原出原始系数。它本质上体现了单位根向量之间的“正交性”。

4.4 应用:计算单位根与验证性质

我们可以用C++的 std::complex 来方便地计算和验证这些性质。

题意概括:
给定一个偶数 和一个整数 ()。计算并输出 , , ,以验证对称性和折半引理。

C++ 代码示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using cd = complex<double>;

const double PI = acos(-1.0);

// 辅助函数,用于打印复数
void print_cd(const string& s, cd z) {
cout << s << fixed << setprecision(6) << z.real()
<< " " << z.imag() << "i\n";
}

void solve() {
int n, k;
cin >> n >> k;

if (n % 2 != 0 || k < 0 || k >= n / 2) {
cout << "Invalid input\n";
return;
}

// 计算 omega_n^k
cd wnk = polar(1.0, 2 * PI * k / n);
print_cd("w_n^k : ", wnk);

// 验证对称性
cout << "\n--- Verifying Symmetry Property ---\n";
cd wnk_half = polar(1.0, 2 * PI * (k + n / 2.0) / n);
print_cd("w_n^(k+n/2) : ", wnk_half);
print_cd("-w_n^k : ", -wnk);

// 验证折半引理
cout << "\n--- Verifying Halving Lemma ---\n";
cd wnk_sq = wnk * wnk;
print_cd("(w_n^k)^2 : ", wnk_sq);
cd wn2k = polar(1.0, 2 * PI * k / (n / 2.0));
print_cd("w_{n/2}^k : ", wn2k);
}

int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);

solve();

return 0;
}

输入样例:
8 1

输出样例:

1
2
3
4
5
6
7
8
9
w_n^k         : 0.707107 0.707107i

--- Verifying Symmetry Property ---
w_n^(k+n/2) : -0.707107 -0.707107i
-w_n^k : -0.707107 -0.707107i

--- Verifying Halving Lemma ---
(w_n^k)^2 : 0.000000 1.000000i
w_{n/2}^k : 0.000000 1.000000i

代码解释:

  • 我们再次使用 polar(1.0, angle) 来生成单位根,这非常直观。
  • 对于对称性验证,我们分别计算 。可以看到后者等于前者的相反数。
  • 对于折半引理验证,我们计算 。可以看到两者是相等的(由于浮点精度问题,可能存在极小的误差,但数值上是吻合的)。
  • 这个例子具体地展示了单位根的性质,将抽象的公式与可计算的数值联系起来,加深了理解。

复杂度分析:

  • 时间复杂度:
  • 空间复杂度:

五、多项式与卷积

我们已经建立了对复数和单位根的深刻理解。现在,是时候将这些工具应用到一个核心的计算问题上:多项式乘法。这不仅是一个常见问题,其背后所蕴含的卷积(Convolution)思想,在信号处理、图像处理、概率论和字符串匹配等众多领域都有着广泛的应用。

5.1 多项式的两种表示法

一个 次多项式 可以用其系数来表示,这被称为系数表示法(Coefficient Representation)。

这个多项式由一个系数向量 唯一确定。

如果我们取 个不同的点 ,计算出多项式在这些点上的值 ,那么我们就得到了一组点值对 。这被称为多项式的点值表示法(Point-value Representation)。一个重要的定理是,这 个点值对同样可以唯一确定一个 次的多项式。

让我们比较一下在这两种表示法下,多项式运算的复杂度:

  • 加法

    • 系数表示:将对应系数相加,。复杂度为
    • 点值表示:若 在相同的点集 上求值,则 。复杂度也为
  • 乘法。假设 都是 次,则 是一个 次的多项式。

    • 系数表示:根据定义,。其第 项的系数 是由所有 项相加而成:

      计算所有系数需要 的时间。
    • 点值表示:若 在相同的 个(或更多)点上求值,则 。计算所有点值对只需要 的时间。

这个对比揭示了一个惊人的事实:在点值表示下,多项式乘法异常简单。这启发我们设计一个全新的算法流程:

  1. **求值 (Evaluation)**:将系数表示的 转换为点值表示。
  2. **点积 (Pointwise Product)**:用 时间计算出 的点值表示。
  3. **插值 (Interpolation)**:将 的点值表示转换回系数表示。

这个流程的瓶颈在于第一步和第三步。如果我们能快速地完成求值和插值,就能实现快速的多项式乘法。快速傅里叶变换(FFT)正是为此而生。

5.2 卷积

我们再看一眼系数表示法下的乘法公式:

这个形式在数学和工程中极为常见,它被称为两个序列 离散卷积(Discrete Convolution)。记作

因此,多项式乘法本质上是其系数向量的卷积。FFT 作为一个强大的工具,其核心用途就是计算离散卷积。

5.3 宏伟蓝图:FFT 加速卷积

FFT 的角色就是高效地完成“求值”和“插值”。具体来说:

  • 求值过程被称为**离散傅里叶变换 (DFT)**。FFT 是一种在 时间内计算 DFT 的算法。它通过巧妙地选择 次单位根 作为求值点 来实现加速。
  • 插值过程被称为**逆离散傅里叶变换 (IDFT)**。IDFT 也可以通过稍作修改的 FFT 算法在 时间内完成。

因此,利用 FFT,我们可以将多项式乘法(卷积)的复杂度从 降至

六、快速傅里叶变换

现在,我们正式进入 FFT 的核心地带。我们将揭示,前文所述的单位根性质是如何被分治思想所利用,从而创造出这个高效的算法。

6.1 离散傅里叶变换 (DFT)

给定一个 次多项式 (或者说,一个长度为 的序列 ),其 DFT 就是将 次单位根 () 代入 得到的 个值的序列


这个过程可以看作是一个矩阵向量乘法:
$$



$$
这个矩阵被称为范德蒙德矩阵。直接计算这个矩阵乘法需要 时间。FFT 的目标就是优化这个计算。

6.2 Cooley-Tukey 算法:分治的魔术

FFT 最经典的算法是 Cooley-Tukey 算法,它基于分治思想。为简化讨论,我们假设 是 2 的幂。如果不是,我们可以将多项式的高次项系数补零,使其长度达到最近的 2 的幂。

我们的目标是计算
第一步:奇偶分离
我们将多项式 按系数的下标是奇数还是偶数,拆分成两个新的多项式:

这两个多项式的次数都是
通过简单的代数变形,我们可以将 表示为:

这个恒等式是整个 FFT 算法的基石。

第二步:分治计算
我们想要求 的值。利用上面的恒等式:

根据折半引理 ,上式变为:

现在考虑计算 DFT 的后半部分,即 ():

我们再次使用折半引理:
并使用对称性
代入得到:

第三步:合并结果 (蝶形变换)
我们来比较一下 的表达式,其中 :


这揭示了一个惊人的结构:

  • 我们只需要递归地计算出 这些点上的值。
  • 那么, 可以通过一次乘法、一次加法和一次减法,由 计算出来:

    • 这个计算过程被称为蝶形变换(Butterfly Transform),因为其数据流图形状酷似蝴蝶。

复杂度分析
我们把一个规模为 的 DFT 问题,转化为了两个规模为 的 DFT 子问题,以及 的合并操作(对每个 from to 做一次蝶形变换)。
递归关系式为 。根据主定理,这个递归式解为

6.3 逆离散傅里叶变换 (IDFT)

现在我们有了从系数到点值的快速算法。如何从点值 快速地插值回系数 呢?
IDFT 的公式可以被证明为:

这个公式与 DFT 的公式 形式上惊人地相似。唯一的区别是:

  1. 单位根的指数从 变成了
  2. 最后的结果需要除以

这意味着,我们可以复用 FFT 的代码来计算 IDFT:

  1. 将 FFT 算法中使用的单位根 替换为它的共轭 。或者等价地,将计算单位根时使用的角度 替换为
  2. 将 FFT 的输入从系数向量 换成点值向量
  3. 对计算出的结果,每个元素都除以

因此,IDFT 也可以在 时间内完成。

七、FFT 的实现:从递归到迭代

虽然递归实现能够清晰地反映分治思想,但在实际编码中,由于函数调用开销和栈深度问题,通常采用更高效的迭代实现

7.1 递归 FFT (仅为理解)

一个简单的递归 FFT 实现框架如下,这有助于我们理解其逻辑,但不推荐在需要高性能的场合使用。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
// 伪代码/逻辑示意
void fft_recursive(vector<cd>& a) {
int n = a.size();
if (n <= 1) return;

vector<cd> a_even(n/2), a_odd(n/2);
for (int i = 0; i < n/2; ++i) {
a_even[i] = a[2*i];
a_odd[i] = a[2*i+1];
}

fft_recursive(a_even);
fft_recursive(a_odd);

double ang = 2 * PI / n;
cd w(1), wn(cos(ang), sin(ang));

for (int k = 0; k < n/2; ++k) {
cd t = w * a_odd[k];
a[k] = a_even[k] + t;
a[k+n/2] = a_even[k] - t;
w *= wn;
}
}

7.2 蝴蝶变换与位逆序

递归的奇偶分离过程,如果我们追踪到底层,会发现一个有趣的规律。以 为例:

  • 第0层:[0,1,2,3,4,5,6,7]
  • 第1层(奇偶分离):[0,2,4,6][1,3,5,7]
  • 第2层(再次分离):[0,4], [2,6][1,5], [3,7]
  • 第3层(最终分离):[0], [4], [2], [6], [1], [5], [3], [7]

最终的序列 [0,4,2,6,1,5,3,7] 是什么?我们来看它们的二进制表示:

  • 000 -> 000 (0)
  • 001 -> 100 (4)
  • 010 -> 010 (2)
  • 011 -> 110 (6)
  • 100 -> 001 (1)
  • 101 -> 101 (5)
  • 110 -> 011 (3)
  • 111 -> 111 (7)
    可以发现,最终位置的下标,恰好是原始位置下标的二进制位翻转(Bit-Reversal)后的结果。

这个发现让我们可以设计一个迭代的 FFT 算法:

  1. 位逆序置换:首先,将输入的系数向量按照位逆序规则重新排列。这样,最底层的蝶形变换所需要的元素就被放在了一起。
  2. 自底向上计算:然后,我们从底层开始,逐层向上进行蝶形变换。
    • 第1层:对长度为 2 的小块进行蝶形变换。
    • 第2层:对长度为 4 的小块进行蝶形变换。
    • 层:对整个长度为 的数组进行一次大的蝶形变换。

这样就避免了递归,将算法转化为一系列循环,大大提高了效率。

7.3 迭代 FFT 实现

题意概括:
实现一个函数,输入一个复数向量 a,和一个代表变换方向的整数 op (1 表示 FFT, -1 表示 IFFT),原地修改向量 a 为其傅里叶变换或逆变换的结果。

C++ 代码示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using cd = complex<double>;

const double PI = acos(-1.0);

// rev: 存储位逆序置换后的下标
// a: 进行FFT的复数向量
// op: 操作类型, 1 for FFT, -1 for IFFT
void fft(vector<cd>& a, int op) {
int n = a.size();
vector<int> rev(n);

// 1. 位逆序置换预处理
for (int i = 0; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]); // 保证每个对只交换一次
}

// 2. 自底向上进行蝶形变换
for (int len = 2; len <= n; len <<= 1) { // len: 当前合并的块长度
double ang = op * 2 * PI / len;
cd wn = polar(1.0, ang); // len次单位根
for (int i = 0; i < n; i += len) { // i: 每个块的起始位置
cd w(1);
for (int k = 0; k < len / 2; ++k) { // k: 块内蝶形变换
cd x = a[i + k];
cd y = w * a[i + k + len / 2];
a[i + k] = x + y;
a[i + k + len / 2] = x - y;
w *= wn; // 更新旋转因子
}
}
}

// 3. IFFT需要额外除以n
if (op == -1) {
for (int i = 0; i < n; ++i) {
a[i] /= n;
}
}
}

// 示例主函数
void solve_poly_mul() {
// 假设A(x) = 1 + 2x, B(x) = 3 + 4x
// C(x) = (1+2x)(3+4x) = 3 + 10x + 8x^2
vector<cd> a = {1, 2}, b = {3, 4};

int n = 1;
// C(x)的次数是 (a.size()-1) + (b.size()-1) = 1+1=2, 需要3个点
// n 必须是2的幂且 >= 3, 所以n=4
while (n < a.size() + b.size() - 1) n <<= 1;

a.resize(n);
b.resize(n);

fft(a, 1);
fft(b, 1);

vector<cd> c(n);
for (int i = 0; i < n; ++i) c[i] = a[i] * b[i];

fft(c, -1);

cout << "Result coefficients:\n";
for (int i = 0; i < n; ++i) {
cout << round(c[i].real()) << " ";
}
cout << "\n"; // 输出应为 3 10 8 0
}

int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);

solve_poly_mul();

return 0;
}

代码解释:

  • 位逆序置换: rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0); 这是一个递推计算位逆序的巧妙方法。rev[i]rev[i/2] 计算得到。rev[i/2]i/2 的翻转,右移一位后得到 i 除了最高位之外的翻转。然后根据 i 的最低位(i&1)来确定 i 翻转后的最高位。
  • 主循环: 外层循环 for (int len ...) 控制合并块的长度,从2, 4, 8, … 一直到 n。中层循环 for (int i ...) 遍历所有块。内层循环 for (int k ...) 在一个块内执行蝶形变换。
  • 旋转因子: wn 是当前块长度对应的主单位根w 在块内迭代,从 1 开始,每次乘以 wn
  • op 参数: 通过 op 巧妙地统一了 FFT 和 IFFT。op=-1 时,角度为负,相当于取共轭,符合 IFFT 的要求。

复杂度分析:

  • 时间复杂度: 。位逆序置换 。主循环有 层,每层内部的两个循环总共执行 次操作。
  • 空间复杂度: 。需要一个 rev 数组和存储复数的向量。算法是原地进行的,没有额外的递归栈开销。

八、应用:高精度乘法

FFT最直接、最经典的应用就是加速大整数乘法。

题意概括:
给定两个大整数 ,它们的位数可能达到 或更高。计算它们的乘积

解题思路:

  1. 转换为多项式: 将大整数看作一个多项式在 处的求值。例如,整数 123 可以看作多项式 。我们将大整数的每一位数字(从低位到高位)作为多项式的系数。
  2. 补零: 设两个数位数分别为 。乘积的位数最多为 。我们需要选择一个 2 的幂 ,将两个多项式的系数向量都补零到长度
  3. FFT: 对两个系数向量进行 FFT。
  4. 点积: 将变换后的结果逐点相乘。
  5. IFFT: 对点积结果进行 IFFT,得到乘积多项式的系数向量。
  6. 处理进位: IFFT 的结果是浮点数,需要四舍五入为整数。由于 ,这些系数可能会大于9,所以需要从低位到高位处理进位,模拟手工乘法的过程。

C++ 代码示例:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <bits/stdc++.h>

using namespace std;
using ll = long long;
using cd = complex<double>;

const double PI = acos(-1.0);

// 复用上一节的 fft 函数
void fft(vector<cd>& a, int op) {
int n = a.size();
vector<int> rev(n);
for (int i = 0; i < n; ++i) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) ? (n >> 1) : 0);
if (i < rev[i]) swap(a[i], a[rev[i]]);
}
for (int len = 2; len <= n; len <<= 1) {
double ang = op * 2 * PI / len;
cd wn = polar(1.0, ang);
for (int i = 0; i < n; i += len) {
cd w(1);
for (int k = 0; k < len / 2; ++k) {
cd x = a[i + k], y = w * a[i + k + len / 2];
a[i + k] = x + y;
a[i + k + len / 2] = x - y;
w *= wn;
}
}
}
if (op == -1) {
for (auto& x : a) x /= n;
}
}

void solve() {
string s1, s2;
cin >> s1 >> s2;

int n1 = s1.length(), n2 = s2.length();
int lim = 1;
while (lim < n1 + n2) lim <<= 1;

vector<cd> a(lim), b(lim);
for (int i = 0; i < n1; ++i) a[i].real(s1[n1 - 1 - i] - '0');
for (int i = 0; i < n2; ++i) b[i].real(s2[n2 - 1 - i] - '0');

fft(a, 1);
fft(b, 1);

for (int i = 0; i < lim; ++i) a[i] *= b[i];

fft(a, -1);

vector<int> res(lim);
int carry = 0;
for (int i = 0; i < lim; ++i) {
// 加0.5是为了四舍五入
int val = (int)(a[i].real() + 0.5) + carry;
res[i] = val % 10;
carry = val / 10;
}

// 处理最高位的进位
int top = lim - 1;
while (carry > 0) {
if(top >= res.size()) res.push_back(0);
res[top] += carry % 10;
carry /= 10;
top++;
}

// 找到最高非零位并输出
int pos = res.size() - 1;
while (pos > 0 && res[pos] == 0) pos--;

for (int i = pos; i >= 0; --i) {
cout << res[i];
}
cout << '\n';
}

int main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
solve();
return 0;
}

复杂度分析:

  • 时间复杂度: ,其中 是两个输入大数位数之和。主要开销在三次 FFT/IFFT 调用上。
  • 空间复杂度: ,用于存储系数向量和 FFT 计算。

九、总结与展望

本文带领我们完成了一次从基础到前沿的复数之旅。我们始于复数作为二维向量的几何直觉,领略了其乘法如何优雅地实现旋转与伸缩。随后,我们深入探索了单位根——这一类镶嵌在单位圆上的特殊复数,揭示了其背后深刻的周期性、对称性和分治性质。

这些看似抽象的数学性质,最终在快速傅里叶变换(FFT)中融为一体,构成了一个能够以 复杂度解决多项式乘法和卷积问题的强大算法。我们不仅推导了其分治核心思想,还一步步实现了工业级的迭代版本,并将其应用于高精度乘法这一经典问题上。

复数与 FFT 的世界远未结束:

  • 数论变换 (NTT): 在处理需要精确整数运算(例如,对大质数取模)的卷积问题时,FFT 的浮点数精度会成为瓶颈。NTT 使用模意义下的“原根”替代单位根,在有限域中实现了与 FFT 相同的代数结构,从而避免了精度误差。
  • 快速沃尔什-阿达玛变换 (FWHT): 用于解决位运算(AND, OR, XOR)卷积问题。
  • 更多应用: FFT 的思想和卷积的应用渗透在字符串匹配(通配符匹配)、动态规划优化、生成函数等诸多高级算法之中。

希望本文为你打开了一扇新的大门。理解复数,不仅仅是掌握一个数学工具或一个算法,更是学会一种全新的视角——通过变换来简化问题。从系数表示到点值表示的切换,正是这一思想的绝佳体现。这种思维模式,将在你攀登算法高峰的道路上,提供源源不断的动力与启发。

  • Title: 复数
  • Author: YZY
  • Created at : 2025-07-05 01:52:30
  • Updated at : 2025-07-05 01:52:30
  • Link: https://blog.dtoj.team/2025/07/05/复数/
  • License: 三金家的Y同学