#include<iostream>
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
using namespace std;

float** makematrix(int n, int m);
void freematrix(int n, float **a);


float powermethod(float** A, int n, float tolerance, long* idum, int k)
{
  //our current and previous guess for the eigenvalue
  float guess1 = 1, guess2 = 10;
  //our guess for the eigenvector
  float *eigenvec1 = (float*)calloc(n, sizeof(float));
  //vectors to work with
  float *eigenvec2 = (float*)calloc(n, sizeof(float));
  float *vec3 = (float*)calloc(n, sizeof(float));
  //initialize our first guess
  for(int i = 0; i < n; i++)
    eigenvec1[i] = (float)rand()/RAND_MAX;
  
  while(guess1-guess2 <= -tolerance || guess1-guess2 >= tolerance)
	{
	  
		float c = 0, d = 0, normconst = 0,test = 0;
		// set the value of c and zero out vec3
		for(i = 0; i < n; i++)
		{
		  c += eigenvec1[i];
		  vec3[i] = 0;
		}
		c = c/n;
		//set the value of d
		for(i = 0; i < n/2; i++)
			d += eigenvec1[i];
		for(i = n/2; i < n; i++)
			d -= eigenvec1[i];
		d = d/n;
		//set the value of eigenvec2
		for(i = 0; i < n/2; i++)
			eigenvec2[i] = eigenvec1[i]-c-d;
		for(i = n/2; i < n; i++)
			eigenvec2[i] = eigenvec1[i]-c+d;
		//normalize eigenvec2 with L2 norm
		for(i = 0; i < n; i++)
			normconst += eigenvec2[i]*eigenvec2[i];
		normconst = sqrt(normconst);
		for(i = 0; i < n; i++)
			eigenvec2[i] = eigenvec2[i]/normconst;
		//compute A(A(eigenvec2)), place results in eigenvec1.
		for(i = 0; i < n; i++)
			for(int j = 0; j < k; j++)
			    vec3[i] += eigenvec2[(int)A[i][j]];
		for(i = 0; i < n; i++)
			eigenvec1[i] = 0;
		for(i = 0; i < n; i++)
			for(int j = 0; j < k; j++)
			    eigenvec1[i] += vec3[(int)A[i][j]];
		guess2 = guess1;
		guess1 = 0;
		//compute eigenvalue guess that corresponds to our
		//eigenvector guess
		for(i = 0; i < n; i++)
			guess1 += eigenvec2[i]*eigenvec1[i];
		guess1 = sqrt(guess1);
   }
  free(eigenvec1);
  free(eigenvec2);
  free(vec3);
  return guess1;
}