#include "supercon.h"

#define SIZE 65536 /* (256*256) */

typedef struct
{
	int p1;
	int p2;
}ypair;

typedef struct
{
	int sign;
	bignum num;
}_bignum;

typedef struct
{
	int x;
	int image;
	bignum mimage;
}recv_msg;

int image[SIZE], eimage[SIZE];
bignum rsaN, rsaE;
int flgE[128], lenE;
bignum hint[SIZE];
int hardness[SIZE];
int s[SIZE][16];
bignum mimage[SIZE];

int isprime[SIZE];

int hokan[SIZE], outimage[SIZE];

MPI_Datatype mpi_bignum;
MPI_Datatype mpi_recv_msg;

ypair yp_data[800000];
ypair *yp[65536];
int yp_len[65536];

void hokan4(int x)
{
	int color;
	color = image[x] ? 1 : -1;
	if(x-Nx >= 0)
		hokan[x-Nx] += color;
	if(x+Nx < SIZE)
		hokan[x+Nx] += color;
	if(x-1 >= 0)
		hokan[x-1] += color;
	if(x+1 < SIZE)
		hokan[x+1] += color;
}

_bignum _bmul(_bignum x, _bignum y)
{
	_bignum a;

	a.sign = x.sign ^ y.sign;
	a.num = bmul(x.num, y.num);
	return a;
}

_bignum _bquo(_bignum x, _bignum y)
{
	_bignum a;

	a.sign = x.sign ^ y.sign;
	a.num = bquo(x.num, y.num);
	return a;
}

_bignum _bsub(_bignum x, _bignum y)
{
	_bignum a;

	if(x.sign ^ y.sign)
	{
		if(x.sign)
		{
			a.sign = 1;
			a.num = badd(x.num, y.num);
		}
		else
		{
			a.sign = 0;
			a.num = badd(x.num, y.num);
		}
	}
	else
	{
		if(x.sign)
		{
			if(bge(x.num, y.num))
			{
				a.sign = 1;
				a.num = bsub(x.num, y.num);
			}
			else
			{
				a.sign = 0;
				a.num = bsub(y.num, x.num);
			}
		}
		else
		{
			if(bge(y.num, x.num))
			{
				a.sign = 1;
				a.num = bsub(y.num, x.num);
			}
			else
			{
				a.sign = 0;
				a.num = bsub(x.num, y.num);
			}
		}
	}

	return a;
}

_bignum _badd(_bignum x, _bignum y)
{
	_bignum a;

	if(x.sign ^ y.sign)
	{
		if(x.sign)
		{
			if(bge(x.num, y.num))
			{
				a.sign = 1;
				a.num = bsub(x.num, y.num);
			}
			else
			{
				a.sign = 0;
				a.num = bsub(y.num, x.num);
			}
		}
		else
		{
			if(bge(y.num, x.num))
			{
				a.sign = 1;
				a.num = bsub(y.num, x.num);
			}
			else
			{
				a.sign = 0;
				a.num = bsub(x.num, y.num);
			}
		}
	}
	else
	{
		if(x.sign)
		{
			a.sign = 0;
			a.num = badd(x.num, y.num);
		}
		else
		{
			a.sign = 1;
			a.num = badd(x.num, y.num);
		}
	}

	return a;
}

bool _bgt(_bignum x, _bignum y)
{
	if(x.sign == 0)
	{
		if(y.sign == 0)
		{
			return bgt(x.num, y.num);
		}
		else
		{
			return true;
		}
	}
	else
	{
		if(y.sign == 0)
		{
			return false;
		}
		else
		{
			return bgt(y.num, x.num);
		}
	}
}

bignum ltob(unsigned long x)
{
	bignum a;
	a.h = 0;
	a.l = x;
	return a;
}


_bignum lto_b(int sign, unsigned long x)
{
	_bignum a;
	a.sign = sign;
	a.num.h = 0;
	a.num.l = x;
	return a;
}

bignum _btob(_bignum x)
{
	bignum a;
	a = x.num;
	return a;
}

_bignum bto_b(bignum x)
{
	_bignum a;
	a.sign = 0;
	a.num = x;
	return a;
}

int yakusu()
{
	int i, p1, p2, len;
	ypair *w;

	w = yp_data;
	for(i=1; i<65536; i++)
	{
		yp[i] = w;

		p1 = 2;
		p2 = i >> 1; /* p2 = i / 2; */
		len = 0;

		while(p1 < p2)
		{
			if(i - p1*p2 == 0 && p2 != 1)
			{
				w->p1 = p1;
				w->p2 = p2;
				w++;
				w->p2 = p1;
				w->p1 = p2;
				w++;
				len+=2;
			}
			p1++;
			p2 = i / p1;
		}

		yp_len[i] = len;
	}

	return w - yp_data;
}

int cmp(const void *a, const void *b)
{
	if( *(int *)a > *(int *)b )
		return 1;
	else if( *(int *)a == *(int *)b )
		return 0;
	else
		return -1;
}

void soinsuu(int s[], int n)
{
	int len;
	int m, i, add;

	/* init */
	len = 0;
	m = n;
	i = 2;

	if(n == 0)
	{
		s[0] = 0;
		return;
	}

	/* i = 2, 3 */
	while((m & 1) == 0) /* while(m % 2 == 0) */
	{
		m >>= 1; /* m /= 2; */
		s[len++] = 2;
	}
	while(m % 3 == 0)
	{
		m /= 3;
		s[len++] = 3;
	}

	/* i = 5, 7, 11, 13, 17, 19...... */
	i = 5;
 	add = 2;
	while(m > 1)
	{
		if(m % i == 0)
		{
			m /= i;
			s[len++] = i;
		}
		else
		{
			i += add;
			if(add == 2)
			{
				add = 4;
			}
			else if(add == 4)
			{
				add = 2;
			}
		}
		if(i <= m/5)
		{
			s[len++] = m;
			break;
		}
	}

	s[len] = 0;
}

bignum exmod(bignum b)
{
	int i;
	bignum a, c;

	a = brem(b, rsaN);
	c.h = 0;
	c.l = 1;

	for(i=0; i<lenE; i++)
	{
		if(flgE[i])
			c = mmul(c, a, rsaN);
		a = mmul(a, a, rsaN);
	}
	return c;
}

void solve(int x)
{
	long range;
	bignum a, b, c;

	a.h = 0;
	range = 1<<hardness[x];
	for(a.l=0; a.l<range; a.l++)
	{
		c = b = madd(hint[x], a, rsaN);
		b = mmul(b, b, rsaN);
		b = mmul(b, b, rsaN);
		b = mmul(b, b, rsaN);
		b = mmul(b, b, rsaN);
		c = mmul(c, b, rsaN);
		if(c.h == 0 && c.l == x)
			break;
	}
	image[x] = (a.l & 1) ^ eimage[x];
	mimage[x] = madd(hint[x], a, rsaN);
}

void EuclidEx(_bignum a, _bignum b, _bignum *x, _bignum *y)
{
	_bignum a0 = a, a1 = b, a2, x0 = lto_b(0,1), x1 = lto_b(0,0), x2, y0 = lto_b(0,0), y1 = lto_b(0,1), y2, q;

	while(a1.num.h != 0 || a1.num.l != 0)
	{
		q = _bquo(a0, a1);
		a2 = _bsub(a0, _bmul(q, a1));
		x2 = _bsub(x0, _bmul(q, x1));
		y2 = _bsub(y0, _bmul(q, y1));

		a0 = a1; a1 = a2;
		x0 = x1; x1 = x2;
		y0 = y1; y1 = y2;
	}
	if(a0.sign == 0)
	{
		*x = x0;
		*y = y0;
	}
	else
	{
		*x = x0;
		*y = y0;
		x->sign ^= 1;
		y->sign ^= 1;
	}
}

void solve_div(int ab, int a)
{
	int b;
	_bignum c1, alp, beta, N;
	bignum cb;

	c1.sign = 0;
	c1.num = mimage[a];
	N.sign = 0;
	N.num = rsaN;
	b = ab / a;

	EuclidEx(c1, N, &alp, &beta);

	if(alp.sign == 1)
		alp.num = bsub(badd(rsaN, bmul(rsaN, bquo(alp.num, rsaN))), alp.num);

	cb = mmul(mimage[ab], alp.num, rsaN);
	mimage[b] = cb;
	image[b] = (cb.l & 1) ^ eimage[b];
}

void solve_prod(int a, int b)
{
	int c;
	c = a*b;
	mimage[c] = mmul(mimage[a], mimage[b], rsaN);
	image[c] = (mimage[c].l & 1) ^ eimage[c];
}


void init_rsaE()
{
	int i, j;

	j = rsaE.l;
	for(i=0; i<64; i++, j>>=1)
	{
		if(flgE[i] = j&1)
			lenE = i+1;
	}
	j = rsaE.h;
	for(i=64; i<128; i++, j>>=1)
	{
		if(flgE[i] = j&1)
			lenE = i+1;
	}
}

void is_prime()
{
	int i, j;

	for(i=0; i<SIZE; i++)
	{
		isprime[i] = 1;
	}
	isprime[0] = 0;
	isprime[1] = 0;

	for(i=0; i<SIZE; i++)
	{
		if(isprime[i] == 1)
		{
			for(j=i*2; j<SIZE; j+=i)
			{
				isprime[j] = 0;
			}
		}
	}
}

void make_mpi_type()
{
	/* declaration for a new MPI type */
	int block[3];
	MPI_Datatype type[3];
	MPI_Aint disp[3];
	/* --- end of the declaration */

	/* making the new MPI type */
	type[0] = MPI_UNSIGNED_LONG;
	type[1] = MPI_UNSIGNED_LONG;
	disp[0] = 0;
	disp[1] = sizeof(unsigned long);
	block[0] = 1;
	block[1] = 1;
	MPI_Type_struct(2, block, disp, type, &mpi_bignum);
	MPI_Type_commit(&mpi_bignum);

	/* making the new MPI type */
	type[0] = MPI_INT;
	type[1] = MPI_INT;
	type[2] = mpi_bignum;
	disp[0] = 0;
	disp[1] = sizeof(int);
	disp[2] = sizeof(int)*2;
	block[0] = 1;
	block[1] = 1;
	block[2] = 1;
	MPI_Type_struct(3, block, disp, type, &mpi_recv_msg);
	MPI_Type_commit(&mpi_recv_msg);
}

void client_process()
{
	int x;
	MPI_Status status;
	recv_msg message;

	while(1)
	{
		MPI_Recv(&x, 1, MPI_INT, 0, 99, MPI_COMM_WORLD, &status);
		solve(x);
		message.x = x;
		message.image = image[x];
		message.mimage = mimage[x];
		MPI_Send(&message, 1, mpi_recv_msg, 0, 99, MPI_COMM_WORLD);
	}
}
void prod_search(int now)
{
	int i, end, k;

	k = SIZE / now;
	end = now * k;
	for(i=end; i>now; i-=now, k--)
	{
		if(image[k] != -1 && image[i] == -1)
		{
			solve_prod(now, k);
			hokan4(i);
			prod_search(i);
		}
	}
}

void div_search(int now)
{
	int i, end;
	ypair *pyp;

	pyp = yp[now];
	end = yp_len[now];
	for(i=0; i<end; i++)
	{
		if(image[pyp[i].p1] != -1 && image[pyp[i].p2] == -1)
		{
			solve_div(now, pyp[i].p1);
			hokan4(pyp[i].p2);
			div_search(pyp[i].p2);
		}
	}
}

void output()
{
	int i;
	for(i=0; i<SIZE; i++)
	{
		if(image[i] == -1)
			outimage[i] = hokan[i]>=0 ? 1 : 0;
		else
			outimage[i] = image[i];
	}
	outresult(outimage);
}

int main(int argc, char **argv)
{
	FILE *fp;
	int pri[SIZE];
	int send, i, now;
	recv_msg message;
	MPI_Status status;

	int size, myrank;

	double oldtime;

	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &size);
	MPI_Comm_rank(MPI_COMM_WORLD, &myrank);

	outstarttime();

	if(argc < 2)
	{
		MPI_Finalize();
		exit(1);
	}
	
	make_mpi_type();

	if(myrank == 0)
	{
		fp=fopen(argv[1], "r");
		getimage(fp, eimage);
		fclose(fp);

		fp=fopen(argv[2], "r");
		gethint(fp, &rsaN, &rsaE, hardness, hint);
		fclose(fp);

		is_prime();

		init_rsaE();

		for(i=0; i<SIZE; i++)
		{
			pri[i] = (hardness[i]<<16) + i;
			soinsuu(s[i], i);
			image[i] = -1;
			hokan[i] = 0;
		}
		qsort(pri, SIZE, sizeof(int), cmp);

		for(i=0; i<SIZE; i++)
		{
			pri[i] &= SIZE-1;
		}

		yakusu();
	}

	MPI_Bcast(hardness, SIZE, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(hint, SIZE, mpi_bignum, 0, MPI_COMM_WORLD);
	MPI_Bcast(eimage, SIZE, mpi_bignum, 0, MPI_COMM_WORLD);
	MPI_Bcast(&rsaN, 1, mpi_bignum, 0, MPI_COMM_WORLD);
	MPI_Bcast(&rsaE, 1, mpi_bignum, 0, MPI_COMM_WORLD);
	MPI_Bcast(flgE, 128, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&lenE, 1, MPI_INT, 0, MPI_COMM_WORLD);

	if(myrank == 0)
	{
		send = 0;
		oldtime = MPI_Wtime();

		image[0] = eimage[0];
		image[0] = 1 - eimage[0];

		for(send = 0; send < size-1; send++)
		{
			MPI_Send(&pri[send], 1, MPI_INT, send+1, 99, MPI_COMM_WORLD);
		}

		while(1)
		{
			MPI_Recv(&message, 1, mpi_recv_msg, MPI_ANY_SOURCE, 99, MPI_COMM_WORLD, &status);
			image[message.x] = message.image;
			mimage[message.x] = message.mimage;
			now = message.x;
			hokan4(now);

			if((send & 8) == 0) /* send % 8 == 0 */
			{
				for(i=0; i<SIZE; i++)
				{
					if(i < SIZE/2 && i >= 2 && image[i] != -1)
					{
						prod_search(i);
					}
					if(!isprime[i] && image[i] != -1)
					{
						div_search(i);
					}
				}
			}

			while(image[pri[send]] != -1) send++;
			MPI_Send(&pri[send], 1, MPI_INT, status.MPI_SOURCE, 99, MPI_COMM_WORLD);
			send++;

			if(MPI_Wtime() > oldtime + 5)
			{
					/*outresult(image);*/
					output();
					oldtime += 5;
			}
		}
	}
	else
	{
		client_process();
	}

	MPI_Finalize();
	exit(0);
}