#include"lib.h"

#define tmat uint8_t
#define N_T 4

struct node { int nb; struct node *t[2]; };
typedef struct node node;

mpf_t theta, tmp, tmp2;
mpz_t entier;
node root;
tmat **mat;
int m[1000], c[1000], i, j, sp, k, sm = 0, lb = 15, bidon[N_T], bound[N_T+1];
double *y, *v, t, ftmp;

double produit_scal_d(double *a, double *b, int it){
  double res = *a**b;
  for(; --it > 0; ) res += a[it]*b[it];
  return res;
}//*/

double produit_scal_d_tmat(double *a, tmat *b, int it){
  double res = 0;
  for(; --it >= 0; ) if(b[it>>3]&(1<<(it&7))) res += a[it];
  return res;
}//*/

double norm(double *a, int it){ return sqrt(produit_scal_d(a, a, it)); }

void produit_scal_mpf(mpf_t res, mpf_t *a, mpf_t *b, int it){
  mpf_mul(res, *a, *b);
  for(; --it > 0; ){
    mpf_mul(tmp, a[it], b[it]);
    mpf_add(res, res, tmp);
  }
}//*/

void produit_scal_mpf_tmat(mpf_t res, mpf_t *a, tmat *b, int it){
  mpf_set_ui(res, 0);
  for(; --it >= 0; ) if(b[it>>3]&(1<<(it&7))) mpf_add(res, res, a[it]);
}//*/

mpf_t *alloc_vector_mpf(int it){
  mpf_t *tmp = (mpf_t*)malloc(it*sizeof(mpf_t));
  for(; it > 0; ) mpf_init(tmp[--it]);
  return tmp;
}

void mpf_v_set_prec(mpf_t* v, int it, int prec){
  for(; it > 0; ) mpf_set_prec(v[--it], prec);
}

void display(){
  mpf_out_str(stdout, 10, 30, theta);
  mpf_mul_2exp(tmp, theta, 1000*(k-sp));
  mpz_set_f(entier, tmp);
  mpz_root(entier, entier, k-sp);
  mpf_set_z(tmp, entier);
  mpf_div_2exp(tmp, tmp, 1000);
  printf(" ");
  mpf_out_str(stdout, 10, 30, tmp);
  //mpf_pow_ui(theta, tmp, k-sp);
  //printf(" ");
  //mpf_out_str(stdout, 10, 30, theta);
  printf(" nb_it=%d PREC=%ld", j, mpf_get_default_prec()); date("");
}

int sq(){
  int xx, jj, kk = 1;
  for(; (xx = 2*kk) <= i+1; kk++){
    for(jj = i-kk; (jj >= 0) && (m[jj] == m[jj+kk]); jj--);
    if(i-jj >= xx) return 1;
  }
  return 0;
}

int sq2(int *m, int i){
  int xx, jj, kk = 3;
  for(; (xx = 2*kk) <= i+1; kk++){
    for(jj = i-kk; (jj >= 0) && (m[jj] == m[jj+kk]); jj--);
    if(i-jj >= xx) return 1;
  }
  return 0;
}

int code(int *m){
  int j = 0;
  node *cur = &root;
  for(; j < sp; j++) if(cur->t[m[j]]) cur = cur->t[m[j]]; else return -1;
  return cur->nb;
}

void add(int *m){
  node *cur = &root;
  for(j = 0; j < sp; j++) if(!cur->t[m[j]]){
    cur->t[m[j]] = (node*)malloc(sizeof(node));
    cur = cur->t[m[j]];
    cur->t[0] = cur->t[1] = 0;
  }else cur = cur->t[m[j]];
  cur->nb = 0;
}

void num(node *nd, int lvl){
  if(!nd) return;
  if(lvl == sp && (nd->t[0] || nd->t[1])) printf("alouf\n");
  if(lvl == sp) nd->nb = sm++; else{
    num(nd->t[0], lvl+1);
    num(nd->t[1], lvl+1);
  }
}

void make_tree(){
  printf("sp=%d k=%d lb=%d tree construction", sp, k, lb); date("");
  root.t[0] = root.t[1] = 0;
  *c = *m = m[i = 2] = 0; m[1] = 1;
  for(; ; ) if(sq()){
    there:while(m[i] == 2) if(--i < 2) goto end;
    m[i]++;
    c[i-2] = (m[i] != m[i-2]);
  }else if(i == sp+lb+1 && !code(c+lb)) goto there;
  else if(i == sp+2*lb){
    add(c+lb);
    i = sp+lb;
   goto there;
  }else{
    m[++i] = 0;
    c[i-2] = (m[i] != m[i-2]);
  }
  end:num(&root, 0);
}

void *morceau(void *ptr){
  int feve = *((int *) ptr), c_i, c_o, i, m[1000], c[1000];

  *c = *m = m[2] = i = 0; m[1] = 1;
  for(; ; ) if(sq2(m, i+2)){
    here:while(c[i]) if(!i--) return NULL;
    c[i] = 1;
    m[i+2] = 3-m[i]-m[i+1];
  }else if(i == k){
    if((c_o = code(c+k-sp)) < 0) goto here;
    if(mat[c_o][c_i>>3]&(1<<(c_i&7))) {printf("overflow : %d\n", mat[c_o][c_i]); exit(0);}
    mat[c_o][c_i>>3] |= 1<<(c_i&7);
    i = k-1;
    goto here;
  }else{
    if(i == sp-1)
      if((c_i = code(c)) < bound[feve]) goto here;
      else if(c_i >= bound[feve+1]) return NULL;
    i++;
    c[i] = !c[i-1];
    m[i+2] = (c[i] ? 3-m[i]-m[i+1] : m[i]);
  }
}

void make_matrix(){
  pthread_t thread_id[N_T];

  printf("sp=%d sm=%d k=%d matrix construction", sp, sm, k); date("");
  y = (double*)malloc(sm*sizeof(double)); v = (double*)malloc(sm*sizeof(double));
  mat = (tmat**)malloc(sm*sizeof(tmat*));
  for(i = 0; i < sm; i++) if(!(mat[i] = (tmat*)malloc((sm/8+1)*sizeof(tmat)))) // 
    {printf("malloc failed sm=%d i=%d\n", sm, i); date(""); exit(0);}
  printf("allocation memoire matrice OK\n");
  for(i = 0; i < sm; i++) for(j = 0; j <= sm/8; j++) mat[i][j] = 0;
  for(i = 0; i < N_T; i++) bidon[i] = i;
  for(i = 0; i < N_T; i++) bound[i] = i*sm/N_T;
  bound[N_T] = sm;
  for(i = 0; i < N_T; i++) pthread_create(&thread_id[i], NULL, morceau, &bidon[i]);
  for(i = 0; i < N_T; i++) pthread_join(thread_id[i], NULL);
}

void power_method(){
  date("first iterations (double)");
  for(i = 0; i < sm; i++) y[i] = 1;
  for(j = 0; j < 50; j++){ // j < 50 si ca ne converge pas
    ftmp = norm(y, sm);
    for(i = 0; i < sm; i++) v[i] = y[i]/ftmp;
    for(i = 0; i < sm; i++) y[i] = produit_scal_d_tmat(v, mat[i], sm);
    t = produit_scal_d(y, v, sm);
    for(i = 0; i < sm; i++) v[i] = y[i]-t*v[i];
    ftmp = norm(v, sm);
    printf("%.30f nb_it=%d\n", t, j);
    if(ftmp*((uint64_t)1<<45) < t) break;
  }
  date("last iterations (mpf_t)");
  //mpf_init(theta); mpf_init(tmp); mpf_init(tmp2); mpz_init(entier);
  mpf_init2(theta, 256); mpf_init2(tmp, 256); mpf_init2(tmp2, 256); mpz_init(entier);
  mpf_t *vect = alloc_vector_mpf(sm), *unit = alloc_vector_mpf(sm);
  for(i = 0; i < sm; i++) mpf_set_d(vect[i], y[i]); //mpf_set_ui(vect[i], 1);
  for(j = 0; ; j++){
    produit_scal_mpf(theta, vect, vect, sm);
    mpf_sqrt(theta, theta);
    for(i = 0; i < sm; i++) mpf_div(unit[i], vect[i], theta);
    for(i = 0; i < sm; i++) produit_scal_mpf_tmat(vect[i], unit, mat[i], sm);
    produit_scal_mpf(theta, unit, vect, sm);
    mpf_set_ui(tmp, 0);
    for(i = 0; i < sm; i++){
      mpf_mul(tmp2, theta, unit[i]);
      mpf_sub(tmp2, vect[i], tmp2);
      mpf_mul(tmp2, tmp2, tmp2);
      mpf_add(tmp, tmp, tmp2);
    }
    mpf_sqrt(tmp, tmp);
    mpf_mul_2exp(tmp2, tmp, 100); // calcul de espilon
    display();
    //if(mpf_cmp(tmp2, theta) < 0) return;
    if(mpf_cmp(tmp2, theta) < 0)
      if(mpf_get_default_prec() == 256) return;
      else{
        mpf_v_set_prec(vect, sm, 256);
        mpf_v_set_prec(unit, sm, 256);
        mpf_set_default_prec(256);
      }
  }
}

void usage(){
  printf("Usage: matrix k l\n");
  exit(0);
}

int main(int ac, char **av){
  if(ac != 3) usage();
  sp = atoi(av[1]);
  k  = atoi(av[2])+sp;
  if(sp >= k) usage();
  mpf_set_default_prec(128);
  make_tree();
  make_matrix();
  power_method();
  //iterations(mat, sm, sp, k-sp);
}