#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define SIZE 1001
/* largest eigenvalue Lanczos, first Householder to triadiag */
/* example matrix=hard-square transfer matrix */
int i,j,k,l,n;
float A[SIZE][SIZE],V[SIZE],tmp[SIZE], norm, lambda;
float M[SIZE][SIZE], TARG[SIZE][SIZE];
int N,x,y,Z,Dim, CYCLES;
float sum, v[SIZE], tmp[SIZE];
void get_bin(int N, int *n, int size);
float get_HS(int row, int col, int size);
int get_dim(int a);
void Tridiag (float Aptr[][SIZE], float targ[][SIZE], int N);
main(int argc, char *argv[]){
    if(argc!=3){
	printf("./tridiag N CYCLES\n");
	exit(1);
    }

    Dim=atoi(argv[1]);
    CYCLES=atoi(argv[2]);
    Z=get_dim(Dim);
    for(x=0;x<Z;x++){
	for(y=0;y<Z;y++){
	    M[x][y]=get_HS(x,y,Z);
	}
       
    }
    Tridiag(M,A,Z);


srand(131);
 
 for(i=0;i<Z;i++){
     v[i]=1.0/sqrt((float)Z);
 }

 for(j=0;j<CYCLES;j++){
     /* ten cycles ? */
    sum=0.0;
    tmp[0]=A[0][0]*v[0]+A[0][1]*v[1];
    sum=sum+tmp[0]*tmp[0];
    for(i=1;i<Z-1;i++){
	tmp[i]=A[i][i-1]*v[i-1]+A[i][i]*v[i]+A[i][i+1]*v[i+1];
	sum=sum+tmp[i]*tmp[i];
    }
    tmp[Z-1]=A[Z-1][Z-2]*v[Z-2]+A[Z-1][Z-1]*v[Z-1];
    sum=sum+tmp[Z-1]*tmp[Z-1];
    sum=sqrt(sum);

    for(i=0;i<Z;i++)
	v[i]=tmp[i]/sum;

    printf("#");
    fflush(NULL);
 }/* end of ten cycles */



 /*for(x=0;x<Z;x++){
	for(y=0;y<Z;y++){
	    printf("%d ", (int)M[x][y]);
	}

	printf("\n");
    }

    printf("\n"); */
	/* for(x=0;x<Z;x++){
	for(y=0;y<Z;y++){
	   	printf("%f ", A[x][y]);
	}

	printf("\n");
	} */
/* one more time */

    tmp[0]=A[0][0]*v[0]+A[0][1]*v[1];
    for(i=1;i<Z-1;i++){
	tmp[i]=A[i][i-1]*v[i-1]+A[i][i]*v[i]+A[i][i+1]*v[i+1];
    }
    tmp[Z-1]=A[Z-1][Z-2]*v[Z-2]+A[Z-1][Z-1]*v[Z-1];
    if(v[0]!=0.0)lambda=tmp[0]/v[0];
    printf("\nlambda_{max}=%f\n", lambda); 
}

void get_bin(int N, int *n, int size){
    int i, tmp;
    tmp=N;
    for(i=0;i<size;i++){
	n[i]=tmp-2*(tmp/2);
	tmp=(tmp-n[i])/2;
    }
}



float get_HS(int row, int col, int size){
    int n[SIZE], m[SIZE], factor1, factor2;
    int i,j;
    get_bin(row,n,size);
    get_bin(col,m,size);

    factor1=1,factor2=1;
    
    for(i=0;i<size-1;i++){
	factor1=factor1*(1-n[i]*n[i+1])*(1-m[i]*m[i+1]);
    }
    for(i=0;i<size;i++){
        factor2=factor2*(1-n[i]*m[i]);
    }

    if(factor2*factor1==0)
	return(0.0);
    else
	return(1.0);

}


int get_dim(int a){
    int i,tmp;
    tmp=1;
    for(i=0;i<a;i++)
	tmp=tmp*2;
    return(tmp);
}


void 
Tridiag (float Aptr[][SIZE], float targ[][SIZE], int N)
{
  /* Routine (successfully) reduces a symmetric matrix to tridiagonal */
  /* by the Householder method                                        */
  float B[SIZE][SIZE], q[SIZE],p[SIZE],u[SIZE],x[SIZE],one[SIZE][SIZE];
  int col;
  float H ,K, Lx;
  int m, n;
/* copy A[][] into B[][] as a local copy, point Bptr to it */
  /* the original matrix is not lost                       */
  for (n = 0; n < N; n++)
    {
      for (m = 0; m < N; m++)
	{
          if(n==m)
	      one[n][m]=1.0;
	  else
	      one[n][m]=0.0;
	  B[n][m] = Aptr[n][m];
	}

    }

 /* This is where over columns begins */
 for(col=0;col<N-2;col++){
for(n=0;n<=col;n++)
x[n]=0.0;
Lx=0.0;
 for(n=col+1;n<N;n++){



  x[n]=B[n][col];    /* construct x[], H */
  Lx=Lx+x[n]*x[n];}
for(n=0;n<N;n++)
u[n]=x[n];
u[col+1]=u[col+1]-sqrt(Lx);
H=Lx-sqrt(Lx)*x[col+1];
if(Lx==0.0|| H==0.0)break;


 for(n=0;n<N;n++){ /* construct p[n] */
p[n]=0.0;
   for(m=0;m<N;m++){
     p[n]=p[n]+B[n][m]*u[m];}
p[n]=p[n]/H;
 }

K=0.0;
 for(n=0;n<N;n++){ /* construct K */
   K=K+u[n]*p[n];}
K=K/(2.0*H);

 for(n=0;n<N;n++){
   q[n]=p[n]-K*u[n];} /* construct q[] */

 for(n=0;n<N;n++){
   for(m=0;m<N;m++){
     B[n][m]=B[n][m]-q[n]*u[m]-u[n]*q[m];}} /* transforms matrix */


 } /* terminates looping over columns */



 for(n=0;n<N;n++){
   for(m=0;m<N;m++){
       targ[n][m]=B[n][m];
   }
 }
}