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

void gs(double ** v, int M, int N, double * norm);
double scalar (double *a, double *b, int N);
double normalize(double *v, int N);
void product(double ** A, double ** B, int M, int N, double ** C);



int main(int argc, char ** argv) {
  double t, T, TRANS, deltat;
  double x, y, z, x0, y0, z0, x1, y1, z1;
  double s, b, r;
  double interval; // intervallo di calcolo
  
  
  s = 10;
  b = 8./3;
  r = 28;
  deltat = 0.001;
  interval = 1*deltat;
  x0 = .2;
  y0 = .5;
  z0 = .3;
  T = 1000;
  TRANS=1;
  
  int N;
  int M;
  int i, j;
  
  double ** v;
  double ** v1;
  double ** v2;
  double ** A;
  double * norm;
  double * slnorm;
  int slcount=0;
  
  srand48(123345);

  N = 3;
  M = 3;
  v = calloc(M, sizeof(double *));
  v1 = calloc(M, sizeof(double *));
  A = calloc(N, sizeof(double *));
  
  for (i=0; i<M; i++) {
    v[i] = calloc(N, sizeof(double));
    v1[i] = calloc(N, sizeof(double));
  }
  
  for (i=0; i<N; i++) {
    A[i] = calloc(N, sizeof(double));
  }
 
  
  norm = calloc(M, sizeof(double));
  slnorm = calloc(M, sizeof(double));

  A[0][0] = 1-s*deltat;
  A[0][1] = s*deltat;
  A[0][2] = 0;
  A[1][0] = 0; // (r-z)*deltat
  A[1][1] = 1-deltat;
  A[1][2] = 0; // -x*deltat
  A[2][0] = 0; // y*deltat
  A[2][1] = 0; // x*deltat
  A[2][2] = 1-b*deltat;
 
  for (i=0; i<M; i++) {
    for (j=0; j<N; j++) {
      v[i][j]= 2*drand48()-1;
    }
  }
  gs(v, M, N, norm);
  
  s = 10;
  b = 8./3;
  r = 28;
  deltat = 0.001;
  
  x = x0;
  y = y0;
  z = z0;
  for (t=0; t<T; 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;
        
    A[1][0] = (r-z)*deltat;
    A[1][2] = -x*deltat;
    A[2][0] = y*deltat;
    A[2][1] = x*deltat;

    product(A, v, M, N, v1);
    v2=v1;
    v1=v;
    v=v2;
    gs(v, M, N, norm);
    
    if (t>TRANS ){ //&& fmod(t, interval) < deltat) {
      printf("%f %f %f (%lf %lf %lf)\n", x, y, z, slnorm[0], slnorm[1], slnorm[2]);
      for (i=0; i<M; i++) {
        slnorm[i] += log(norm[i]);
      }
      slcount ++;
    }
  }
  printf("#------------ count = %d\n", slcount);
  for (i=0; i<M; i++) {
    printf("#lyap[%d]=%lf\n", i, slnorm[i]/slcount/deltat);
  }

}



void product(double ** A, double ** B, int M, int N, double ** C) {
  int i, j;
  
  for (i=0; i<M; i++) {
    for (j=0; j<N; j++) {
      C[i][j] = scalar(B[i], A[j], N);
    }
  }
}
      
    

void gs(double ** v, int M, int N, double * norm) {
  int i, j, k;
  double s;
  
  for (i=0; i<M; i++) {
    for (j=0; j<i; j++) {
      s = scalar(v[i], v[j], N);
      for (k=0; k<N; k++) {
        v[i][k] -= s*v[j][k];
      }
    }      
    norm[i] = normalize(v[i], N);
  }
}

double scalar (double *a, double *b, int N){
  int i;
  double s;
  
  s=0;
  for (i=0; i<N; i++) {
    s += a[i]*b[i];
  }
  return (s);
}

double normalize(double *v, int N) {
  int i; 
  double norm;
  
  norm = sqrt(scalar(v, v, N));
  for (i=0; i<N; i++) {
    v[i] /= norm;
  }
  return(norm);
}