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

typedef struct {
  int l; // length
  double * d; // data
} vector;

typedef struct {
  int l; // length
  double ** d; // data
} operator;


vector * new_vector(int length);
operator * new_operator(int length);
double scalar (vector * x, vector * y);
void scale (vector * x, double a);
double norm(vector * x);
double normalize(vector * x);
void combine(vector * v, double a, vector *x, double b, vector * y);
void applyoperator(vector * w, operator * A, vector * v);
void copyvector(vector * w , vector *v);
void randomize(vector * v);
void gs(vector ** v, int m, double * norm);
void printoperator(operator *A);
void printvector(vector *v);


int main(int argc, char ** argv) {
  double t, T, deltat, TRANS;
  double x, y, z, x0, y0, z0, x1, y1, z1;
  double s, b, r;
  double vis;
  int L=3; // dimension of vectors
  int m=3; // number of vectors 
  double *norm;
  double *slognorm;
  int scount;
  vector **v;
  vector **w;
  int i;
  operator * J; // jacobian

  norm = calloc(m, sizeof(double));
  slognorm = calloc(m, sizeof(double));
  scount = 0;

  J = new_operator(L);
  
  v = calloc(m, sizeof(vector *));
  for (i=0; i<m; i++) {
    v[i] = new_vector(L);
    randomize(v[i]);
  }
  
  w = calloc(m, sizeof(vector *));
  for (i=0; i<m; i++) {
    w[i] = new_vector(L);
  }
  
  gs(v, m, norm); // ornonormali

  s = 10;
  b = 8./3;
  r = 28;
  x0=0.3;
  T =100;
  TRANS=100;
  y0=1;
  z0=1;
  deltat = 0.001;
  vis = 0.01;
  
  J->d[0][0] = 1-s*deltat;
  J->d[0][1] = s*deltat;
  J->d[0][2] = 0;
  J->d[1][0] = 0;// (r-z)*deltat
  J->d[1][1] = 1-deltat;
  J->d[1][2] = 0; // -x * deltat
  J->d[2][0] = 0; // y * deltat
  J->d[2][1] = 0; // x * deltat
  J->d[2][2] = 1-b*deltat;
  
  x = x0;
  y = y0;
  z = z0;
  for (t=0; t<T+TRANS; t+= deltat) {
    x1 = x + s*(y-x)*deltat;
    y1 = y + (x*(r-z)-y)*deltat;
    z1 = z + (x*y - b*z)*deltat;
    x = x1;
    y = y1;
    z = z1;
    
    if (t > TRANS) {
      J->d[1][0] = (r-z)*deltat;
      J->d[1][2] =  -x * deltat;
      J->d[2][0] =  y * deltat;
      J->d[2][1] =  x * deltat;
      for (i=0; i<m; i++) {
        applyoperator(w[i], J, v[i]);
         copyvector(v[i], w[i]);
      }
      gs(v, m, norm);
      for (i=0; i<m; i++) {
        slognorm[i] += log(norm[i]);
      }
      scount ++;
    }
  }
  for (i=0; i<m; i++) {
    printf("lyap[%d] = %f\n", i, slognorm[i]/scount/deltat);
  }
}


vector * new_vector(int length) {
  vector * v;
  v = calloc(1, sizeof(vector));
  v->l = length;
  v->d = calloc (length, sizeof(double));
  return v;
}

operator * new_operator(int length) {
  operator * A;
  int i;
  
  A = calloc(1, sizeof(vector));
  A->l = length;
  A->d = calloc (length, sizeof(double *));
  for (i=0; i<length; i++) {
    A->d[i]= calloc (length, sizeof(double));
  }
  return A;
}

double scalar (vector * x, vector * y) {
  double s;
  int i;
  
  if (x->l != y->l) {
    fprintf(stderr, "error in scalar: vectors have different length.\n");
    exit(1);
  }
  s = 0;
  for (i=0; i< x->l; i++) {
    s += x->d[i] * y->d[i];
  }
  return s;
}

void scale (vector * x, double a) {
  int i; 
  for (i=0; i< x->l; i++) {
    x->d[i] *= a;
  }
}

double norm(vector * x) {
  double norm;
  
  norm = sqrt(scalar(x,x));
  
  return(norm);
}

double normalize(vector * x) {
  double n;
  
  n = norm(x);
  scale(x, 1./n);
  
  return(n);
}

// v = a*x+b*y
void combine(vector * v, double a, vector *x, double b, vector * y) {
  int i; 
  if (x->l != y->l || x->l != v->l) {
    fprintf(stderr, "error in combine: vectors have different length.\n");
    exit(1);
  }
  for (i=0; i< v->l; i++) {
    v->d[i] = a * x->d[i] + b * y->d[i];
  }
}


void applyoperator(vector * w, operator * A, vector * v) {
  int i, j; 
  if (v->l != w->l || v->l != A->l) {
    fprintf(stderr, "error in applyoperator: vectors have different length.\n");
    exit(1);
  }
  for (i=0; i< v->l; i++) {
    w->d[i] = 0;
    for (j=0; j< v->l; j++) {
      w->d[i] += A->d[i][j] * v->d[j];
    }
  }
}

void copyvector (vector * w , vector *v) {
  int i;
  
  for (i=0; i < v->l; i++) {
    w->d[i] = v->d[i];
  }
}

void randomize (vector * v) {
  int i;

  for (i=0; i< v->l; i++) {
    v->d[i] = drand48();
  }
}
 
void gs(vector ** v, int m, double * norm) {
  double a;
  int i, j;
  
  for (i=0; i<m; i++) {
    for (j=0; j<i ; j++) {
      a = scalar(v[i], v[j]);
      combine(v[i], 1, v[i], -a, v[j]);
    }
    norm[i] = normalize(v[i]);
  }
}
  
void printoperator(operator *A) {
  int i, j;
  for (i=0; i<A->l; i++) {
    for (j=0; j<A->l; j++) {
      printf("%f ", A->d[i][j]);
    }
    printf("\n");
  }
}

void printvector(vector *v) {
  int i; 
  for (i=0; i< v->l; i++) {
    printf("%f ", v->d[i]);
  }
  printf("\n");
}