每日算法 — 使用java实现FFT：快速傅里叶变换与多项式乘法加速

算法概述

快速傅里叶变换（Fast Fourier Transform，FFT）是离散傅里叶变换（DFT）的高效算法，由Cooley和Tukey于1965年提出。FFT将DFT的时间复杂度从O(n²)降至O(n log n)，是20世纪最重要的数值算法之一。

在算法竞赛和工程实践中，FFT的核心价值在于多项式乘法的加速：两个n次多项式相乘，朴素算法需要O(n²)时间，而借助FFT可在O(n log n)内完成。这一特性使其成为大整数乘法、信号处理、卷积计算和密码学中的关键工具。

核心原理

多项式表示法

一个n-1次多项式有两种等价表示：
– 系数表示：A(x) = a_0 + a_1x + a_2x² + … + a_{n-1}x^{n-1}
– 点值表示：n个不同点(x_0, A(x_0)), …, (x_{n-1}, A(x_{n-1}))

卷积定理

两个多项式相乘等价于系数的离散卷积。在频域中，卷积变为逐点乘法：
– Coeff → DFT → 逐点乘 → IDFT → Coeff

单位复根与分治

取n次单位复根ω_n = e^{2πi/n}，利用其周期性：
– ω_n^k = ω_{n/2}^{k/2}（当k为偶数）
– ω_n^{k+n/2} = -ω_n^k

将多项式按奇偶项拆分，递归求解两个半规模子问题，然后O(n)合并。

Java完整实现

/**
 * 快速傅里叶变换(FFT)实现
 * 用于多项式乘法和卷积计算的O(n log n)算法
 */
public class FFT {

    /**
     * 复数类，支持FFT所需的基本运算
     */
    public static class Complex {
        public double real;
        public double imag;

        public Complex(double real, double imag) {
            this.real = real;
            this.imag = imag;
        }

        public Complex add(Complex other) {
            return new Complex(this.real + other.real, this.imag + other.imag);
        }

        public Complex subtract(Complex other) {
            return new Complex(this.real - other.real, this.imag - other.imag);
        }

        public Complex multiply(Complex other) {
            return new Complex(
                this.real * other.real - this.imag * other.imag,
                this.real * other.imag + this.imag * other.real
            );
        }

        public Complex divide(double scalar) {
            return new Complex(this.real / scalar, this.imag / scalar);
        }

        @Override
        public String toString() {
            return String.format("%.4f %+.4fi", real, imag);
        }
    }

    /**
     * 迭代版FFT（Cooley-Tukey算法，蝶形运算）
     * @param a 输入复数数组，长度必须是2的幂
     * @param invert true表示逆变换（IDFT），false表示正变换（DFT）
     */
    public static void fft(Complex[] a, boolean invert) {
        int n = a.length;

        // 位逆序置换：将数组按二进制位反转重新排列
        for (int i = 1, j = 0; i < n; i++) {
            int bit = n >> 1;
            for (; j >= bit; bit >>= 1) {
                j -= bit;
            }
            j += bit;
            if (i < j) {
                Complex temp = a[i];
                a[i] = a[j];
                a[j] = temp;
            }
        }

        // 蝶形运算，自底向上合并
        for (int len = 2; len <= n; len <<= 1) {
            double ang = 2 * Math.PI / len * (invert ? -1 : 1);
            Complex wlen = new Complex(Math.cos(ang), Math.sin(ang));

            for (int i = 0; i < n; i += len) {
                Complex w = new Complex(1, 0);
                for (int j = 0; j < len / 2; j++) {
                    Complex u = a[i + j];
                    Complex v = a[i + j + len / 2].multiply(w);

                    a[i + j] = u.add(v);
                    a[i + j + len / 2] = u.subtract(v);

                    w = w.multiply(wlen);
                }
            }
        }

        // 逆变换需要除以n
        if (invert) {
            for (int i = 0; i < n; i++) {
                a[i] = a[i].divide(n);
            }
        }
    }

    /**
     * 多项式乘法：计算两个多项式的乘积
     * @param a 第一个多项式的系数（a[0]是常数项）
     * @param b 第二个多项式的系数
     * @return 乘积多项式的系数
     */
    public static long[] multiply(long[] a, long[] b) {
        // 确定结果多项式的次数并找到最近的2的幂
        int resultSize = a.length + b.length - 1;
        int n = 1;
        while (n < resultSize) {
            n <<= 1;
        }

        // 填充到长度n
        Complex[] fa = new Complex[n];
        Complex[] fb = new Complex[n];
        for (int i = 0; i < n; i++) {
            fa[i] = new Complex(i < a.length ? a[i] : 0, 0);
            fb[i] = new Complex(i < b.length ? b[i] : 0, 0);
        }

        // 正变换：系数表示 → 点值表示
        fft(fa, false);
        fft(fb, false);

        // 点值逐点相乘
        for (int i = 0; i < n; i++) {
            fa[i] = fa[i].multiply(fb[i]);
        }

        // 逆变换：点值表示 → 系数表示
        fft(fa, true);

        // 四舍五入取整（处理浮点误差）
        long[] result = new long[resultSize];
        for (int i = 0; i < resultSize; i++) {
            result[i] = Math.round(fa[i].real);
        }

        return result;
    }

    /**
     * 大整数乘法：利用FFT加速大数相乘
     * 将大整数按位拆分为多项式系数
     */
    public static String bigIntegerMultiply(String num1, String num2) {
        // 处理符号
        boolean negative = false;
        if (num1.startsWith("-")) {
            negative = !negative;
            num1 = num1.substring(1);
        }
        if (num2.startsWith("-")) {
            negative = !negative;
            num2 = num2.substring(1);
        }

        // 去除前导零
        num1 = num1.replaceFirst("^0+", "");
        num2 = num2.replaceFirst("^0+", "");
        if (num1.isEmpty() || num2.isEmpty()) return "0";

        // 将数字字符串转为系数数组（每位一个系数）
        int n1 = num1.length(), n2 = num2.length();
        long[] a = new long[n1];
        long[] b = new long[n2];

        for (int i = 0; i < n1; i++) {
            a[i] = num1.charAt(n1 - 1 - i) - '0';
        }
        for (int i = 0; i < n2; i++) {
            b[i] = num2.charAt(n2 - 1 - i) - '0';
        }

        // FFT多项式乘法
        long[] result = multiply(a, b);

        // 处理进位
        long carry = 0;
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < result.length || carry > 0; i++) {
            long sum = carry;
            if (i < result.length) sum += result[i];
            sb.append(sum % 10);
            carry = sum / 10;
        }

        // 反转并添加符号
        String ans = sb.reverse().toString().replaceFirst("^0+", "");
        if (ans.isEmpty()) ans = "0";
        return negative && !ans.equals("0") ? "-" + ans : ans;
    }

    /**
     * 计算两个序列的循环卷积
     */
    public static long[] convolution(long[] a, long[] b) {
        return multiply(a, b);
    }

    /**
     * 主程序：测试FFT多项式乘法
     */
    public static void main(String[] args) {
        System.out.println("===== 快速傅里叶变换(FFT)多项式乘法 =====\n");

        // 测试1：简单多项式乘法
        System.out.println("--- 测试1：多项式乘法 ---");
        long[] poly1 = {1, 2, 3};      // 1 + 2x + 3x^2
        long[] poly2 = {4, 5};         // 4 + 5x
        long[] product = multiply(poly1, poly2);

        System.out.print("A(x) = ");
        printPoly(poly1);
        System.out.print("B(x) = ");
        printPoly(poly2);
        System.out.print("A(x) * B(x) = ");
        printPoly(product);
        // 期望: 4 + 13x + 22x^2 + 15x^3

        // 测试2：大整数乘法
        System.out.println("\n--- 测试2：大整数乘法 ---");
        String[] testPairs = {
            "123456789", "987654321",
            "999999999999999999", "999999999999999999",
            "3141592653589793238462643383279502884197", "2718281828459045235360287471352662497757"
        };

        for (int i = 0; i < testPairs.length; i += 2) {
            String a = testPairs[i];
            String b = testPairs[i + 1];
            long start = System.currentTimeMillis();
            String result = bigIntegerMultiply(a, b);
            long end = System.currentTimeMillis();

            System.out.println(a + " * " + b);
            System.out.println("= " + result);
            System.out.println("耗时: " + (end - start) + "ms\n");
        }

        // 测试3：序列卷积
        System.out.println("--- 测试3：序列卷积 ---");
        long[] seq1 = {1, 2, 3, 4, 5};
        long[] seq2 = {1, 0, -1};
        long[] conv = convolution(seq1, seq2);
        System.out.print("序列1: ");
        printArray(seq1);
        System.out.print("序列2: ");
        printArray(seq2);
        System.out.print("卷积结果: ");
        printArray(conv);

        // 测试4：性能对比
        System.out.println("\n--- 测试4：性能对比 ---");
        int[] sizes = {100, 1000, 5000, 10000};
        for (int size : sizes) {
            long[] arr1 = new long[size];
            long[] arr2 = new long[size];
            for (int i = 0; i < size; i++) {
                arr1[i] = i % 10;
                arr2[i] = (i * 2) % 10;
            }

            long start = System.currentTimeMillis();
            multiply(arr1, arr2);
            long fftTime = System.currentTimeMillis() - start;

            start = System.currentTimeMillis();
            naiveMultiply(arr1, arr2);
            long naiveTime = System.currentTimeMillis() - start;

            System.out.printf("规模: %5d, FFT: %4dms, 朴素: %6dms, 加速比: %.1fx%n",
                size, fftTime, naiveTime, (double) naiveTime / fftTime);
        }
    }

    /**
     * 朴素O(n^2)多项式乘法（用于性能对比）
     */
    private static long[] naiveMultiply(long[] a, long[] b) {
        long[] result = new long[a.length + b.length - 1];
        for (int i = 0; i < a.length; i++) {
            for (int j = 0; j < b.length; j++) {
                result[i + j] += a[i] * b[j];
            }
        }
        return result;
    }

    private static void printPoly(long[] coeff) {
        StringBuilder sb = new StringBuilder();
        for (int i = coeff.length - 1; i >= 0; i--) {
            if (coeff[i] == 0) continue;
            if (sb.length() > 0 && coeff[i] > 0) sb.append(" + ");
            else if (coeff[i] < 0) sb.append(" - ");
            if (Math.abs(coeff[i]) != 1 || i == 0) sb.append(Math.abs(coeff[i]));
            if (i > 0) sb.append("x");
            if (i > 1) sb.append("^").append(i);
        }
        System.out.println(sb.length() > 0 ? sb.toString() : "0");
    }

    private static void printArray(long[] arr) {
        System.out.print("[");
        for (int i = 0; i < arr.length; i++) {
            System.out.print(arr[i]);
            if (i < arr.length - 1) System.out.print(", ");
        }
        System.out.println("]");
    }
}

蝶形运算详解

FFT的核心是蝶形运算，它将DFT分解为两个半长DFT的组合：

对于k = 0, ..., n/2 - 1:
    A_k = E_k + ω^k · O_k
    A_{k+n/2} = E_k - ω^k · O_k

其中E_k是偶数项子多项式的DFT，O_k是奇数项子多项式的DFT，ω是n次单位根。

这种结构使得FFT可以用迭代方式自底向上实现，避免了递归开销。

复杂度分析

数值稳定性与优化

浮点误差控制

FFT涉及大量复数运算，累积浮点误差可能导致结果偏差。常用策略：
– 结果四舍五入到最近整数
– 对模运算场景，使用NTT（数论变换）替代FFT
– 分块FFT：将大系数拆分为高位和低位分别处理

NTT（数论变换）

在模素数域中寻找原根替代复数单位根，完全避免浮点误差。常用模数：998244353 = 119 × 2^23 + 1，原根为3。

应用场景

1. 大整数乘法：Karatsuba和FFT是密码学中大数运算的基础
2. 信号处理：音频/图像的频域分析和滤波
3. 字符串匹配：利用卷积计算模式匹配（KMP的替代方案）
4. 组合数学：快速计算生成函数乘积
5. 概率论：独立随机变量和的分布计算

总结

FFT通过分治策略将多项式乘法的复杂度从二次降至线性对数级，是算法设计中”变换-运算-逆变换”范式的典范。其Java实现展示了复数运算、位逆序置换和蝶形合并的完整流程，为理解更高级的数论变换（NTT）和工程应用奠定了基础。