/* *
 * Copyright (c) 2014, James S. Plank and Kevin Greenan
 * All rights reserved.
 *
 * Jerasure - A C/C++ Library for a Variety of Reed-Solomon and RAID-6 Erasure
 * Coding Techniques
 *
 * Revision 2.0: Galois Field backend now links to GF-Complete
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 *  - Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 *  - Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in
 *    the documentation and/or other materials provided with the
 *    distribution.
 *
 *  - Neither the name of the University of Tennessee nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
 * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY
 * WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

/* Jerasure's authors:

   Revision 2.x - 2014: James S. Plank and Kevin M. Greenan.
   Revision 1.2 - 2008: James S. Plank, Scott Simmerman and Catherine D. Schuman.
   Revision 1.0 - 2007: James S. Plank.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "jerasure.h"
#include <sys/time.h>


typedef unsigned long mylong;
#define LLUI (long long unsigned int)

double TimeStart()
{
  struct timeval tstart;
  gettimeofday(&tstart,0);
  return( (double) (tstart.tv_sec + tstart.tv_usec*1e-6) );
}

double TimeStop(double t)
{
  struct timeval tend;

  gettimeofday(&tend,0);
  t = (double) (tend.tv_sec + tend.tv_usec*1e-6) - t;
  return (t);
}

void display(mylong *mat, int r, int c) {
 for(int i=0;i<r;i++) {
    for(int j=0;j<c;j++) {
      printf("%016llu ",LLUI mat[i*c+j]);
    }
    printf("\n");
  }
 printf("\n");
}

mylong *matrix_multiply(gf_t *gf, mylong *m1, mylong *m2, int r1, int c1, int r2, int c2, int w)
{
  mylong *product;
  int i, j, k;

  product = (mylong *) malloc(sizeof(mylong)*r1*c2);
  for (i = 0; i < r1*c2; i++) product[i] = 0;

  for (i = 0; i < r1; i++) {
    for (j = 0; j < c2; j++) {
      for (k = 0; k < r2; k++) {
        product[i*c2+j] ^= gf->multiply.w64(gf,m1[i*c1+k], m2[k*c2+j]);
      }
    }
  }
  return product;
}



int invert_matrix(gf_t *gf, mylong *mat, mylong *inv, int rows)
{
  int cols, i, j, k, x, rs2;
  int row_start;
  mylong tmp, inverse;
 
  cols = rows;

  k = 0;
  for (i = 0; i < rows; i++) {
    for (j = 0; j < cols; j++) {
      inv[k] = (i == j) ? 1 : 0;
      k++;
    }
  }
//  display(inv, rows, rows);
//  printf("\n");
  
  /* First -- convert into upper triangular  */
  for (i = 0; i < cols; i++) {
    row_start = cols*i;

    /* Swap rows if we ave a zero i,i element.  If we can't swap, then the
       matrix was not invertible  */

    if (mat[row_start+i] == 0) {
      for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
      if (j == rows) return -1;
      rs2 = j*cols;
      for (k = 0; k < cols; k++) {
        tmp = mat[row_start+k];
        mat[row_start+k] = mat[rs2+k];
        mat[rs2+k] = tmp;
        tmp = inv[row_start+k];
        inv[row_start+k] = inv[rs2+k];
        inv[rs2+k] = tmp;
      }
    }
 
    /* Multiply the row by 1/element i,i  */
    tmp = mat[row_start+i];
    if (tmp != 1) {
      inverse = gf->divide.w64(gf,1, tmp);
      for (j = 0; j < cols; j++) {
        mat[row_start+j] = gf->multiply.w64(gf,mat[row_start+j], inverse);
        inv[row_start+j] = gf->multiply.w64(gf,inv[row_start+j], inverse);
      }
    }

    /* Now for each j>i, add A_ji*Ai to Aj  */
    k = row_start+i;
    for (j = i+1; j != cols; j++) {
      k += cols;
      if (mat[k] != 0) {
        if (mat[k] == 1) {
          rs2 = cols*j;
          for (x = 0; x < cols; x++) {
            mat[rs2+x] ^= mat[row_start+x];
            inv[rs2+x] ^= inv[row_start+x];
          }
        } else {
          tmp = mat[k];
          rs2 = cols*j;
          for (x = 0; x < cols; x++) {
            mat[rs2+x] ^= gf->multiply.w64(gf,tmp, mat[row_start+x]);
            inv[rs2+x] ^= gf->multiply.w64(gf,tmp, inv[row_start+x]);
          }
        }
      }
    }
  }

  /* Now the matrix is upper triangular.  Start at the top and multiply down  */

  for (i = rows-1; i >= 0; i--) {
    row_start = i*cols;
    for (j = 0; j < i; j++) {
      rs2 = j*cols;
      if (mat[rs2+i] != 0) {
        tmp = mat[rs2+i];
        mat[rs2+i] = 0;
        for (k = 0; k < cols; k++) {
          inv[rs2+k] ^= gf->multiply.w64(gf,tmp, inv[row_start+k]);
        }
      }
    }
  }

/*  printf("mat\n");
  display(mat, rows, rows);
  printf("\n");
   printf("inv\n");
  display(inv, rows, rows);
  printf("\n");
*/
  return 0;
}




int invertible_matrix(gf_t *gf, int *mat, int rows, int w)
{
  int cols, i, j, k, x, rs2;
  int row_start;
  ulong tmp, inverse;
 
  cols = rows;

  /* First -- convert into upper triangular  */
  for (i = 0; i < cols; i++) {
    row_start = cols*i;

    /* Swap rows if we ave a zero i,i element.  If we can't swap, then the 
       matrix was not invertible  */

    if (mat[row_start+i] == 0) { 
      for (j = i+1; j < rows && mat[cols*j+i] == 0; j++) ;
      if (j == rows) return 0;
      rs2 = j*cols;
      for (k = 0; k < cols; k++) {
        tmp = mat[row_start+k];
        mat[row_start+k] = mat[rs2+k];
        mat[rs2+k] = tmp;
      }
    }
 
    /* Multiply the row by 1/element i,i  */
    tmp = mat[row_start+i];
    if (tmp != 1) {
      inverse =  gf->divide.w64(gf,1, tmp);
      for (j = 0; j < cols; j++) { 
        mat[row_start+j] =  gf->multiply.w64(gf,mat[row_start+j], inverse);
      }
    }

    /* Now for each j>i, add A_ji*Ai to Aj  */
    k = row_start+i;
    for (j = i+1; j != cols; j++) {
      k += cols;
      if (mat[k] != 0) {
        if (mat[k] == 1) {
          rs2 = cols*j;
          for (x = 0; x < cols; x++) {
            mat[rs2+x] ^= mat[row_start+x];
          }
        } else {
          tmp = mat[k];
          rs2 = cols*j;
          for (x = 0; x < cols; x++) {
            mat[rs2+x] ^=  gf->multiply.w64(gf,tmp,mat[row_start+x]);
          }
        }
      }
    }
  }
  return 1;
}



int main(int argc, char **argv)
{
  unsigned int  w;
  int invert;
  mylong *matS;
  mylong *matG;
  mylong *matC;
  mylong *inverse;
  mylong *identity;


  int size=500000000;
  int t=4;
  int n=8;
  int len=size/t;

  w=64;
  gf_t gf;
  gf_init_easy(&gf, w);  
  
  matS = malloc(sizeof(mylong)*t*len);
  matG = malloc(sizeof(mylong)*t*n);

 

  srand48(time(0));
  for(int i=0;i<t;i++) {
    for(int j=0;j<len;j++) {
      matS[i*len+j]=lrand48()<<32|lrand48();
    }
  }

  for(int i=0;i<n;i++) { 
    for(int j=0;j<t;j++) {
      matG[i*t+j]=lrand48()<<32|lrand48();
    }
  }



//  printf("Matrix S:\n");
//  display(matS, t, len);

  printf("Matrix G:\n");
  display(matG, n, t);

  double tt=TimeStart();
  matC=matrix_multiply(&gf, matG, matS, n, t, t, len, w);
  printf("1st multiply %f\n",TimeStop(tt));
//  printf("Matrix C:\n");
//  display(matC, n, t);

  ulong *matCs = malloc(sizeof(mylong)*t*len);
  ulong *matGs = malloc(sizeof(mylong)*t*t);

  for(int i=0;i<t;i++) {
    for(int j=0;j<len;j++) {
      matCs[i*len+j]=matC[i*len+j];
    }
    for(int j=0;j<t;j++) {
      matGs[i*t+j]=matG[i*t+j];
    }
  }
    
  printf("Matrix Gs:\n");
  display(matGs, t, t);

//  printf("Matrix Cs:\n");
//  display(matCs, t, len);
  
  ulong* matGs_copy = malloc(sizeof(mylong)* t*t);

  //WARNINT invert changes the matrix
//  invert = invertible_matrix(&gf,matGs, t, w);
//  printf("\nInvertible Gs: %s\n", (invert == 1) ? "Yes" : "No");

  memcpy(matGs_copy, matGs, sizeof(mylong)*t*t);
  ulong *matInvGs = malloc(sizeof(mylong)*t*t);
  invert_matrix(&gf, matGs_copy, matInvGs, t);

  /*
  WARNING   this changes matGs
  identity=matrix_multiply(&gf, matInvGs, matGs, t, t, t, t, w);
  printf("Identity:\n");
  display(identity, t, t);
  */

  
  //printf("Matrix Gs:\n");
  //display(matGs, t, t);
  
  ulong *matS2 = malloc(sizeof(mylong)*t*len);
  tt=TimeStart();
  matS2=matrix_multiply(&gf, matInvGs, matCs, t, t, t, len, w);
  printf("2nd multiply %f\n",TimeStop(tt));

  //printf("Matrix S2:\n");
  //display(matS2, t, len);
  
  int equal=1;
  for(int i=0;i<t;i++) {
    for(int j=0;j<len;j++) {
      equal=matS[i*len+j]==matS2[i*len+j];
//      if(!equal)
//	printf("ARGH\n");
    }
  }
  if(equal)
    printf("EQUAL !!!\n");
  
  return 0;
}