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


#define NUMITER 5000
#define odd(x) ((x) & 1)

#define twopower80 (((float) (1LL<<40))*((float) (1LL<<40)))
#define twopowerminus100 1/(((float) (1LL<<50))*((float) (1LL<<50)))


/*no allibero*/

void motzkinpath2setofsets(int *p, int *outset, int n){
  
  int ind, i,j, level;  
  int *size_of_set; 
  int **set;

  size_of_set = (int *)malloc(n * sizeof(int));
  if (size_of_set ==0) {printf("Error malloc size_of_set\n"); exit(1);}

  set=(int **)malloc(n * sizeof(int *));
  if (set==0) {printf("Error malloc set\n"); exit(1);}
  for (j = 0; j < n-1 ; j++){
   set[j] = (int *)malloc(n * sizeof(int));
   if (set[j]==0) {printf("Error malloc set\n"); exit(1);}
  }

  ind= 0;
  level= 0;
  
  for(i= 0; i<=n-1 ; i++){
    switch(p[i]){
    case 1: 
      level++;
      if(odd(level)){                 /*start a new set set[(level-1)/2]*/
	size_of_set[level>>1]= 0;
      }
      break;
    case -1: 
      if(odd(level)){                 /*current set is complete; store set[(level-1)/2]*/
	for(j= 0;j<size_of_set[level>>1];) outset[ind++]= set[level>>1][j++];
	outset[ind++]= 0;        
      }
      level--;
      break;
    } 
    if(odd(level)){                   /*add i+1 (els elements van de 1..n) to current set*/
      set[level>>1][size_of_set[level>>1]++]= i+1;
    }
  }
  outset[ind++]= 0;                   /*add another final terminal marker*/
}

void generate_motzkinpath( long long **stock, int n, int maxstock){

  int i, row, col, height;
 
  /* altura 0..maxstock */ 

  stock[0][n]=1;
  for (i=1 ; i<= maxstock ; i++) {
    stock[i][n-i+1]=0; 
    stock[i][n-i+2]=0;
  }
  
  for (col = n-1; col >= 0 ; col--){
  
    if (col < maxstock)  height = col;                           /* col < n- col */
    else if (col > maxstock +1 ) height = n-col;                 /* col > n- col */
    else if (col == maxstock + 1 && !odd(n))  height = n-col;
    else {                                          /* col = n-col = maxstock, highest stock a part*/
      height = maxstock -1; 
      stock[height+1][col] = stock[height][col+1]+stock[height+1][col+1];
    }

    stock[0][col] = stock[0][col+1]+ stock[1][col+1];                                                     

    for (row = 1 ; row <= height ; row++)
      stock[row][col] = stock[row-1][col+1]+stock[row][col+1]+stock[row+1][col+1]; 
  }
}

void number2motzkinpath(int np, int *p, long long **stock, int n,  int maxstock){

  int h, col, num;

  h = 0;
  num = np;

  for(col=0 ; col <= n-1; col++){ 
    if (h == 0){   
      if (stock[0][col+1]>=num) p[col] = 0;                
      else{
        p[col] = 1;
        num = num-stock[h][col+1]; 
        h++;
      }              
    }                                                                                  
    else{      
      if (stock[h-1][col+1]>=num){
        p[col]=-1;
        h--;    
      }                                                              
      else if (stock[h-1][col+1]+stock[h][col+1]>=num){
        p[col]=0; 
        num = num-stock[h-1][col+1];         
      }                                                            
      else{ 
        p[col]=1; 
        num=num-stock[h-1][col+1]-stock[h][col+1];
        h++;   
      } 
    } 
  }
}

long long motzkinpath2number(int *p, long long **stock){

  int h, i;
  long long min;
 
  h=0;
  min=1;   
  /* la i enlloc de 1 comenca per 0*/

  for(i = 0; stock[h][i] != 1; i++){ 
    if ( p[i] == 1 ){
      if (h>0) min = min + stock[h-1][i+1]+stock[h][i+1];
      else     min = min + stock[h][i+1];
      h++;
    }
    else if (p[i]==-1) h--; 
    else
      if (h>0) min = min + stock[h-1][i+1]; 
  }  
  return(min);
}
 
int *shifting_(int *p, int n){

  int i; 
  
  for(i=n-1 ;i >= 1; i--) p[i] = p[i-1];
  p[0]=0; 
  return(p);
}

int* shiftingw(int *p, int n){

  int i,level,last;
                                                                                                                       
  if (p[n-2]==1){   
    p[n-1]=0;
    for (i = n-2; i>=2; i--) p[i] = p[i-1];
    p[0]=1;
    p[1]=-1;                                                    
  }
  else{  
    last=n-2;
    level=1-p[last];   
    while (level!=1){  
      p[last+1]=p[last];
      last--;
      level=level-p[last];
    }                                                         
    if (last < n-2){ 
       p[n-1]=0;
       p[last+1]=-1;
    }                                                                         
    last--;
    level=level-p[last];  
    while (level != 0){ 
      p[last+1]=p[last];
      last--;
      level=level-p[last];
    }                                       
    if (last==0){
      p[0]=1;
      p[1]=0;
      return(p);                                                                               
    }   
    p[last+1]=-1;    
    for (i=last ;i>= 2 ; i--) p[i]=p[i-1]; 
    p[0]=1;
    p[1]=1;
  }                                                                                                   
  return(p);
}


 int w, n, j, k, first, last, level, maxstock, scaled;
 float vap, min, max;
 long long i, M;
 char ch;

 float *new, *old, *temp;
 long long *succ0, *succ1;
 int *p, *paux, *tempp, *setofsets;

 long long **stock;

#include "printnums.c"

 
int main(int argc, char **argv){

 if (argc != 2){ 
   printf("Error, w?\n");
   exit(0);
 }

 w = atoi(argv[1]);

 n = w+1;
 maxstock = n/2;  /* floor, son enters.  L'altura es maxstock+1  */ 
 
 stock=(long long **)malloc((maxstock + 1) * sizeof(long long *));
 if (stock==0) {printf("Error malloc stock\n"); exit(1);}
 for (j = 0; j < maxstock + 1 ; j++){
   stock[j] = (long long *)malloc((n+1) * sizeof(long long));
   if (stock[j]==0) {printf("Error malloc stock\n"); exit(1);}
 }

 generate_motzkinpath(stock,n,maxstock);
 M=stock[0][0];
 
 printf("n= %d, M= %lld,  long long = %d bytes, pointer = %d bytes,\n",n , M,sizeof(long long),sizeof(succ0));
 /*scanf("%c",&ch);*/

 succ0 = (long long *)malloc((M - 1) * sizeof(long long));
 if (succ0==0) {printf("Error malloc succ0\n"); exit(1);}

 succ1 = (long long *)malloc((M -1) * sizeof(long long));
 if (succ1==0) {printf("Error malloc succ1\n"); exit(1);}

 p = (int *)malloc(n * sizeof(int));
 if (p==0) {printf("Error malloc p\n"); exit(1);}

 paux = (int *)malloc(n * sizeof(int));
 if (paux==0) {printf("Error malloc paux\n"); exit(1);}

 for(i=2; i<=M; i++){

   /* states van de 0 a M-2, i succ0 sera nul quan succ0=-1 */
   /* resto -1 als indexos perque els vectors amb C van de 0..n-1*/

   number2motzkinpath(i, p, stock, n, maxstock);

   if (p[n-1]==0){
     memcpy(paux, p, sizeof(int) * n);
     succ0[i-2]=motzkinpath2number(shifting_(paux,n),stock)-2; 
   }  
   else{   
     memcpy(paux, p, sizeof(int) * n);                         /* paux = p */
     
     if (p[n-2]==0){   
       paux[n-2]=-1;    
       paux[n-1]=0;
       
       succ0[i-2]=motzkinpath2number(shifting_(paux,n),stock)-2;
     }
     else if (p[n-2]==-1){     
       paux[n-1]=0; 
       paux[n-2]=0;  
       last=n-3;
       level=2-paux[last];       
       
       while (level>1){
	 last--;
	 level=level-paux[last];
       }
       paux[last]=-1;
       
       succ0[i-2]=motzkinpath2number(shifting_(paux,n),stock)-2;       
     }
     else{
       succ0[i-2]=-1;          
     }  
   }
   if(p[0]==1){  
     if(p[n-1]==-1){             /* p[1]=1 and p[n]=-1 */     
       last=n-2;
       level=1-p[last];
       
       while (level>0){
	 last--;
	 level=level-p[last];
       }         
       if(last>0){
	 p[last]=-1;
	 first=1; 
	 level=1+p[first];             
	 
	 while(level>0){
	   first++; 
	   level=level+p[first];
	 }                                                       
	 p[first]=1;                                                                                      
       }   
       succ1[i-2]=motzkinpath2number(shiftingw(p,n),stock)-2;                                                         
     }   
     else{                         /* p[1]=1 and p[n]=0 */ 
       shifting_(p,n); 
       p[0]=1; 
       p[1]=0; 
        succ1[i-2]=motzkinpath2number(p,stock)-2;
     }
   }                                                                                                      
   else{          
      if(p[n-1]==-1) succ1[i-2]=motzkinpath2number(shiftingw(p,n),stock)-2;    /* p[1]=0 and p[n]=-1 */ 
      else{                       /* p[1]=0 and p[n]=0 */   
        shifting_(p,n);
        p[0]=1;
        p[1]=-1;
        succ1[i-2]=motzkinpath2number(p,stock)-2;
      }
   }
 } 
 printf("Successor generated, begin iteration \n ");
 
 new = (float *)malloc((M - 1) * sizeof(float));
 if (new==0) {printf("Error malloc new\n"); exit(1);}

 old = (float *)malloc((M - 1) * sizeof(float));
 if (old==0) {printf("Error malloc old\n"); exit(1);}

 for (i=0 ; i <= M-2 ; i++) old[i]=1.;  
 
 setofsets = (int *)malloc(n * sizeof(int));
 if (setofsets==0) {printf("Error malloc setofsets\n"); exit(1);}
 
 scaled=0;
 for (j=0; j < NUMITER ; j++){
   
   for(i=0; i<=M-2; i++){                                                                    
     if (succ0[i] != -1) new[i] = new[succ0[i]]+old[succ1[i]];                                              
     else new[i] = old[succ1[i]];
   } 

   if (j%10==0){

     /* find min and max */
     min = max = new[0]/old[0];
     
     for(i=1 ; i <= M-2 ; i++){
       vap = new[i]/old[i];
       if (vap < min) min = vap;
       else if (vap > max) max = vap;
     } 
     printf("iter %d: min= %0.9f, max= %0.9f\n",j,min,max);

     if (max<1.000001*min){ 
       printf("Scaled %i times\n ", scaled);
       printf("Do you want to see the (scaled) eigenvector (y/n)?\n ");
       scanf("%c",&ch);
       if (ch=='y'){
         printcheck();
       }
       exit(0);
     }

     /* reescalar si cal */
     max = min = new[0];
     for(i=1 ; i <= M-2 ; i++){
       if (max < new[i]) max = new[i];
       if (min > new[i]) min = new[i];
     } 
     if (max > twopower80){
       printf("newmin= %g, newmax= %g, ratio=%g , %f\n ", min, max, max/min, twopower80);
       for(i=0 ; i <= M-2 ; i++) new[i]*=twopowerminus100;
       scaled++;  
     }
   }
   temp=new;                        /*old=new*/
   new=old; 
   old=temp;
 }                               
exit (0);
}