#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/* Data types for matrices and vectors */
typedef char **matrix;
typedef double *vector;

/* Upon input of k, return H(k) */
matrix createHadMat ( unsigned int k )
{
   unsigned int i, j, s, t;
   matrix H;

   /* Allocate memory */
   H = (char **)malloc((1 << k) * sizeof(char *));
   for (i = 0; i < (1 << k); ++i)
      H[i] = (char *)malloc((1 << k) * sizeof(char));

   /* Make a copy of H(0) */
   H[0][0] = 1; t = 1;

   /* Iteratively populate H(1), H(2), ..., H(k) */
   for (s=1; s<=k; ++s) {
      for (i=0; i<t; ++i) {
         for (j=0; j<t; ++j) {
            H[t+i][j] = H[i][t+j] = H[i][j];
            H[t+i][t+j] = -H[i][j];
         }
      }
      t <<= 1;
   }

   return H;
}

/* Function for printing a Hadamard matrix */
void prnHadMat ( matrix H, unsigned int k )
{
   unsigned int i, j;

   printf("H(%u) =\n", k);
   for (i = 0; i < (1 << k); ++i) {
      for (j = 0; j < (1 << k); ++j) printf("%3d", H[i][j]);
      printf("\n");
   }
}

/* Function for generating a random vector of size n with double-precision
   entries in the range [0,1] */
vector genRandVec ( unsigned int n )
{
   unsigned int i;
   vector v;

   v = (double *)malloc(n * sizeof(double));
   for (i=0; i<n; ++i) v[i] = (double)rand() / (double)RAND_MAX;
   return v;
}

/* Function for printing a column vector */
void prnVec ( vector v, unsigned int n )
{
   unsigned int i;

   for (i = 0; i < n; ++i) printf("%20.10lf\n", v[i]);
}

/* Function for computing H(k)v using standard matrix-multiplication formula */
vector matVecProd ( matrix H, vector v, unsigned int n )
{
   unsigned int i, j;
   vector w;

   w = (double *)malloc(n * sizeof(double));
   for (i = 0; i < n; ++i) {
      /* Compute the dot product of the ith row of H and the vector v */
      w[i] = 0; for (j = 0; j < n; ++j) w[i] += H[i][j] * v[j];
   }
   return w;
}

/* Function for freeing a 2-dimensional matrix */
void freeMat ( matrix H, unsigned int n )
{
   unsigned int i;

   for (i=0; i<n; ++i) free(H[i]);
   free(H);
}

/* Divide-and-conquer function for computing H(k)v. This function does
   not require the matrix H(k). Only n = 2^k and v are passed as arguments */
vector HadMatVecProd ( vector v, unsigned int n )
{
   vector w, w1, w2;
   unsigned int i;

   if (n == 1) {
      w = (double *)malloc(sizeof(double));
      w[0] = v[0];
   } else {
      /* Split v in two equal halves v1 and v2. Recursively compute
         H(k-1)v1 and H(k-1)v2 */
      w1 = HadMatVecProd(v, n >> 1);
      w2 = HadMatVecProd(v + (n >> 1), n >> 1);

      w = (double *)malloc(n * sizeof(double));

      /* H(k)v = ( H(k-1)v1 + H(k-1)v2, H(k-1)v1 - H(k-1)v2 ) */
      for (i = 0; i < (n >> 1); ++i) {
         w[i] = w1[i] + w2[i]; w[(n >> 1) + i] = w1[i] - w2[i];
       }

       /* Clean local dynamic memory */
      free(w1); free(w2);
   }
   return w;
}

int main ( int argc, char *argv[] )
{
   unsigned int k, n;
   matrix H;
   vector v, w;

   /* Seed the random number generator with current time, so that we can get
      different random vectors v during different runs of the program. */
   srand((unsigned int)time(NULL));

   if (argc > 1) k = atoi(argv[1]);
   else {
      printf("k = "); scanf("%u", &k);
   }
   n = (1 << k);

   H = createHadMat(k);
   prnHadMat(H,k);
   v = genRandVec(n);
   printf("\nv = \n"); prnVec(v,n);
   w = matVecProd(H,v,n);
   printf("\n+++ Method 1\nHv = \n"); prnVec(w,n);
   free(w); freeMat(H,n);
   w = HadMatVecProd(v,n);
   printf("\n+++ Method 2\nHv = \n"); prnVec(w,n);
   free(w);

   exit(0);
}