fft.c

使用FFT实现的两个多项式相乘的算法。 输入文件： 第一行为(n-1) 第二行为第一个多项式系数序列 第三行为第二个多项式系数序列 系数序列的格式为：an,an-1,an-2 ,…, a1

#include "FFT.h"
#include <malloc.h>
#include <math.h>


void InitializePoly(poly* p, int m, int* co)
{
	int i;
	p->MaxExp = m;
	p->Coefficient = (int*)malloc(sizeof(int) * (m + 1));
	for(i = 0; i < m + 1; i++)
		p->Coefficient[i] = co[i];
}

MyComplex* FFT(int size, MyComplex* Array, MyComplex C)
{
	MyComplex* V = (MyComplex*)malloc(sizeof(MyComplex) * size);
	MyComplex *Array1, *Array2;
	MyComplex *U, *W;
	int i;

	if(size == 1)
	{
		V[0].real = Array[0].real;
		V[0].imagine = Array[0].imagine;
		return V;
	}
	else
	{
		Array1 = (MyComplex*)malloc(sizeof(MyComplex) * size / 2);
		Array2 = (MyComplex*)malloc(sizeof(MyComplex) * size / 2);
		for(i = 0; i < size; i+=2)
		{
			Array1[i / 2].imagine = Array[i].imagine;
			Array1[i / 2].real = Array[i].real;
			Array2[i / 2].imagine = Array[i + 1].imagine;
			Array2[i / 2].real = Array[i + 1].real;
		}
		U = FFT(size / 2, Array1, Exp(C, 2));
		W = FFT(size / 2, Array2, Exp(C, 2));
		for(i = 0; i < size / 2; i++)
		{
			V[i] = Add(U[i], Mult(Exp(C, i), W[i]));
			V[i + size / 2] = Sub(U[i], Mult(Exp(C, i), W[i]));
		}
		return V;
	}
}

MyComplex Exp(MyComplex c, int exp)
{
	MyComplex result;
	double angle = atan(c.imagine / c.real);
	double length = sqrt(c.imagine * c.imagine + c.real * c.real);
	if(c.real < 0)
		angle = angle + PI;
	angle = angle * 180 / PI;
	angle = angle * exp;
	angle = angle + 360;
	while(angle > 360)
		angle = angle - 360;
	//angle = (double)((int)angle % 360);
	angle = angle * PI / 180;
	result.real = length * cos(angle);
	result.imagine = length * sin(angle);
	return result;
}

MyComplex Mult(MyComplex a, MyComplex b)
{
	MyComplex result;
	result.real = a.real * b.real - a.imagine * b.imagine;
	result.imagine = a.real * b.imagine + a.imagine * b.real;
	return result;
}

MyComplex Add(MyComplex a, MyComplex b)
{
	MyComplex result;
	result.real = a.real + b.real;
	result.imagine = a.imagine + b.imagine;
	return result;
}

MyComplex Sub(MyComplex a, MyComplex b)
{
	MyComplex result;
	result.real = a.real - b.real;
	result.imagine = a.imagine - b.imagine;
	return result;
}

MyComplex Inv(MyComplex a)
{
	double length = sqrt(a.real * a.real + a.imagine * a.imagine);
	MyComplex result;
	result.real = a.real / length;
	result.imagine = -a.imagine / length;
	return result;
}

void Equal(MyComplex* a, MyComplex b)
{
	if(a != 0)
	{
		a->imagine = b.imagine;
		a->real = b.real;
	}
}

int GetUpExpo(int a)
{
	int b = 1;
	while(b < a)
		b *= 2;
	return b;
}