• Status: Solved
• Priority: Medium
• Security: Public
• Views: 827

# Polynomial Multiplication with FFT

I need working code in  C to multiply two polinomials using Fast Fourrier Transform (FFT) in recursive way.

Ex.
I have two arrays of n elements:
int *poli1 = { 2, 3, 1 } /* 2x^2 + 3x + 1 */
int *poli2 = { 3, 4, 8 } /* 3x^2 + 4x + 8 */

I need an array R of results:  R = { 6, 17, 31, 28, 8 }  /* 6x^4 + 17x^3 + 31x^2 + 28x + 8 */

with FFT the execution time is something like O(n^1,58)

This is URGENT!
0
Gonella
2 Solutions

Commented:
As this looks like homework, I'm not gonna hand you any code on a platter, but I'll give you an explanation of the mathematical principal behind multiplication (of polynomials eg.) using FFT :

http://en.wikipedia.org/wiki/Multiplication_algorithm

scroll down to "Fourier transform methods". It gives 4 steps of how to do it.

Do you have any specific problems implementing this ? What code do you have already ?
0

Commented:
I wrote this code based on Xavier Gourdon code.
The key function is       "MulBigInt"

The main function is this program shows multiplication of two polinomials using Fast Fourrier Transform (FFT) (recursive) with n=250000

Tussin

/*
* Xavier Gourdon : Sept. 99 (xavier.gourdon.free.fr)
*
* Bigc : Basic file for manipulation of Large Integers.
*            It rely on FFT.c to multiply number with the FFT
*            technique.
*
* Several optimizations can be made :
*  - Implementation of a classic multiplication for small BigInts.
*  - Avoid using UpdateBigInt for the sum of BigInts.
*  - Use a larger base for the internal data storage of BigInts to
*    save space.
*  ...
*
*  Informations can be found on
*    http://xavier.gourdon.free.fr/Constants/constants.html
*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

typedef struct BigIntTag {
double * Coef;
long Size, SizeMax;
} BigInt;

/*
* Decrease the base if the worst error in FFT become too large
*/
#define BASE 10000.
#define invBASE 0.0001
#define NBDEC_BASE 4

extern double FFTSquareWorstError;

#define PI 3.1415926535897932384626

typedef struct ComplexTag {
double R,I;
} Complex;

int FFTLengthMax;
Complex * OmegaFFT;
Complex * ArrayFFT0, * ArrayFFT1;
Complex * ComplexCoef;
Complex * MainCoef,* MainFFT;
double FFTSquareWorstError;

void PrintComplexList(Complex *A,int Size) {
for (int i=0;i<Size;i++) {
printf("i=%d ; R=%f, I=%f\n",i,A[i].R,A[i].I);
}
}

void InitializeFFT(int MaxLength)
{
int i;
double Step;

FFTLengthMax = MaxLength;
OmegaFFT = (Complex *) malloc(MaxLength/2*sizeof(Complex));
ArrayFFT0 = (Complex *) malloc(MaxLength*sizeof(Complex));
ArrayFFT1 = (Complex *) malloc(MaxLength*sizeof(Complex));
ComplexCoef = (Complex *) malloc(MaxLength*sizeof(Complex));
Step = 2.*PI/(double) MaxLength;
for (i=0; 2*i<MaxLength; i++) {
OmegaFFT[i].R = cos(Step*(double)i);
OmegaFFT[i].I = sin(Step*(double)i);
}
//PrintComplexList(OmegaFFT,MaxLength >> 1);
FFTSquareWorstError=0.;
}

int aaa=0;

void RecursiveFFT(Complex * Coef, Complex * FFT, int Length, int Step,
int Sign, char * comment)
{
int i, OmegaStep;
Complex * FFT0, * FFT1, * Omega;
double tmpR, tmpI;

//System.out.println("Recursive "+CoefStartPos+" "+FFTStartPos+" "+Length+ " "+Step+" "+Sign+" "+comment);
//printf("Recursive %d %d %d %d %d %s\n",Coef-MainCoef,FFT-MainFFT,Length,Step,Sign,comment);
if (Length==2) {
FFT[0].R = Coef[0].R + Coef[Step].R;
FFT[0].I = Coef[0].I + Coef[Step].I;
FFT[1].R = Coef[0].R - Coef[Step].R;
FFT[1].I = Coef[0].I - Coef[Step].I;
//printf("zzz=%f : %f : %f : %f\n",FFT[0].R,FFT[0].I,FFT[1].R,FFT[1].I);
//printf("xxx=%f\n",MainFFT[3].R);
aaa++;
return;
}

FFT0 = FFT;
RecursiveFFT(Coef     ,FFT0,Length/2,Step*2,Sign,"upper");
FFT1 = FFT+Length/2;
//printf("step=%d\n",Step);
RecursiveFFT(Coef+Step,FFT1,Length/2,Step*2,Sign,"lower");

//printf("hi\n");
//PrintComplexList(MainFFT,16);

Omega = OmegaFFT;
OmegaStep = FFTLengthMax/Length;
for (i=0; 2*i<Length; i++, Omega += OmegaStep) {
/* Recursion formula for FFT :
FFT[i]          <-  FFT0[i] + Omega*FFT1[i]
FFT[i+Length/2] <-  FFT0[i] - Omega*FFT1[i],
Omega = exp(2*I*PI*i/Length) */
tmpR = Omega[0].R*FFT1[i].R-Sign*Omega[0].I*FFT1[i].I;
tmpI = Omega[0].R*FFT1[i].I+Sign*Omega[0].I*FFT1[i].R;
FFT1[i].R = FFT0[i].R - tmpR;
FFT1[i].I = FFT0[i].I - tmpI;
FFT0[i].R = FFT0[i].R + tmpR;
FFT0[i].I = FFT0[i].I + tmpI;
//if ((FFT0-MainFFT+i >=22 && FFT0-MainFFT+i <=28) ||(FFT1-MainFFT+i >=22 && FFT1-MainFFT+i <=28))
//printf("effect %d %d\n",FFT1-MainFFT+i,FFT0-MainFFT+i);

//      printf("effect %d %d\n",FFT1-MainFFT+i,FFT0-MainFFT+i);
}
//  printf("hi2\n");

//PrintComplexList(MainFFT,64);
// printf("%f %f\n",MainFFT[22].R,MainFFT[22].I);
}

/* Compute the complex Fourier Transform of Coef into FFT */
void FFT(double * Coef, int Length, Complex * FFT, int NFFT)
{
int i;
/* Transform array of real coefficient into array of complex */
for (i=0; i<Length; i++) {
ComplexCoef[i].R = Coef[i];
ComplexCoef[i].I = 0.;
}
for (; i<NFFT; i++)
ComplexCoef[i].R = ComplexCoef[i].I = 0.;

MainCoef=ComplexCoef;
MainFFT=FFT;

RecursiveFFT(ComplexCoef,FFT,NFFT,1,1,"start");
}

/* Compute the inverse Fourier Transform of FFT into Coef */
void InverseFFT(Complex * FFT, int NFFT, double * Coef, int Length)
{
int i;
double invNFFT = 1./(double) NFFT, tmp;

MainCoef=FFT;
MainFFT=ComplexCoef;
RecursiveFFT(FFT, ComplexCoef, NFFT, 1, -1,"inverse");
for (i=0; i<Length; i++) {
/* Closest integer to ComplexCoef[i].R/NFFT */
tmp = invNFFT*ComplexCoef[i].R;
Coef[i] = floor(0.5+tmp);
if ((tmp-Coef[i])*(tmp-Coef[i])>FFTSquareWorstError)
FFTSquareWorstError = (tmp-Coef[i])*(tmp-Coef[i]);
}
}

void Convolution(Complex * A, Complex * B, int NFFT, Complex * C)
{
int i;
double tmpR, tmpI;

for (i=0; i<NFFT; i++) {
tmpR = A[i].R*B[i].R-A[i].I*B[i].I;
tmpI = A[i].R*B[i].I+A[i].I*B[i].R;
C[i].R = tmpR;
C[i].I = tmpI;
}
}

void MulWithFFT(double * ACoef, int ASize,
double * BCoef, int BSize,
double * CCoef)
{
int NFFT = 2;

while (NFFT<ASize+BSize)
NFFT *= 2;

if (NFFT>FFTLengthMax) {
printf("Error, FFT Size is too big in MulWithFFT\n");
return;
}
FFT(ACoef, ASize, ArrayFFT0, NFFT);
//printf("xxx\n");
//PrintComplexList(ArrayFFT0,NFFT);
FFT(BCoef, BSize, ArrayFFT1, NFFT);
//PrintComplexList(ArrayFFT1,NFFT);
Convolution(ArrayFFT0,ArrayFFT1,NFFT,ArrayFFT0);
InverseFFT(ArrayFFT0,NFFT,CCoef, ASize+BSize-1);
}

void InitializeBigInt(BigInt * A, int MaxSize)
{
A->Coef = (double *) malloc(MaxSize*sizeof(double));

if (A->Coef==0)
printf("Not enough memory\n");
A->Size = 0;
A->SizeMax = MaxSize;
}

void FreeBigInt(BigInt * A)
{
free(A->Coef);
}

void PrintBigInt(BigInt * A)
{
int i,j,Digit=0,Dec;
double Pow, Coef;

return;

if (A->Size==0) {
printf("0");
return;
}

for (i=A->Size-1; i>=0; i--) {
Pow = BASE*0.1;
Coef = A->Coef[i];
for (j=0; j<NBDEC_BASE; j++) {
//      printf("c%f p%f",Coef,Pow);
Dec = (int) (Coef/Pow);
Coef -= Dec*Pow;
Pow *= 0.1;
printf("%ld",Dec);
Digit++;
//  if (Digit%10==0) printf(" ");
//  if (Digit%50==0) printf(": %ld\n",Digit);
}
}
}

/*
* Update A to the normal form (0<= A.Coef[i] < BASE)
*/
void UpdateBigInt(BigInt * A)
{
int i;
double carry=0., x;

for (i=0; i<A->Size; i++) {
x = A->Coef[i] + carry;
carry = floor(x*invBASE);
A->Coef[i] = x-carry*BASE;
}
if (carry!=0) {
while (carry!=0.) {
x = carry;
carry = floor(x*invBASE);
A->Coef[i] = x-carry*BASE;
i++;
A->Size=i;
}
if (A->Size > A->SizeMax) {
printf("Error in UpdateBigInt, Size>SizeMax\n");
}
}
else {
while (i>0 && A->Coef[i-1]==0.) i--;
A->Size=i;
}
}

/*
* Compute C = A*B
*/
void MulBigInt(BigInt * A, BigInt * B, BigInt * C)
{
MulWithFFT(A->Coef,A->Size,B->Coef,B->Size,C->Coef);
C->Size = A->Size+B->Size-1;
UpdateBigInt(C);
}

void main() {

BigInt a,b,c;

int i =65536;

InitializeFFT(524288);
InitializeBigInt(&a,262144);
for (int i=0;i<250000;i++)
a.Coef[i]=9999;
a.Size=250000;
InitializeBigInt(&b,262144);
for (int i=0;i<250000;i++)
b.Coef[i]=9999;
b.Size=250000;
InitializeBigInt(&c,524288);
//      printf("\na=");
//      PrintBigInt(&a);
//      printf("\nb=");
//      PrintBigInt(&b);
MulBigInt(&a,&b,&c);
//      printf("\nc=");
PrintBigInt(&c);
//      printf("%d\n",aaa);
}

0
Question has a verified solution.

Are you are experiencing a similar issue? Get a personalized answer when you ask a related question.

Have a better answer? Share it in a comment.