FFT加速多项式乘法C语言版（基2FFT）

本文代码中FFT使用递归版本实现

FFT加速多项式乘法原理不多说了，直接贴代码如下：

在vs2017上测试成功

#include "pch.h"
#define _CRT_SECURE_NO_WARNINGS
#include "stdlib.h"
#include "math.h"
#include "stdio.h"


#define N 8
#define MAXN 100

#define Pi  3.1415927   //定义圆周率Pi
#define LEN sizeof(struct Compx)  //定义复数结构体大小

//-----定义复数结构体-----------------------
typedef struct Compx
{
	double real;
	double imag;
}Compx;

//-----复数乘法运算函数---------------------
struct Compx mult(struct Compx b1, struct Compx b2)
{
	struct Compx b3;
	b3.real = b1.real*b2.real - b1.imag*b2.imag;
	b3.imag = b1.real*b2.imag + b1.imag*b2.real;
	return(b3);
}

//-----复数减法运算函数---------------------
struct Compx sub(struct Compx a, struct Compx b)
{
	struct Compx c;
	c.real = a.real - b.real;
	c.imag = a.imag - b.imag;
	return(c);
}

//-----复数加法运算函数---------------------
struct Compx add(struct Compx a, struct Compx b)
{
	struct Compx c;
	c.real = a.real + b.real;
	c.imag = a.imag + b.imag;
	return(c);
}

void fft(Compx *a, int n, int inv);

int main()
{
	int i;

	int x[N] = { 0 }, y[N] = { 0 };

	printf("
N=%d
" , N);
	printf("
输入两个多项式的系数，输入系数为N多项式长度的一半
");
	printf("
输入第一个多项式的系数
");
	for (i = 0; i < N/2; i++)
	{
		scanf("%d", &x[i]);
	}

	printf("
输入第二个多项式的系数
");
	for (i = 0; i < N/2; i++)
	{
		scanf("%d", &y[i]);
	}

	struct  Compx * a = (struct Compx *)malloc(N*LEN);			//为结构体分配存储空间
	struct  Compx * b = (struct Compx *)malloc(N*LEN);
	struct  Compx * c = (struct Compx *)malloc(N*LEN);

	//初始化======================================= 
	printf("
a多项式初始化：
");
	for (i = 0; i < N; i++)
	{
		a[i].real = x[i];
		a[i].imag = 0;
		printf("%.4f ", a[i].real);
		printf("+%.4fj  ", a[i].imag);
		printf("
");
	}
	printf("
b多项式初始化：
");
	for (i = 0; i < N; i++)
	{
		b[i].real = y[i];
		b[i].imag = 0;
		printf("%.4f ", b[i].real);
		printf("+%.4fj  ", b[i].imag);
		printf("
");
	}
	printf("
结果c多项式初始化：
");
	for (i = 0; i < N; i++)
	{
		c[i].real = 0;
		c[i].imag = 0;
		printf("%.4f ", c[i].real);
		printf("+%.4fj  ", c[i].imag);
		printf("
");
	}

	fft(a, N, 1);

	printf("
第一个多项式FFT计算结果：
");
	for (i = 0; i < N; i++)
	{
		printf("%.4f ", a[i].real);
		printf("+%.4fj  ", a[i].imag);
		printf("
");
	}

	fft(b, N, 1);

	printf("
第二个多项式FFT计算结果：
");
	for (i = 0; i < N; i++)
	{
		printf("%.4f ", b[i].real);
		printf("+%.4fj  ", b[i].imag);
		printf("
");
	}

	for (i = 0; i < N; i++)
		a[i] = mult(a[i] , b[i]);

	fft(a, N, -1);
	for (i = 0; i < N; i++) {
		c[i].real = a[i].real / N;
		c[i].imag = a[i].imag / N;
	}
		
	printf("
乘积多项式结果：
");
	for (i = 0; i < N; i++)
	{
		printf("%.4f ", c[i].real);
		printf("+%.4fj  ", c[i].imag);
		printf("
");
	}

	return 0;
}

void fft(Compx *a, int n, int inv) {
	if (n == 1)return;
	int mid = n / 2;
	static Compx b[MAXN];
	int i;
	for (i = 0; i < mid; i++) {
		b[i] = a[i * 2];
		b[i + mid] = a[i * 2 + 1];
	}
		
	for (i = 0; i < n; i++)
		a[i] = b[i];
	fft(a, mid, inv);
	fft(a + mid, mid, inv);//分治
	for (i = 0; i < mid; i++)
	{
		Compx x;
		x.real = cos(2 * Pi*i / n);
		x.imag = inv * sin(2 * Pi*i / n);

		b[i] = add(a[i], mult(x, a[i + mid]));
		b[i + mid] = sub(a[i] , mult(x , a[i + mid]));
	}
	for (i = 0; i < n; i++)
		a[i] = b[i];
}