/*
     widentt.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 NTT multiplication using very
** wide primes.  The code is written for simplicity, not
** for performance, so it is *EXTREMELY* slow.
**
** To compile this using GCC:
** gcc main.c widentt.c -o widentt.exe
*/

/*
** How many bits should be in our ModInt?
** Choice of 64, 253, 254 and 256
*/
#define MOD_BITS 64

#include 
#include 
#include 
#include 

#include "wmodmath.h"
size_t ModWords; /* How many UINT16 words are our ModInts? */

typedef short int Short;

#define CalcNTTLen(_NumLen) ((((_NumLen)*BASE_DIG)*2)/BASE_DIG)
/* NumLen*BaseDig*ZeroPadding/Dig_Per_NTT */

static ModIntPtr NTTNum1=NULL, NTTNum2=NULL;
static int BASE;
static int BASE_DIG;

ModIntPtr Prime,PrimvRoot,MulInv;




static void
ComputeRoot(ModIntPtr w,size_t Step,int Dir,size_t ModWords)
/*
** Compute the trig root for this step size.
*/
{ModIntPtr Temp;

 Temp=ModMalloc(1,ModWords);

if (Dir > 0)
  {
   CB_Copy(Temp,ModWords,Prime,ModWords);CB_Sub1(Temp,ModWords);
   CB_DivInt(Temp,Step,ModWords);CB_Sub(Temp,Prime,Temp,ModWords);
   CB_Sub1(Temp,ModWords);
   ModPow(w, PrimvRoot,Temp,ModWords);
  }
else
  {
   CB_Copy(Temp,ModWords,Prime,ModWords);CB_Sub1(Temp,ModWords);
   CB_DivInt(Temp,Step,ModWords);
   ModPow(w, PrimvRoot,Temp,ModWords);
  }

/*
 ModPow(NthRoot, PrimvRoot,Prime-1-(Prime-1)/NTTLen);
 ModPow(Nthroot1,PrimvRoot,(Prime-1)/NTTLen);
*/
 ModFree(Temp);
}

static void
NTTReOrder(ModIntPtr Data, size_t Len,size_t ModWords)
{size_t Index,xednI,k;

xednI=0;
for (Index=0;Index Index)
     ModSwap(&Data[Index*ModWords],&Data[xednI*ModWords],ModWords);

   k=Len/2;
   while ((k <= xednI) && (k >=1)) {xednI-=k;k/=2;}
   xednI+=k;
  }
}

static void
NTT(ModIntPtr Data, size_t Len,int Dir,size_t ModWords)
/*
** A Radix 2, decimation in time, iterative NTT.
*/
{size_t j,Step,HalfStep,z;
 ModIntPtr u,w,T,End;

NTTReOrder(Data,Len,ModWords);

Step=1;z=1;
u=ModMalloc(1,ModWords);
w=ModMalloc(1,ModWords);
T=ModMalloc(1,ModWords);
End=&Data[Len*ModWords];

while (Step < Len)
  {
   HalfStep=Step;
   Step*=2;

   ModSet1(u,ModWords);
   ComputeRoot(w,Step,Dir,ModWords);
   for (j=0;j 4)
  {
   printf("Error:  The fft is slightly hardwired for <= 4 digits in the base.\n");
   exit(0);
  }

if ((NTTNum1==NULL) || (NTTNum2==NULL))
   {
    printf("Unable to allocate memory for WideNTTNum1, WideNTTNum2.\n");
    printf("Len=%d Words=%d\n",(int)Len,(int)ModWords );
    exit(0);
   }

Prime=ModMalloc(1,ModWords);
PrimvRoot=ModMalloc(1,ModWords);
MulInv=ModMalloc(1,ModWords);

#if (MOD_BITS==256)
/* 256 bits */
HexStr2Mod(Prime,"CF00000000000000000000000000000000000000000000000000000000000001",ModWords); /* 207*2^248+1 */
HexStr2Mod(PrimvRoot,"5",ModWords); /* primitive root=5 */
#elif (MOD_BITS==254)
/* 254 bits */
HexStr2Mod(Prime,"4083000000000000000000000000000000000000000000000000000000000001",ModWords); /* 16515*2^240+1 */
HexStr2Mod(PrimvRoot,"d",ModWords); /* primitive root=13 */
#elif (MOD_BITS==253)
/* 253 bits */
HexStr2Mod(Prime,"1FE00000000000000000000000000000000000000000000000000000000001",ModWords); /* 255*2^237+1 */
HexStr2Mod(PrimvRoot,"13",ModWords); /* primitive root=19 (0x13) */
#elif (MOD_BITS==64)
/* 64 bits */
HexStr2Mod(Prime,"ffffffff00000001",ModWords); /* 18446742974197923841 */
HexStr2Mod(PrimvRoot,"7",ModWords);
#else
#error Unknown number of MOD_BITS
#endif
}

void DeInitFFT(unsigned long int Len)
{
ModFree(NTTNum1);
ModFree(NTTNum2);
ModFree(Prime);
ModFree(PrimvRoot);
ModFree(MulInv);
DeInitModMath();
}