/*
     fgt.c
*/
/*
** This file is placed into the public domain by its author,
** Carey Bloodworth (Carey@Bloodworth.org) on July 16, 2001
**
** This multiplication demo is not designed for high performance.
** It's a tutorial program designed to be used with the information
** on my web site at www.Bloodworth.org
*/
/*
** This file demonstrates a very basic FGT multiply.  It uses a single
** 31 bit prime.  That means you can only put a single decimal in there
** without things overflowing quickly.  It does, however, mean that
** everybody can experiment with this.
**
** To compile this using GCC:
** gcc main.c fgt.c -o fgt.exe
*/
#include 
#include 
#include 
#include 
#include 

#define CalcFGTLen(_NumLen) ((((_NumLen)*BASE_DIG)*2)/(BASE_DIG*2))
/* NumLen*BaseDig*ZeroPadding/Dig_Per_FFT */
/* It's a complex transform, so two data pieces per FGT element */

typedef short int Short;
typedef signed long    INT32; /* 32/31 bit signed int */
typedef unsigned long UINT32; /* 32 bit unsigned int */
typedef INT32 ModInt;
typedef struct {ModInt r,i;} CmplxModInt;

static CmplxModInt *FGTNum1=NULL, *FGTNum2=NULL;
static int BASE;
static int BASE_DIG;

ModInt Prime,MulInv;
CmplxModInt NthRoot,NthRoot1,PrimvRoot;
static double RecipPrime=0.0;
int PrimeQ;

ModInt
ModMul(ModInt a, ModInt b)
/* Limited to 31 bits. */
{ModInt rem;
rem = a * b;
rem = rem - Prime * ((ModInt) floor(0.5+RecipPrime * ((double) a) * ((double) b)));
if (rem < 0) rem +=Prime;
return rem;
}

ModInt ModAdd(ModInt a,ModInt b)
{double Sum; /* Big enough to hold sum */
Sum=((double)a)+((double)b);
if (Sum >= Prime) Sum-=Prime;
return (ModInt)Sum;
}

ModInt ModSub(ModInt a,ModInt b)
{double Dif; /* Big enough to hold signed difference */
Dif=((double)a)-((double)b);
if (Dif < 0) Dif+=Prime;
return (ModInt)Dif;
}


ModInt
ModPow(ModInt Base,ModInt Expon)
{ModInt prod,b;

if (Expon<=0) return 1;

b=Base;
while (!(Expon&1)) {b=ModMul(b,b);Expon>>=1;}
prod=b;

while (Expon>>=1)
  {
   b=ModMul(b,b);
   if (Expon&1) prod=ModMul(prod,b);
  }
return prod;
}

CmplxModInt CmplxModConj(CmplxModInt M)
{
M.i=(Prime-M.i)%Prime;
return M;
}

CmplxModInt CmplxModMul(CmplxModInt Num1,CmplxModInt Num2)
{CmplxModInt P;
P.r=ModSub(ModMul(Num1.r,Num2.r),ModMul(Num1.i,Num2.i));
P.i=ModAdd(ModMul(Num1.r,Num2.i),ModMul(Num1.i,Num2.r));
return P;
}

CmplxModInt CmplxModAdd(CmplxModInt Num1,CmplxModInt Num2)
{CmplxModInt S;
S.r=ModAdd(Num1.r,Num2.r);
S.i=ModAdd(Num1.i,Num2.i);
return S;
}

CmplxModInt CmplxModSub(CmplxModInt Num1,CmplxModInt Num2)
{CmplxModInt D;
D.r=ModSub(Num1.r,Num2.r);
D.i=ModSub(Num1.i,Num2.i);
return D;
}

CmplxModInt
CmplxModPow(CmplxModInt base,UINT32 expon)
{CmplxModInt prod,b;

prod.r=1;prod.i=0;
if (expon<=0) return prod;

b=base;
while (!(expon&1))
  {
   b=CmplxModMul(b,b);
   expon>>=1;
  }
prod=b;

while (expon>>=1)
  {
   b=CmplxModMul(b,b);
   if (expon&1) prod=CmplxModMul(prod,b);
  }
return prod;
}

int Log2(int Num)
{int x=-1;
if (Num==0) return 0;
while (Num) {x++;Num/=2;}
return x;
}


static void
FGTReorder(CmplxModInt *Data, int Len)
{int Index,xednI,k;

xednI = 0;
for (Index = 0;Index < Len;Index++)
  {
   if (xednI > Index)
     {CmplxModInt Temp;
      Temp=Data[xednI];
      Data[xednI]=Data[Index];
      Data[Index]=Temp;
     }
   k=Len/2;
   while ((k <= xednI) && (k >=1)) {xednI-=k;k/=2;}
   xednI+=k;
  }
}

void FGT(CmplxModInt *Data, int Len, int Dir)
/* A simple minded, generic transform */
{int j,step,halfstep;
 int index,index2;
 CmplxModInt temp;
 CmplxModInt u,w;
 int power;

FGTReorder(Data,Len);

power=Len;
step=1;
while (step < Len)
  {
   halfstep=step;
   step*=2;power/=2;

   u.r=1;u.i=0;
   if (Dir > 0) w=NthRoot; else w=NthRoot1;
   w=CmplxModPow(w,power);
   for (j=0;j 0; x--)
    {
      Pyramid = FGTNum1[FGTLen - x].r + Carry;
      Carry = floor(Pyramid / BASE);
      Prod[x*2 - 1] = Pyramid - Carry * BASE;

      Pyramid = FGTNum1[FGTLen - x].i + Carry;
      Carry = floor(Pyramid / BASE);
      Prod[x*2 - 2] = Pyramid - Carry * BASE;
    }
}
}

void
InitFFT(unsigned long int Len,int Base,int BaseDig)
{int Bytes;

BASE=Base;
BASE_DIG=BaseDig;

Bytes=sizeof(CmplxModInt)*CalcFGTLen(Len);

if ( (BaseDig > 1) || (BASE > 10))
  {
   printf("Error:  The FGT is hardwired for just 1 digit in the base.\n");
   printf("In 'main.c' please set BASE to 10 and BASE_DIG to 1\n");
   exit(0);
  }

FGTNum1=(CmplxModInt*)malloc(Bytes);
FGTNum2=(CmplxModInt*)malloc(Bytes);

if ((FGTNum1==NULL) || (FGTNum2==NULL))
   {
    printf("Unable to allocate memory for NTTNum.\n");
    printf("Len=%d Bytes=%d\n",(int)Len,(int)Bytes);
    exit(0);
   }
}

void DeInitFFT(unsigned long int Len)
{
free(FGTNum1);free(FGTNum2);
}