/*
 *  ec_stats.cpp
 *  
 *
 *  Created by Christopher Evans on 6/19/11.
 *  Copyright 2011 University of Arkansas. All rights reserved.
 *
 */


/*
 * @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>
#include <iostream>
#include <string>
#include <fstream>
#include <iomanip>
#include <sstream>

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

//int PRIME=59;

//std::string INPUT="59_curve_list.txt";
//std::string OUTPUT="stats_59";

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;
	}
*/
	while(b > 0) {
		rval = (rval * a) % p;
		b -= 1;
	}
	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) {
	long delta=0;
	Point R;
	if(equals(P,Q)){
//		delta=inverse((2*P.y),p);
		delta= (3*(P.x)*(P.x) + a)*inverse((2*P.y),p);
		R.x= trueMod(((delta*delta)-2*P.x),p);
		R.y= trueMod((delta*(P.x-R.x)-P.y),p);
		R.z=1;
	}else if(P.x==0 && P.y==1 && P.z==0){
		return Q;
	}else if(Q.x==0 && Q.y==1 && Q.z==0){
		return P;
	}else if(P.x==Q.x){
		Point inf;
		inf.x=0;
		inf.y=1;
		inf.z=0;
		return inf;
	}else{
		delta= (Q.y-P.y)*(inverse(trueMod((Q.x-P.x),p),p));
		R.x= trueMod(((delta*delta)-P.x-Q.x),p);
		R.y= trueMod((-P.y+(delta*(P.x-R.x))),p);
		R.z=1;
	}
	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:
 * 1-cycs 2-cycs 3-cycs 5-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 */
	
	//Uncesessary calculations
	
	/* 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, "%ld %ld %ld %ld %ld %f %ld %ld %f %f %ld %ld : ", /*terminalNodes, imageNodes,*/
			fixedPoints, twoCycles, threeCycles, fiveCycles, components, avgCycleSize, cyclicNodes,
			tailNodes, avgCycle, avgTail, maxCycle, maxTail);
	fprintf(out, "%ld, %ld, %ld, %ld, %ld \n", a, b, p, B.x, B.y);
//	printf("y^2 = x^3 + %ldx + %ld mod %ld with %ld points, B = (%ld, %ld)\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("(%ld, %ld, 1) is on E\n", x, y);
	printf("(0, 1, 0) is on E \n");
}


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);
	printf("B= %ld, %ld \n", B.x, B.y);
	for(; i <= order; i++) {
		Q = add(B, Q, a, p);
		printf("%ldB = (%ld, %ld)\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=0; 

//	std::cout << "Enter value for N" << std::endl;
	std::cin >> N;
	
	std::string s;
	std::stringstream out1;
	out1 << N <<"_curve_list.txt";
	s = out1.str();

	std::string filename=s;

	std::stringstream out2;
	out2<< "stats_" << N;
	s= out2.str();
	std::string fileout=s;
	
	

	
//	printBmults(B, 1, 6, 11);
/*
	std::cin >> N;
	std::cout << "Enter input file name: "<< std::endl;
	std::cin >> filename;
	std::cout << "Enter output filename: " << std::endl;
	std::cin >> fileout;
*/
//	filename= "/"+filename;
	
//	static const char filename[] = "997curves.txt";
//	static const char fileout[] = "997out.txt";
//	FILE *file = fopen(filename.c_str(), "r"); /* read only */


	
	std::ifstream file;
	file.open(filename.c_str(), std::ifstream::in);
	FILE *out;
	out = fopen(fileout.c_str(), "w"); // 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)
		while(!file.eof()){
			//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);
			
			file >> a;
			file.ignore(1);
			file >> b;
			file.ignore(1);
			file >> p;
			file.ignore(1);
			file >> B.x;
			file.ignore(1);
			file >> B.y;
			file.ignore(1);

//			printECpts(a,b,p);

//			printBmults(B, a, b, p);

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

	
}