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

extern double gasdev(long *idum);

float** makegraph(int n, int k);
long idum = 0;
void freematrix(int n, float** a);
float powermethod(float ** A, int n, float a, long *idum, int k);

void main()
{
	int n = 6; int k = 7;
	long idum = rand();
	char filename[20] = "600data.txt";
	for(n = 600; n < 1001; n+=100)
	{
		ofstream fout;
		int regular = 0, ramanujan = 0;
		// counters to keep track of number of regular and ramanujan
		fout.open(filename);
		fout<<n<<"  nodes"<<endl;
		for(int z = 0; z < 5000; z++)
		{
			float ** A = makegraph(n, k);
			float eigen = powermethod(A, n, .00001, &idum, k);
			fout<<eigen<<'\t';
			if(eigen <=2*sqrt(6))
				ramanujan++;
			regular++;
			fout<<regular<<'\t'<<ramanujan<<endl;
			freematrix(n, A);
		}
		if(filename[0] == '9')
			filename[0] = 'a';
		filename[0]++;
		fout.close();
	}
}

float** makematrix(int n, int k)
//return an nxk matrix of floats
{
	float** graph;
	graph = (float**)calloc(n, sizeof(int*));
	for(int i = 0; i < n; i++)
		graph[i] = (float*)calloc(k,sizeof(float));
	return graph;
}

int compare(const void* e1, const void* e2)
//comparison for call to qsort in makegraph
{
	int* v1 = *(int**)(e1), *v2 = *(int**)(e2);
	if(v1[0] < v2[0]) return -1;
	if(v1[0] == v2[0]) return 0;
	return 1;
}

int comparefloat(const void* e1, const void* e2)
{
	float v1 = *(float*)(e1), v2 = *(float*)(e2);
	if(v1 < v2) return 1;
	if(v1 == v2) return 0;
	return -1;
}

void freematrix(int n, float** a)
//free a matrix defined by makematrix()
{
	for(int i = 0; i < n; i++)
		free(a[i]);
	free(a);
}

float** makegraph(int n, int k)
//makes a random k-regular bipartite graph of size n
{
	int i, j; float** graph;
	if(n%2 != 0)
		fprintf(stderr, "n needs to be EVEN\n");
	graph = makematrix(n, n);
	float** b = makematrix(n/2, 2);
	bool flag = 0;
	while(flag == 0)
	{
		for(j = 0; j < k; j++)
		{
			for(i = 0; i < n/2; i++)
			{
				b[i][0] = rand();
				b[i][1] = n/2+i;
				//assign each node label between n/2 and n a random number
				//stick everything in b
			}
			qsort((void*)b, n/2, sizeof(int*), compare);
			// sort b based on the rand's.node labels 
			// along for the ride" are put into random order
			for(i = 0; i < n/2; i++)
				// insert these random edges into the graph
				// represented by an adjacency list.
			{
				graph[i][j] = b[i][1];
				graph[(int)b[i][1]][j] = i;
			}
		}
		flag = 1;
// remove the comments below if you only want single bonds, using the throw away method
	/*	for(i = 0; i < n; i++)
			for(j = 0; j < k-1; j++)
			{
				float test = graph[i][j];
				for(int l = j+1; l < k; l++)
					if(graph[i][l] == test)
						flag = 0;
			}
		if(flag == 0)
			for(i = 0; i < n;i++)
				for(j = 0; j < k; j++)
					graph[i][j] = 0;*/

	}
		freematrix(n/2, b);
		return graph;
}


float** makegraph2(int n, int k)
//makes a random k-regular bipartite graph of size n
{
	int i, j; float** graph;
	if(n%2 != 0)
		fprintf(stderr, "n needs to be EVEN\n");
	graph = makematrix(n, n);
	float** b = makematrix(n/2, 2);
	int* nodes = (int*)calloc(n/2, sizeof(int));
	for(j = 0; j < k; j++)
	{
		int bad;
		for(i = 0; i < n/2; i++)
			nodes[i] = n/2+i;
		for(i = 0; i < n/2; i++)
		{
			bad = 1;
			while(bad == 1)
			{
				bad = 0;
				int place = (int)((float)rand()/RAND_MAX*(n/2-i-1));
				//cout<<place<<endl;
				for(int l = 0; l < j; l++)
					if(nodes[place] == graph[i][l])
					{
						if(i < n/2-j)
							bad = 1;
						else bad = 3;
						continue;
					}
				if(bad == 3)
					break;
				if(bad == 1)
					continue;
				b[i][1] = nodes[place];
				for(l = place; l < n/2-1; l++)
					nodes[l] = nodes[l+1];
			}
			if(bad == 3)
			{j--; break;}
		}
		if(bad != 3)
			for(i = 0; i < n/2; i++)
				// insert these random edges into the graph
				// represented by an adjacency list.
			{
				graph[i][j] = b[i][1];
				graph[(int)b[i][1]][j] = i;
			}
	}
	free(nodes);
	freematrix(n/2, b);
	return graph;
}