/*
     ss2.h
*/
/*
** 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 is a support file for SS2.c
*/

typedef signed long    INT32; /* 32/31 bit signed int */
typedef INT32 ModInt31;


ModInt31 Prime31=2130706433;
ModInt31 PrimvRoot31=3;
ModInt31 MulInv31;
static double RecipPrime31=1.0/2130706433.0;


ModInt31
ModMul31(ModInt31 a, ModInt31 b)
/* Limited to 31 bits. */
{ModInt31 rem;
rem = a * b;
rem = rem - Prime31 * ((ModInt31) floor(0.5+RecipPrime31 * ((double) a) * ((double) b)));
if (rem < 0) rem +=Prime31;
return rem;
}

ModInt31 ModAdd31(ModInt31 a,ModInt31 b)
{UINT32 Sum; /* Big enough to hold sum */
Sum=a+b;
if (Sum >= Prime31) Sum-=Prime31;
return (ModInt31)Sum;
}

ModInt31 ModSub31(ModInt31 a,ModInt31 b)
{INT32 Dif; /* Big enough to hold signed difference */
Dif=a-b;
if (Dif < 0) Dif+=Prime31;
return (ModInt31)Dif;
}


ModInt31
ModPow31(ModInt31 Base,ModInt31 Expon)
{ModInt31 prod,b;

if (Expon<=0) return 1;

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

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

ModInt31
FindInverse31(ModInt31 Num)
{ModInt31 i;
i=ModPow31(Num,Prime31-2);
return i;
}

static void
NTTReorder31(ModInt31 *Data, int Len)
{int Index,xednI,k;

xednI = 0;
for (Index = 0;Index < Len;Index++)
  {
   if (xednI > Index)
     {ModInt31 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 NTT31(ModInt31 *Data, int Len, int Dir)
/* A simple minded, generic transform */
{int j,step,halfstep;
 int index,index2;
 ModInt31 u,w,temp;

NTTReorder31(Data,Len);

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

   u=1;
   if (Dir > 0) w=ModPow31(PrimvRoot31,Prime31-1-((Prime31-1)/step));
   else         w=ModPow31(PrimvRoot31,(Prime31-1)/step);

   for (j=0;j ((UINT32)Prime31) )
    {
     printf("SmallNTTMul called with too big of a base: %d\n",BASE);
     exit(0);
    }

  MulInv31  =FindInverse31(Len);
  Bytes=sizeof(ModInt31)*Len;

  NTTNum1=(ModInt31*)malloc(Bytes);
  NTTNum2=(ModInt31*)malloc(Bytes);

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

  for (x = 0; x < Len; x++)
    NTTNum1[x] = ModGet32(Num1+x*ModWords,ModWords) % BASE;
  NTT31(NTTNum1, Len, 1);

  if (Num1 == Num2)
    {
      /*
      ** Now do the convolution
      */
      for (x=0;x< Len;x++)
        NTTNum1[x]=ModMul31(NTTNum1[x],NTTNum1[x]);
    }
  else
    {
      for (x = 0; x < Len; x++)
        NTTNum2[x] = ModGet32(Num2+x*ModWords,ModWords) % BASE;
      NTT31(NTTNum2, Len, 1);

      /*
      ** Now do the convolution
      */
      for (x=0;x< Len;x++)
        NTTNum1[x]=ModMul31(NTTNum1[x],NTTNum2[x]);
    }

  NTT31(NTTNum1, Len, -1);

/* Don't release our carries.  Just make it modulo oure BASE */
  for (x = 0; x < Len; x++)
      ModSet(Prod+x*ModWords,
             ModMul31(NTTNum1[x],MulInv31) % BASE,
             ModWords);

free(NTTNum1);free(NTTNum2);
}