?
Solved

Polynomial Multiplication with FFT

Posted on 2006-05-28
4
Medium Priority
?
804 Views
Last Modified: 2012-05-05
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!
Please Help!!!
0
Comment
Question by:Gonella
2 Comments
 
LVL 53

Accepted Solution

by:
Infinity08 earned 1000 total points
ID: 16779807
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
 
LVL 2

Assisted Solution

by:Tussin
Tussin earned 1000 total points
ID: 16788345
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

Featured Post

Concerto Cloud for Software Providers & ISVs

Can Concerto Cloud Services help you focus on evolving your application offerings, while delivering the best cloud experience to your customers? From DevOps to revenue models and customer support, the answer is yes!

Learn how Concerto can help you.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

This tutorial is posted by Aaron Wojnowski, administrator at SDKExpert.net.  To view more iPhone tutorials, visit www.sdkexpert.net. This is a very simple tutorial on finding the user's current location easily. In this tutorial, you will learn ho…
Windows programmers of the C/C++ variety, how many of you realise that since Window 9x Microsoft has been lying to you about what constitutes Unicode (http://en.wikipedia.org/wiki/Unicode)? They will have you believe that Unicode requires you to use…
Video by: Grant
The goal of this video is to provide viewers with basic examples to understand and use while-loops in the C programming language.
The goal of this video is to provide viewers with basic examples to understand and use switch statements in the C programming language.
Suggested Courses

840 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question