/*
 * fft.c
 * File revision 0
 * Multiply two polynomials via Fast Fourier Transform.
 * (c) 2000 Jacob Lundberg, jacob@chaos2.org
 *
 * This program must be linked against the math library (gcc -lm).
 */


#include <sys/time.h>
#include <stdlib.h>
#include <unistd.h>
#include <stdio.h>
#include <math.h>


/* Are we in profiling mode (timer output only)? */
#define  PROFILE                0


int _mul(__complex__ double *coeff, __complex__ double *left, __complex__ double *right, int degree);
int _out(__complex__ double *coeff, int degree);
int _fft(__complex__ double *new, __complex__ double *old, int size, int inverse);
int _brc(__complex__ double *new, __complex__ double *old, int size, int bits);
int _rev(unsigned int val, int bits);
int _nxt(int val);


int main(int argc, char **argv) {
/*
 * main()
 * Slurp up the args, run the multiplication, print the result.
 */

   int counter, degree;
   __complex__ double *coeff, *left, *right;
   struct timeval before, after;

   /* There must be at least a degree and two polynomials. */
   if(argc < 4) {
      printf("The minimal input is `fft 0 x y'\n");
      return(1);
   }

   /* Arg 1 is specified to be the degree of the polynomials. */
   degree = atoi(argv[1]);

   /* Make sure the degree given is consistant with the input. */
   if(2 * degree + 2 != argc - 2) {
      printf("Arguments disagree (degree != size).\n");
      return(1);
   }

   /* Allocate the memory for the arrays (and zero it). */
   left  = (__complex__ double *)malloc(2 * _nxt(2 * degree) * sizeof(__complex__ double));
   right = (__complex__ double *)malloc(2 * _nxt(2 * degree) * sizeof(__complex__ double));
   coeff = (__complex__ double *)malloc(2 * _nxt(2 * degree) * sizeof(__complex__ double));

   /* Zero the arrays. */
   for(counter = 0; counter < 2 * _nxt(2 * degree); ++counter)
      left[counter] = 0 + 0i;
   for(counter = 0; counter < 2 * _nxt(2 * degree); ++counter)
      right[counter] = 0 + 0i;
   for(counter = 0; counter < 2 * _nxt(2 * degree); ++counter)
      coeff[counter] = 0 + 0i;

   /* Read the left polynomial into its input array. */
   for(counter = 0; counter <= degree; ++counter)
      left[counter] = atof(argv[1 + 1 + degree - counter]);

   /* Read the right polynomial into its input array. */
   for(counter = 0; counter <= degree; ++counter)
      right[counter] = atof(argv[1 + 2 + 2 * degree - counter]);

   /* Do the multiplication. */
   gettimeofday(&before, NULL);
   _mul(coeff, left, right, degree);
   gettimeofday(&after, NULL);

   /* Print the results. */
   if(PROFILE)
      printf("Two polynomials of degree %i multiplied in %lu microseconds.\n", degree,
            1000000 * (after.tv_sec - before.tv_sec) + after.tv_usec - before.tv_usec);
   else
      _out(coeff, 2 * degree);

   free(left);
   free(right);
   free(coeff);

   return(0);

}


int _mul(__complex__ double *coeff, __complex__ double *left, __complex__ double *right, int degree) {
/*
 * _mul()
 * Multiply two polynomials.
 */

   int counter;
   __complex__ double *new_left, *new_right;

   /* Allocate memory for the clones. */
   new_left  = (__complex__ double *)malloc(2 * _nxt(1 + degree) * sizeof(__complex__ double));
   new_right = (__complex__ double *)malloc(2 * _nxt(1 + degree) * sizeof(__complex__ double));

   /* Do an FFT on each input. */
   _fft(new_left, left, 2 * (1 + degree), 0);
   _fft(new_right, right, 2 * (1 + degree), 0);

   /* Multiply the inputs into the left input. */
   for(counter = 0; counter < 2 * (1 + degree); ++counter)
      new_left[counter] *= new_right[counter];

   /* Reverse the FFT into the new array. */
   _fft(coeff, new_left, 2 * (1 + degree), 1);

   /* Free the extra mamory. */
   free(new_left);
   free(new_right);

   return(0);

}


int _out(__complex__ double *coeff, int degree) {
/*
 * _out()
 * Display a polynomial.
 */

   int degcnt = degree;

   for(; degcnt >= 0; --degcnt)
      if(__real__ coeff[degcnt])
         printf("%gx^%i ", __real__ coeff[degcnt], degcnt);

   printf("\n");

   return(0);

}


int _fft(__complex__ double *new, __complex__ double *old, int size, int inverse) {
/*
 * _fft()
 * Fast Fourier Transform.
 */

   double angle = 2.0 * M_PI, da, w;
   __complex__ double list[3], first, second, temp;
   int counter = size, bits = 0;
   int i, j, k, n, bs, be;

   /* Count bits. */
   for(counter = size; counter > 0; ++bits)
      counter >>= 1;

   if(1 << (bits - 1) == size)
      --bits;

   /* Reverse the rotation of inverses. */
   if(inverse) angle *= -1;

   /* Organize the new array. */
   _brc(new, old, size, bits);

   /* Sorry it's so dense.  :( */
   be = 1;
   for(bs = 2; bs <= size; bs <<= 1) {
      da = angle / (double)bs;
      __real__ first = cos(-da);
      __imag__ first = sin(-da);
      __real__ second = cos(-2 * da);
      __imag__ second = sin(-2 * da);
      w = 2 * __real__ first;
      for(i = 0; i < size; i += bs) {
         list[1] = first;
         list[2] = second;
         for (j = i, n = 0; n < be; ++j, ++n) {
            list[0] = w * list[1] - list[2];
            list[2] = list[1];
            list[1] = list[0];
            k = j + be;
            temp = list[0] * new[k];
            new[k] = new[j] - temp;
            new[j] += temp;
         }
      }
      be = bs;
   }

   /* Normalize the inverse case (that pesky 1/n). */
   if(inverse)
      for(counter = 0; counter < size; ++counter)
         new[counter] /= (double)size;

   return(0);

}


int _brc(__complex__ double *new, __complex__ double *old, int size, int bits) {
/*
 * _brc()
 * Perform a ``bit-reverse-copy''.
 */

   unsigned int temp;

   /* Reverse those bits! */
   for(temp = 0; temp <= size - 1; ++temp)
      new[_rev(temp, bits)] = old[temp];

   return(0);

}


int _rev(unsigned int val, int bits) {
/*
 * _rev()
 * Reverse the significant (lower) bits in val.
 */

   int counter = 0;

   if(val <= 0 || bits <= 1)
      return(val);

   for(counter = 0; counter < (bits >> 1); ++counter)
      val = /* Purge the old bits. */
            (val & ~(1 << counter | 1 << (bits - counter - 1))) |
            /* Shift the right bit left. */
            (val & (1 << counter)) << (bits - (counter << 1) - 1) |
            /* Shift the left bit right. */
            (val & (1 << (bits - counter - 1))) >> (bits - (counter << 1) - 1);

   return(val);

}


int _nxt(int val) {
/*
 * _nxt()
 * Find the next smallext value in 2^n : n memb of I.
 */

   int counter = 0;

   if(val <= 0)
      return(0);

   --val;
   for(; val; ++counter)
      val >>= 1;

   return(1 << counter);

}