#include "supercon.h"

#define Maxhardness 20
#define MaxTime 40
int TAG = 870619;

int myrank, numprocs;
bignum rsaN, rsaE;
bignum hint[Nxy];
int eimage[Nxy], image[Nxy];
int hardness[Nxy];

unsigned long own_data[Nxy*(3 + 1)]; int own_data_n, arr_data[129], displs[129];
unsigned long new_data[Nxy*3*(64 + 1)];
int list[Nxy * 2], list_n, last_list_n;
int hash[Nxy * 2];


//********************************************************
void signedbsub(bignum* a, bignum* b, bignum* c, int fa, int fb, int* fc) {
	if (fa > 0 && fb > 0) {
		if (bgt(*a, *b)) {
			*c = bsub(*a, *b);
			*fc = 1;
		}
		else {
			*c = bsub(*b, *a);
			*fc = -1;
		}
	}
	else if (fa > 0 && fb < 0) {
		*c = badd(*a, *b);
		*fc = 1;
	}
	else if (fa < 0 && fb > 0) {
		*c = badd(*a, *b);
		*fc = -1;
	}
	else if (fa < 0 && fb < 0) {
		if (bgt(*b, *a)) {
			*c = bsub(*b, *a);
			*fc = 1;
		}
		else {
			*c = bsub(*a, *b);
			*fc = -1;
		}
	}
}	

void signedbadd(bignum* a, bignum* b, bignum* c, int fa, int fb, int* fc) {
	if (fa > 0 && fb < 0) {    
		if (bgt(*a, *b)) {
			*c = bsub(*a, *b);
			*fc = 1;
		}
		else {
			*c = bsub(*b, *a);
			*fc = -1;
		}
	}
	else if (fa > 0 && fb > 0) {
		*c = badd(*a, *b);
		*fc = 1;
	}
	else if (fa < 0 && fb > 0) { 
		*c = badd(*a, *b);
		*fc = -1;
	}
	else if (fa < 0 && fb > 0) { 
		if (bgt(*b, *a)) {
			*c = bsub(*b, *a);
			*fc = 1;
		}
		else {
			*c = bsub(*a, *b);
			*fc = -1;
		}
	}
}	

void signedmul(bignum* a, bignum* b, bignum* c, int fa, int fb, int* fc, bignum* n) {
	*c = mmul(*a, *b, *n);
	*fc = fa * fb;
	if (fc < 0) 
		*c = bsub(*n, *c);
}

void bset(bignum* a, bignum* b, int* fa, int fb) {
	*a = *b;
	*fa = fb;
}

void gcd(bignum* a, bignum* b, bignum* r, bignum* s, bignum* t, int* fs, int* ft, bignum* n) {
	bignum a0, b0, s0, t0, q, QT, QS, temp;
	int fs0, ft0, fq, ftemp, fQT, fQS;
	a0 = *a; 
	b0 = *b; 
	t0.h = 0, t0.l = 0, ft0 = 1; 
	t->h = 0, t->l = 1, *ft = 1;   
	s0.h = 0, s0.l = 1, fs0 = 1; 
	s->h = 0, s->l = 0, *fs = 1;  
	q = bquo(a0, b0), fq = 1;   
	*r = brem(a0, b0);          
	while (r->l > 0 || r->h > 0) { 
		signedmul(&q, t, &QT, fq, *ft, &fQT, n); signedbsub(&t0, &QT, &temp, ft0, fQT, &ftemp);
		bset(&t0, t, &ft0, *ft); 
		bset(t, &temp, ft, ftemp); 
		signedmul(&q, s, &QS, fq, *fs, &fQS, n); signedbsub(&s0, &QS, &temp, fs0, fQS, &ftemp);
		bset(&s0, s, &fs0, *fs); 
		bset(s, &temp, fs, ftemp); 
		a0 = b0;
		b0 = *r;
		q = bquo(a0, b0), fq = 1; 
		*r = brem(a0, b0);         
	}
	
	*r = b0;
}

void inverse(bignum* a, bignum* b, bignum* c) {
	bignum r, s, t;
	int fs, ft;
	gcd(&rsaN, b, &r, &s, &t, &fs, &ft, &rsaN);
	
	if (ft < 0) {
		r = brem(t, rsaN);
		t = bsub(rsaN, r);
	}

	*c = brem(t, rsaN);
	*c = mmul(*a, *c, rsaN);
}
//********************************************************




void Ersa (bignum m, bignum *ans)
{
	*ans = m;
	m = mmul(m,m,rsaN);
	m = mmul(m,m,rsaN);
	m = mmul(m,m,rsaN);
	m = mmul(m,m,rsaN);
	*ans = mmul(*ans,m,rsaN);
}

bignum enumerate(int ij, int hardness, bignum hint)
{
	long upper; MPI_Request request;
	bignum ans; MPI_Status status;
	int i;
	
	ans.h = hint.h;
	if (!myrank)
		MPI_Irecv(&ans.l, 1, MPI_UNSIGNED_LONG, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &request);
	else
		MPI_Irecv(&ans.l, 1, MPI_UNSIGNED_LONG, 0, 0, MPI_COMM_WORLD, &request);
	
	upper = hint.l+(1<<hardness);
	if (hint.l<upper-myrank&&hint.l+myrank<upper) for (hint.l+=myrank; ; hint.l+=numprocs)
	{
		bignum testans; int flag;
		
//		if (myrank) fprintf(stderr,"Proccessor %d: hint.l = %ld",myrank,hint.l);
		
		Ersa(hint, &testans);
//		if (myrank) fprintf(stderr," after Ersa that is"); bputnum10(testans); fprintf(stderr,"\n");
		if (beq(testans,btransnum(0,ij)))
			MPI_Send(&hint.l, 1, MPI_UNSIGNED_LONG, 0, 0, MPI_COMM_WORLD);
		MPI_Test(&request, &flag, &status);
		
		if (flag)
			if (!myrank)
			{
				for (i=1; i<numprocs; i++) MPI_Send(&ans.l, 1, MPI_UNSIGNED_LONG, i, 0, MPI_COMM_WORLD);
				return ans;
			}
			else return ans;
		if (hint.l+numprocs<upper&&hint.l<upper-numprocs)
		{}
		else break;
	}
	MPI_Wait(&request, &status);
	if (!myrank)
	{
		for (i=1; i<numprocs; i++) MPI_Send(&ans.l, 1, MPI_UNSIGNED_LONG, i, 0, MPI_COMM_WORLD);
		return ans;
	}
	else return ans;
}

void xx()
{
	int bfs[Nxy], t,w,k,i,j, hash[Nxy];
	const int d[8][2] = {{0,1},{1,0},{0,-1},{-1,0},{1,1},{1,-1},{-1,1},{-1,-1}};

	memset(hash,0,sizeof(hash));
	t = 0, w = 0;
	for (i=0; i<Nxy; i++) if (!hardness[i]) bfs[w++] = i;
	while (t<w)
	{
		for (k=0; k<4; k++)
		{
			i = bfs[t]/256+d[k][0];
			j = bfs[t]%256+d[k][1];
			if (i>=0&&i<256&&j>=0&&j<256)
				if (hardness[i*256+j]&&!hash[i*256+j])
				{
					hash[i*256+j] = 1;
					image[i*256+j] = image[bfs[t]];
					bfs[w++] = i*256+j;
				}
		}
		t++;
	}
}

int main(int argc, char **argv)
{
	FILE *fp1, *fp2; int i,j,x, total;
	
 	total = 0;
 	fp1 = fopen(argv[1], "r");	
	fp2 = fopen(argv[2], "r");
//	freopen(stdout,"xx.txt","wr");
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
	MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
	
	if (myrank == 0)
	{
		outstarttime();

		getimage(fp1, eimage);
		gethint(fp2, &rsaN, &rsaE, hardness, hint);
		MPI_Bcast(&rsaN.h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaN.l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaE.h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaE.l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		list_n = 0;
		for (i = 0; i < Nxy; i++)
		{
//			if (i%100==0) fprintf(stderr,"Proccessor %d start on %dth place\n",myrank,i);
			if (i%10000==0)
			{
				for (j = 0; j < Nxy; j++)
					image[j] = eimage[j] ^ (hint[j].l&1);
				xx();
				outresult(image);
			}
			if (hardness[i] < Maxhardness || i<=20 && hardness[i]<=29)
			{
				int now;
				now = 0;
				if (hardness[i]>=22)
					now = 1<<hardness[i]-22;
				if (hardness[i]<22||total+now<=MaxTime)
				{
					total += now;
//					if (hardness[i]>=Maxhardness)
//						fprintf(stderr,"Proccessor %d start on special number %d with hardness %d\n",myrank,i,hardness[i]);
					MPI_Bcast(&hardness[i], 1, MPI_INT, 0, MPI_COMM_WORLD);
					MPI_Bcast(&hint[i].l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
					MPI_Bcast(&hint[i].h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
					hint[i] = enumerate(i, hardness[i], hint[i]);
					hardness[i] = 0;
					if (i>1) list[list_n++] = i;
				}
				else
					MPI_Bcast(&TAG, 1, MPI_INT, 0, MPI_COMM_WORLD);
			}	
			else
				MPI_Bcast(&TAG, 1, MPI_INT, 0, MPI_COMM_WORLD);
		}
		hint[1].l = 1, hardness[0] = hardness[1] = 0;
		MPI_Bcast(hardness, Nxy, MPI_INT, 0, MPI_COMM_WORLD);
		
		memset(hash,0,sizeof(hash));
		last_list_n = 0;
		for (/*x=1*/;;/*x++*/)
		{
//			fprintf(stderr,"Proccessor %d: list = %d\n", myrank, list_n);
//			fprintf(stderr,"Proccessor %d: rsaN = %lx %lx\n", myrank, rsaN.h, rsaN.l);
			for (i = 0; i < Nxy; i++)
				image[i] = eimage[i] ^ (hint[i].l&1);
			xx();
			outresult(image);
			own_data_n = 0;
			for (i = myrank; i < list_n; i += numprocs)
			{
//				fprintf(stderr,"Proccessor %d: searching i list[i] own_data_n = %d %d %d\n",myrank,i,list[i],own_data_n);
				for (j = last_list_n>i ? last_list_n : i ; j < list_n; j++)
				{
					unsigned long ij,large,small;
//					if (x==2) fprintf(stderr, "Proccessor %d: searching i j list[i] list[j] = %d %d %d %d\n",myrank,i,j,list[i],list[j]);
//					if (x==2) fprintf(stderr, "%lx %lx  :  %lx %lx\n", hint[list[i]].h, hint[list[i]].l, hint[list[j]].h, hint[list[j]].l);
					ij = list[i]*list[j];
					if ( ij<Nxy && hardness[ij] && !hash[ij])
					{
						bignum temp;
						
						hash[ij] = 1;
						own_data[own_data_n++] = ij;
						temp = mmul(hint[list[i]],hint[list[j]],rsaN);
						own_data[own_data_n++] = temp.l;
						own_data[own_data_n++] = temp.h;
					}

					if (list[i]<list[j])
						large = list[j], small = list[i];
					else
						large = list[i], small = list[j];
					if (large%small==0)
					{
						ij = large/small;
						if ( ij>1 && hardness[ij] && !hash[ij] )
						{
							bignum temp;
							
//							fprintf(stderr,"Proccessor %d INVERSE: i j = %d %d\n",myrank, large, small);
							hash[ij] = 1;
							own_data[own_data_n++] = ij;
//							fprintf(stderr,"ij = %d\n", own_data[own_data_n-1]);
							inverse(&hint[large], &hint[small], &temp);
//							fprintf(stderr,"%lx %lx\n",temp.h,temp.l);
							own_data[own_data_n++] = temp.l;
							own_data[own_data_n++] = temp.h;
						}
					}
				}
			}
//			fprintf(stderr,"Proccessor %d: own_data_n = %d\n", myrank, own_data_n);
			MPI_Allgather(&own_data_n, 1, MPI_INT, arr_data, 1, MPI_INT, MPI_COMM_WORLD);
			displs[0] = 0;
			for (i=1; i<=numprocs; i++) displs[i] = displs[i-1] + arr_data[i-1];
			if (!displs[numprocs]) break;
			MPI_Allgatherv(own_data, own_data_n, MPI_UNSIGNED_LONG, new_data, arr_data, displs, MPI_UNSIGNED_LONG, MPI_COMM_WORLD);

			last_list_n = list_n;
			own_data_n = 0;		//***********Another meaning of "own_data_n"
			for (; own_data_n<displs[numprocs]; )
			{
				int ij;
				ij = new_data[own_data_n++];
				if (hardness[ij])
				{
					hardness[ij] = 0;
					hint[ij].l = new_data[own_data_n++];
					hint[ij].h = new_data[own_data_n++];
					list[list_n++] = ij;
				}
				else own_data_n+=2;
			}
		}
		
		for (i = 0; i < Nxy; i++)
			image[i] = eimage[i] ^ (hint[i].l&1);
		xx();
		outresult(image);
	}
	else
	{
		MPI_Bcast(&rsaN.h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaN.l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaE.h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		MPI_Bcast(&rsaE.l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
		list_n = 0;
		for (i = 0; i < Nxy; i++)
		{
			int len;
//			if (i%100==0) fprintf(stderr,"Proccessor %d start on %dth place\n",myrank,i);
			MPI_Bcast(&len, 1, MPI_INT, 0, MPI_COMM_WORLD);
			if (len != TAG) {
				bignum hint_now;
				MPI_Bcast(&hint_now.l, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
				MPI_Bcast(&hint_now.h, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
				hint[i] = enumerate(i, len, hint_now);
				if (i>1) list[list_n++] = i;
			}
		}
		hint[1].l = 1;
		MPI_Bcast(hardness, Nxy, MPI_INT, 0, MPI_COMM_WORLD);

/*		if (myrank==1)
			for (i=0; i<list_n; i++)
				fprintf(stderr,"list[%d] = %d\n",i,list[i]);*/
		memset(hash,0,sizeof(hash));
		last_list_n = 0;
		for (/*x=1*/;;/*x++*/)
		{
//			fprintf(stderr,"Proccessor %d: list = %d\n", myrank, list_n);
			own_data_n = 0;
			for (i = myrank; i < list_n; i += numprocs)
			{
//				fprintf(stderr,"Proccessor %d: searching i list[i] own_data_n = %d %d %d\n",myrank,i,list[i],own_data_n);
				for (j = last_list_n>i ? last_list_n : i ; j < list_n; j++)
				{
					unsigned long ij,large,small;
//					if (x==2) fprintf(stderr, "Proccessor %d: searching i j = %d %d\n",myrank,i,j);
//					fprintf(stderr,"Proccessor %d: searching i j = %d %d\n",myrank,i,j);
					ij = list[i]*list[j];
					if ( ij<Nxy && hardness[ij] && !hash[ij])
					{
						bignum temp;
						
						hash[ij] = 1;
						own_data[own_data_n++] = ij;
						temp = mmul(hint[list[i]],hint[list[j]],rsaN);
						own_data[own_data_n++] = temp.l;
						own_data[own_data_n++] = temp.h;
					}
					
					if (list[i]<list[j])
						large = list[j], small = list[i];
					else
						large = list[i], small = list[j];
					if (large%small==0)
					{
						ij = large/small;
						if ( ij>1 && hardness[ij] && !hash[ij] )
						{
							bignum temp;
							
//							fprintf(stderr,"Proccessor %d INVERSE: i j = %d %d\n",myrank, large, small);
							hash[ij] = 1;
							own_data[own_data_n++] = ij;
							inverse(&hint[large], &hint[small], &temp);
//							fprintf(stderr,"%lx %lx\n",temp.h,temp.l);
							own_data[own_data_n++] = temp.l;
							own_data[own_data_n++] = temp.h;
						}
					}
				}
			}
//			fprintf(stderr,"Proccessor %d: own_data_n = %d\n", myrank, own_data_n);
			MPI_Allgather(&own_data_n, 1, MPI_INT, arr_data, 1, MPI_INT, MPI_COMM_WORLD);
			displs[0] = 0;
			for (i=1; i<=numprocs; i++) displs[i] = displs[i-1] + arr_data[i-1];
			if (!displs[numprocs]) break;
			MPI_Allgatherv(own_data, own_data_n, MPI_UNSIGNED_LONG, new_data, arr_data, displs, MPI_UNSIGNED_LONG, MPI_COMM_WORLD);

			last_list_n = list_n;
			own_data_n = 0;		//***********Another meaning of "own_data_n"
			for (; own_data_n<displs[numprocs]; )
			{
				int ij;
				ij = new_data[own_data_n++];
				if (hardness[ij])
				{
					hardness[ij] = 0;
					hint[ij].l = new_data[own_data_n++];
					hint[ij].h = new_data[own_data_n++];
					list[list_n++] = ij;
				}
				else own_data_n+=2;
			}
		}

	}
	
	MPI_Finalize();
	return 0;
}