FFT参考文章1

FFT参考文章2

快速傅里叶变换实现将两个多项式的乘法从O(n2)O(n^2)时间复杂度降低到O(nlogn)O(nlogn)

主要思想为分治，实现为DFT(离散傅里叶变换)和IDFT(离散傅里叶逆变换)，将多项式的系数表达式和点值表达式建立联系

递归实现

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
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <stdlib.h>

#define max(a, b) a>b?a:b
#define MAX 1<<22
#define eps 1e-6
#define pi acos(-1.0)
typedef struct complex {
	double x;
	double y;
} Complex, *Complexprt;
Complex a[MAX], b[MAX];

void mulComplex(Complexprt w1, Complexprt w2) {//两个复数相乘给第一个数
	//(x+iy)*(x'+iy') = x*x'-y*y'+i(xy'+yx')
	double x = w1->x * w2->x - w1->y * w2->y;
	double y = w1->x * w2->y + w1->y * w2->x;
	w1->x = x;
	w1->y = y;
}

void FFT(Complex a[], int n, int inv) {//inv为虚部符号，inv=1为FFT，inv=-1为IFFT
	if (n == 1) {//只有一项就返回
		return;
	}
	int mid = n / 2;
	Complex A1[mid + 1];
	Complex A2[mid + 1];
	for (int i = 0; i < n; i += 2) {//多项式拆分
		A1[i / 2] = a[i];
		A2[i / 2] = a[i + 1];
	}
	FFT(A1, mid, inv);
	FFT(A2, mid, inv);//求解子多项式
	Complex w0 = (Complex) {
		1, 0
	};
	Complex wn = (Complex) {
		cos(2 * pi / n), inv * sin(2 * pi / n)
	};
	for (int i = 0; i < mid; i++, mulComplex(&w0, &wn)) {
		//a[i] = A1[i] + w0 * A2[i];
		a[i].x = A1[i].x + w0.x * A2[i].x - w0.y * A2[i].y;
		a[i].y = A1[i].y + w0.x * A2[i].y + w0.y * A2[i].x;
		//a[i + n/2] = A1[i] -w0 * A2[i]
		a[i + n / 2].x = A1[i].x - w0.x * A2[i].x + w0.y * A2[i].y;
		a[i + n / 2].y = A1[i].y - w0.x * A2[i].y - w0.y * A2[i].x;
	}
}

int main() {
	//freopen("in.txt", "r", stdin);
	//freopen("out.txt", "w", stdout);
	int n, m;
	scanf("%d%d", &n, &m);
	for (int i = 0; i < n; i++) {//输入第一个系数多项式
		double x;
		scanf("%lf", &x);
		a[i].x = x;
	}
	for (int i = 0; i < m; i++) {//输入第二个系数多项式
		double x;
		scanf("%lf", &x);
		b[i].x = x;
	}
    //ceil为向上取整
	int len = max(1, (int) ceil(log2(n + m)));//FFT需要项数为2的整数次幂，所以多项式次数len为第一个大于n+m的2的整数次幂
	len = 1 << len;
	FFT(a, len, 1);//系数表达式转点值表达式
	FFT(b, len, 1);//系数表达式转点值表达式
	for (int i = 0; i <= len; i++) {//O(n)乘法
		mulComplex(&a[i], &b[i]);
	}
	FFT(a, len, -1);//点值表达式转系数表达式，逆变换之后需要除以n才是系数
	for (int i = 0; i < n + m - 1; i++) { //最高n+m-2次幂
		printf("%.0f ", a[i].x / len + eps);//输出系数除以n
	}
	return 0;
}

迭代实现

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
#define eps 1e-7
#define pi acos(-1.0)
typedef long double ld;
typedef struct complex {
    double x;
    double y;
} Complex, *Complexprt;
Complex a[4000005], b[4000005];
char A[1000005], B[1000005];
int rev[4000005];
int ans[4000005];

void mulComplex(Complexprt w1, Complexprt w2) {//两个复数相乘给第一个数
    //(x+iy)*(x'+iy') = x*x'-y*y'+i(xy'+yx')
    double x = w1->x * w2->x - w1->y * w2->y;
    double y = w1->x * w2->y + w1->y * w2->x;
    w1->x = x;
    w1->y = y;
}

void swap(Complex *a, Complex *b) {
    Complex p = *a;
    a->x = b->x;
    a->y = b->y;
    b->x = p.x;
    b->y = p.y;
}

void FFT(Complex a[], int n, int inv) {
    for (int i = 0; i < n; i++) {
        if (i < rev[i])swap(&a[i], &a[rev[i]]);
    }
    for (int h = 1; h < n; h <<= 1) {

        Complex wn = (Complex) {
                cos(pi / h), inv * sin(pi / h)
        };
        for (int j = 0; j < n; j += h << 1) {
            Complex w0 = (Complex) {
                    1, 0
            };
            for (int k = j; k < j + h; k++) {
                Complex x = a[k];
                Complex y = (Complex) {w0.x * a[k + h].x - w0.y * a[k + h].y, w0.x * a[k + h].y + w0.y * a[k + h].x};
                a[k] = (Complex) {x.x + y.x, x.y + y.y};
                a[k + h] = (Complex) {x.x - y.x, x.y - y.y};
                mulComplex(&w0, &wn);
            }
        }
    }
    if (inv == -1) {
        for (int i = 0; i < n; i++) {
            a[i].x = a[i].x / n + eps;
        }
    }
}
int main(){
    int len = 1;
        int l = 0;
        while (len < n + m - 1) {//len是项数,n,m也是项数
            len <<= 1;
            l++;//l是幂次
        }
        for (int i = 0; i < len; i++) {
            rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
        }
    FFT(a,len,1);
    FFT(b,len,1);
    for (int i = 0; i < len; i++) {
            mulComplex(&a[i], &b[i]);
        }
    FFT(a, len, -1);
}