/*
 * @author Aaron Blumenfeld
 * The following program implements Elliptic Curve arithmetic
 * over prime fields F_p.  F_2 and F_3 are not supported.
 * Projective coordinates are used, not affine.  Therefore,
 * to understand the points, (0 : 1 : 0) represents the point
 * at infinity; otherwise, the z-coordinate is 1, and
 * (x : y : 1) represents (x, y) in affine coordinates.  Points
 * are represented as structs.
 *
 * jacobi quickly computes Jacobi (i.e., Legendre) symbols.
 * addhelper adds EC points without reducing to unique (x : y : z)
 * add reduces to unique (x : y : z)
 * mult calls addhelper and reduces at end, using repeated doubling
 *
 * All the other methods have plenty of comments in there, so scroll down
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define bool char
#define false 0
#define FALSE 0
#define true 1
#define TRUE 1

struct Point { /* representation of Points (projective coordinates) for ECs */
    long x;
    long y;
    long z;
};

long trueMod(long a, long p) { /* does mod, not remainder */
	if(a >= 0)
	    return a % p;
	else
		return (-a * (p-1)) % p;
}

long modExp(long a, long b, long p) {
	long rval = 1;
	while(b > 0) {
		if(!(b & 1)) /* if b is odd */
			rval = (rval * a) % p;
		b >>= 1;
		a = (a * a) % p;
	}
	return rval;
}

long inverse(long a, long p) { /* uses Fermat */
	return modExp(a, p-2, p);
}

bool equals(struct Point P, struct Point Q) {
	return P.x == Q.x && P.y == Q.y && P.z == Q.z;
}

int jacobi(long aa, long p) { /* somewhat faster than a^(p-1)/2 (mod p) */
	long a = trueMod(aa, p);
	if(a == 0)
		return 0;
	int rval = 1;
	int mod8;
	long temp;

	while(a != 0) {
		while(a % 2 == 0) { /* pull out factors of 2 and compute (2/n) */
			a >>= 1;
			mod8 = (int) trueMod(p, 8);
			if(mod8 == 3 || mod8 == 5) {
				if(rval == 1)
					rval = -1;
				else
					rval = 1;
			}
		}
	temp = a;
	a = p; /* swap a and p */
	p = temp;

		if(trueMod(a, 4) == 3 && trueMod(p, 4) == 3) { /* apply quadratic reciprocity */
			if(rval == 1)
				rval = -1;
			else
				rval = 1;
		}
		a = trueMod(a, p);
	}
	return rval;
}

/* this routine separate so mult doesn't invert the z-coordinate O(logk) times */
struct Point addhelper(struct Point P, struct Point Q, long a, long p) {
	struct Point R, negQ;
	negQ.x = Q.x; negQ.y = p - Q.y; negQ.z = Q.z;
	long u, v, w, t;
	struct Point inf;
	inf.x = 0; inf.y = 1; inf.z = 0;
	if(equals(P, inf)) /* if P = infinity ==> P + Q = Q */
		return Q;
	else if(equals(Q, inf)) /* if Q = infinity ==> P + Q = P */
	    return P;
	else if(!equals(P, Q) && !equals(P, negQ)) { /* if P != +-Q */
		u = trueMod(trueMod(Q.y*P.z, p) - trueMod(P.y*Q.z, p), p);
		v = trueMod(trueMod(Q.x*P.z, p) - trueMod(P.x*Q.z, p), p);
		w = trueMod(trueMod(trueMod(u*u, p)*trueMod(P.z*Q.z, p), p) - trueMod(trueMod(v*v, p)*v, p) - trueMod(trueMod(2*v*v, p), p)*trueMod(P.x*Q.z, p), p);
		R.x = trueMod(v*w, p);
		R.y = trueMod(u*(trueMod(v*v, p)*trueMod(P.x*Q.z, p) - w), p) - trueMod(trueMod(trueMod(v*v, p)*v, p)*trueMod(P.y*Q.z, p), p);
		R.z = trueMod(trueMod(trueMod(v*v, p)*v, p)*trueMod(P.z*Q.z, p), p);
		return R;
	}
	else if(equals(P, Q)) { /* if P = Q */
		t = trueMod(trueMod(trueMod(a*P.z, p)*P.z, p) + trueMod(3*P.x*P.x, p), p);
		u = trueMod(P.y*P.z, p);
		v = trueMod(trueMod(u*P.x, p)*P.y, p);
		w = trueMod(trueMod(t*t, p) - trueMod(8*v, p), p);
		R.x = trueMod(2*u*w, p);
		R.y = trueMod(trueMod(t*(4*v - w), p) - trueMod(trueMod(trueMod(8*P.y, p)*trueMod(P.y*u, p), p)*u, p), p);
		R.z = trueMod(trueMod(8*u, p)*trueMod(u*u, p), p);
		return R;
	}
	else
		return inf; /* if P = -Q ==> P + Q = infinity */
}

/* call addhelper and reduce the coordinate at the end, since the result is
an equivalence class, so it must be scaled back down to reduced coordinates */
struct Point add(struct Point P, struct Point Q, long a, long p) {
	struct Point R;
	R = addhelper(P, Q, a, p);
	if(R.z%p == 0) {
		R.y = 1; /* scale infinity down to unique (0, 1, 0) */
	}
	else {
		long inv = inverse(R.z, p);
		R.x *= inv;
		R.y *= inv; /* scale z coordinate back down to 1 */
		R.z *= inv;
	}
	R.x = trueMod(R.x, p);
	R.y = trueMod(R.y, p);
	R.z = trueMod(R.z, p); /* reduce coordinates mod p */
	return R;
}

/* repeated doubling implemented from Washington's pseudocode */
struct Point mult(long k, struct Point P, long a, long p) {
	long A = k;
	struct Point B, C;
	B.x = 0; B.y = 1; B.z = 0;
	C.x = P.x; C.y = P.y; C.z = P.z;
	while(A > 0) {
		if(!(A&1)) { /* if A is even */
			A >>= 1;
			C = addhelper(C, C, a, p);
			C.x = trueMod(C.x, p);
			C.y = trueMod(C.y, p); /* reduce mod p O(logk) times to prevent overflow */
			C.z = trueMod(C.z, p);
		}
		else { /* if A is odd */
			A -= 1;
			B = addhelper(B, C, a, p);
			B.x = trueMod(B.x, p);
			B.y = trueMod(B.y, p); /* reduce mod p O(logk) times to prevent overflow */
			B.z = trueMod(B.z, p);
		}
	}
	if(trueMod(B.z, p) == 0)
		B.y = 1; /* scale infinity down to unique (0 : 1 : 0) */
	else {
		long inv = inverse(B.z, p);
		B.x *= inv;
		B.y *= inv; /* scale z coordinate back down to 1 */
		B.z *= inv;
	}
	B.x = trueMod(B.x, p);
	B.y = trueMod(B.y, p);
	B.z = trueMod(B.z, p); /* reduce coordinates mod p */
	return B;
}

long ecOrder(long a, long b, long p) { /* O(plogp) algorithm, maybe fix
	     later using Shoof's algorithm for O(log^8 p) run-time */
	long x;
	long order = 1;
	for(x = 0; x < p; x++)
		order += 1 + jacobi(trueMod(trueMod(x*x, p)*x, p) + trueMod(a*x, p) + b, p);
	return order;
}

long termNodes(long a, long b, long p, long N) {
	long x;
	long count = 0;
	for(x = 0; x < p; x++)
		if(jacobi(trueMod(trueMod(x*x, p)*x, p) + trueMod(a*x, p) + b, p) == -1)
		    count++;
	count += (N - p); /*(k, ?) can't get mapped to if k >= p... */
	return count;
}

/* Brent's algorithm, which is supposedly faster than Floyd's algorithm */
void tailCycleLength(long a, long b, long p, long k, struct Point B, long res[]) {
	long power, lambda; power = lambda = 1;
	long tortoise, hare;
	tortoise = k; hare = mult(k, B, a, p).x;
	while(tortoise != hare) { /* search successive powers of 2 */
		if(power == lambda) {
			tortoise = hare;
			power *= 2;
			lambda = 0;
		}
		hare = mult(hare, B, a, p).x;
		lambda++;
	}

	/* find the position of the first repetition of length lambda */
	long mu = 0;
	tortoise = hare = k;
	long i = 0;
	for(; i < lambda; i++) {
		hare = mult(hare, B, a, p).x;
	}
	while(tortoise != hare) {
		tortoise = mult(tortoise, B, a, p).x;
		hare = mult(hare, B, a, p).x;
		mu++;
	}
	res[0] = mu; res[1] = lambda; /* this is void and modifies an array of length 2
									 so that it can "return" 2 values */
}

/* terminal/image nodes: compute with Jacobi symbols
 * fixedPoints, twoCycles, threeCycles, fiveCycles: loop through using an eZ if statement
 * components = cycles, avg cycle size: #cyclic nodes/#components
 * cyclicNodes/tailNodes: uses Brent's algorithm for cycle finding
 * avgCycle/avgTail: just sum over the array and divide by (order+1)
 * maxCycle, maxTail: just update when finding a new largest tail or cycle
 *
 * When outputting a line to a file, have a line look like the following:
 * termNodes imNodes 1-cycs 2-cycs 3-cycs comps avgCycSize cycNodes tailNodes avgCyc avgTail maxCyc maxTail
 */
void compute(long a, long b, long p, long order, struct Point B, long terminalNodes, long imageNodes, long fixedPoints, long twoCycles, long threeCycles, long fiveCycles,
long components, long cyclicNodes, long tailNodes, long cycleSum, long tailSum, long maxCycle, long maxTail,
double avgCycle, double avgTail, double avgCycleSize, FILE *out) {
	long cLength[order+1];
	long tLength[order+1];
	long k = 0; /* use as an index for iteration */

	/* compute # cyclic/tail nodes, maxCycle/maxTail length */
	tailNodes = 1; /* 0 is a tail node */
	cyclicNodes = 1; /* infinity is a cyclic node */
	maxCycle = 1;
	maxTail = 0;
	long temp[2]; temp[0] = -1; temp[1] = -1; /* use to store tail and cycle length of a node */
	cLength[0] = 1; tLength[0] = 1; /* use 0th index for 0 */
	cLength[order] = 1; tLength[0] = 0; /* use order'th index for infinity */
	for(k = 1; k < order; k++) {
		tailCycleLength(a, b, p, k, B, temp);
		tLength[k] = temp[0]; cLength[k] = temp[1];
		if(tLength[k] > maxTail)
			maxTail = tLength[k];
		if(cLength[k] > maxCycle)
			maxCycle = cLength[k];
		if(tLength[k] == 0)
			cyclicNodes++;
	}
	tailNodes = order+1 - cyclicNodes;
	/* end compute # cyclic/tail nodes, maxCycle/maxTail length */

	/* compute avg cycle / avg tail */
	cycleSum = 2;
	tailSum = 1;
	for(k = 1; k < order; k++) {
		cycleSum += cLength[k];
		tailSum += tLength[k];
	}
	avgCycle = (double)cycleSum; avgCycle /= (double)(order+1);
	avgTail = (double)tailSum; avgTail /= (double)(order+1);
	/* end compute avg cycle / avg tail */

	/* count terminal/image nodes */
	terminalNodes = termNodes(a, b, p, order);
	imageNodes = order + 1 - terminalNodes;
	/* end count terminal/image nodes */

	/* find various small cycles */
	fixedPoints = 1;
	twoCycles = 0;
	threeCycles = 0;
	fiveCycles = 0;
	long temp1 = 0;
	for(k = 1; k < order; k++) {
		temp1 = mult(k, B, a, p).x;
		if(temp1 == k)
			fixedPoints++;
		else if((temp1 = mult(temp1, B, a, p).x) == k)
			twoCycles++;
		else if((temp1 = mult(temp1, B, a, p).x) == k)
			threeCycles++;
		else if((temp1 = mult(mult(temp1, B, a, p).x, B, a, p).x) == k)
			fiveCycles++;
	}
	twoCycles >>= 1; /* since 2 representations of 2-cycles */
	threeCycles /= 3; /* since 3 representations of 3-cycles */
	fiveCycles /= 5; /* since 5 representations of 5-cycles */
	/* end find various small cycles */

	/* compute # components, average cycle */
	components = 0;
	long cSize[maxCycle+1];
	cSize[0] = 0; cSize[1] = 1;
	for(k = 2; k <= maxCycle; k++) /* initialize array */
		cSize[k] = 0;
	for(k = 1; k < order; k++)
		if(tLength[k] == 0) /* if k is a cyclic node */
			cSize[cLength[k]]++; /* increment # cyclic nodes by its cycle length */
	for(k = 1; k <= maxCycle; k++) {
		cSize[k] /= k; /* since there are k representations of any k-cycle */
		components += cSize[k]; /* increment components by the number of k-cycles for each k */
	}
	avgCycleSize = (double)cyclicNodes; avgCycleSize /= (double)components; /* avg. cycle length = #cyclic nodes/#components */
	/* end compute # components, average cycle */

	fprintf(out, "%d %d %d %d %d %d %d %f %d %d %f %f %d %d\n", terminalNodes, imageNodes,
	fixedPoints, twoCycles, threeCycles, fiveCycles, components, avgCycleSize, cyclicNodes,
	tailNodes, avgCycle, avgTail, maxCycle, maxTail);

	printf("y^2 = x^3 + %dx + %d mod %d with %d points, B = (%d, %d)\n", a, b, p, order, B.x, B.y);
	/*printf("%d terminal nodes, %d image nodes\n", terminalNodes, imageNodes);
	printf("%d fixed points, %d 2-cycles, %d 3-cycles, %d 5-cycles\n", fixedPoints, twoCycles, threeCycles, fiveCycles);
	printf("%d components, %f avg cycle size\n", components, avgCycleSize);
	printf("%d cyclic nodes, %d tail nodes\n", cyclicNodes, tailNodes);
	printf("%f avg cycle length, %f avg tail length\n", avgCycle, avgTail);
	printf("%d max cycle length, %d max tail length\n\n", maxCycle, maxTail);*/
}

void printECpts(long a, long b, long p) { /* this function was only used for testing purposes */
	long x, y;
	for(x = 0; x < p; x++)
		for(y = 0; y < p; y++)
			if(trueMod(y*y - x*x*x - a*x - b, p) == 0)
				printf("(%d, %d, 1) is on E\n", x, y);
	printf("(0, 1, 0) is on E");
}


void printBmults(struct Point B, long a, long b, long p) { /* this function was only used for testing purposes */
	struct Point Q; Q.x = 0; Q.y = 1; Q.z = 0;
	long i = 1; long order = ecOrder(a, b, p);
	for(; i <= order; i++) {
		Q = add(B, Q, a, p);
		printf("%dB = (%d, %d)\n", i, Q.x, Q.y);
	}
}

int main(void) {
	long a, b, p;
	struct Point B; B.z = 1;
	long order = ecOrder(a, b, p);
	long N = 997; /* change file names by hand based on N */
	static const char filename[] = "997curves.txt";
	static const char fileout[] = "997out.txt";
	FILE *file = fopen(filename, "r"); /* read only */
	FILE *out;
	out = fopen(fileout, "a+"); /* append file (add text to a file) or create a file if it does not exist */
	if(file != NULL) {
		char line[32]; /* 32 chars long should be fine for numbers < 1 million */
		while(fgets(line, sizeof line, file) != NULL) {
			a = atol(line); fgets(line, sizeof line, file);
			b = atol(line); fgets(line, sizeof line, file);
			p = atol(line); fgets(line, sizeof line, file);
			B.x = atol(line); fgets(line, sizeof line, file);
			B.y = atol(line);

			long terminalNodes, imageNodes, fixedPoints, twoCycles, threeCycles, fiveCycles, components, cyclicNodes,
			tailNodes, cycleSum, tailSum, maxCycle, maxTail;
			double avgCycle, avgTail, avgCycleSize;
			compute(a, b, p, N, B, terminalNodes, imageNodes, fixedPoints, twoCycles, threeCycles,
			fiveCycles, components, cyclicNodes, tailNodes, cycleSum, tailSum, maxCycle, maxTail,
			avgCycle, avgTail, avgCycleSize, out);
		}
		fclose(file);
	}
	else
		perror(filename);
}