/**
 * @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.
 *
 * Elliptic curves are of the form y^2 = x^3 + ax + b.  They have
 * no repeated roots on the right side of the equation, which is
 * equivalent to 4a^3 + 27b^2 != 0 (mod p).  The isEC method
 * tests if a given equation is an elliptic curve or not.
 *
 * One other concern is that integers modulo p are represented as
 * longs, not as BigIntegers, so for primes hundreds of digits long,
 * this program will not work, but for the purpose of moderately large
 * primes, 64-bit numbers will suffice.
 */

import java.io.BufferedReader;
import java.io.FileOutputStream;
import java.io.FileReader;
import java.io.PrintStream;
import java.util.LinkedList;
import java.util.ArrayList;

public class EC {

	public static long trueMod(long a, long b) { /*fixes Java % feature when dealing
	        with negative numbers*/
		if (a >= 0)
			return a % b;
		else
			return (-a * (b-1)) % b;
	}

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

	public static long inverse(long a, long p) { /* lazy way to compute a^(-1) (mod p) */
		return modExp(a, p-2, p);
	}

	public static Point addhelper(Point P, Point Q, long a, long p) { /* uses formulas
	        for projective coordinates from Washington */
		/* this method separate so mult doesn't invert the z-coordinate O(logk) times */
		long x3, y3, z3;
		if(!P.equals(Q) && !P.equals(new Point(Q.x, trueMod(-Q.y, p), Q.z))) { /* if P != +-Q */
			long u = trueMod(Q.y*P.z, p) - trueMod(P.y*Q.z, p);
			long v = trueMod(Q.x*P.z, p) - trueMod(P.x*Q.z, p);
			long w = trueMod(trueMod(u*u, p)*trueMod(P.z*Q.z, p), p) - trueMod(v*v*v, p) - trueMod(trueMod(2*v*v, p), p)*trueMod(P.x*Q.z, p);
			x3 = trueMod(v*w, p);
			y3 = trueMod(u*(trueMod(v*v, p)*trueMod(P.x*Q.z, p) - w), p) - trueMod(trueMod(v*v*v, p)*trueMod(P.y*Q.z, p), p);
			z3 = trueMod(trueMod(v*v*v, p)*trueMod(P.z*Q.z, p), p);
		}
		else if(P.equals(Q)) { /* if P = Q */
			long t = trueMod(a*P.z*P.z, p) + trueMod(3*P.x*P.x, p);
			long u = trueMod(P.y*P.z, p);
			long v = trueMod(u*P.x*P.y, p);
			long w = trueMod(t*t, p) - trueMod(8*v, p);
			x3 = trueMod(2*u*w, p);
			y3 = trueMod(t*(4*v - w), p) - trueMod(8*P.y*P.y*u*u, p);
			z3 = trueMod(8*u*u*u, p);
		}
		else { /* if P = -Q ==> P + Q = infinity */
			x3 = 0;
			y3 = 1;
			z3 = 0;
		}
		if(P.equals(new Point(0, 1, 0))) { /* if P = infinity ==> P + Q = Q */
			x3 = Q.x; y3 = Q.y; z3 = Q.z;
		}
		if(Q.equals(new Point(0, 1, 0))) { /* if Q = infinity ==> P + Q = P */
			x3 = P.x; y3 = P.y; z3 = P.z;
		}
		return new Point(x3, y3, z3);
	}

	public static Point add(Point P, Point Q, long a, long p) { /* 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 */
		Point R = addhelper(P, Q, a, p);
		if(trueMod(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;
	}

	public static Point mult(long k, Point P, long a, long p) { /* compute kP
	    using algorithm from Washington */
		long A = k; /* should fix later to reduce mod N, the order of E */
		Point B = new Point(0, 1, 0);
		Point C = new Point(P);
		while(A > 0) {
			if(A % 2 == 0) { /* if A is even */
				A /= 2;
				C = new Point(addhelper(C, C, a, p));
				C.x = trueMod(C.x, p);
				C.y = trueMod(C.y, p); /* have to reduce mod p O(logk) times to prevent overflow */
				C.z = trueMod(C.z, p);
			}
			else { /* if A is odd */
				A -= 1;
				B = new Point(addhelper(B, C, a, p));
				B.x = trueMod(B.x, p);
				B.y = trueMod(B.y, p); /* have to 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;
	}

	public static boolean isEC(long a, long b, long p) {
		if(trueMod(4*a*a*a + 27*b*b, p) != 0) /* make sure no multiple roots */
			return true;
		return false;
	}

	public static boolean soluble(long a, long b, long p) {
		for(long x = 0; x < p; x++)
			if(trueMod(x*x*x + a*x + b, p) == 0)
				return true;
		return false;
	}

	public static void printECinfo(Point B, long a, long b, long p) {
		System.out.println("E : y^2 = x^3 + " + a + "x + " + b + " with base point " + B);
	}

	public static int legendre(long a, long p) { /* compute the Legendre symbol (a/p) */
		if(trueMod(a, p) == 0)
			return 0;
		long rval = modExp(a, (p-1)/2, p);
		if(rval == 1)
			return 1;
		else
			return -1;
	}

	public static long jacobi(long a, long p){
		if((a % p) == 0)
			return 0;
		long rval = 1;
		long mod8;
		long temp;

		a = (a % p);
		while(a != 0) {
			while(a % 2 == 0) { /* pull out factors of 2 and compute (2/n) */
				a >>= 1;
				mod8 = (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((a % 4) == 3 && (p % 4) == 3) { /* apply quadratic reciprocity */
				if(rval == 1)
					rval = -1;
				else rval = 1;
			}
			a = a % p;
		}
		return rval;
	}

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

	public static void printECpts(EllipticCurve E) {
		for(long x = 0; x < E.p; x++)
			for(long y = 0; y < E.p; y++)
				if(trueMod(y*y, E.p) == trueMod(x*x*x + E.a*x + E.b, E.p))
					System.out.println(new Point(x, y, 1) + " is on E");
		System.out.println(new Point(0, 1, 0) + " is on E");
	}

	public static void printBmults(Point B, EllipticCurve E) {
		Point Q = new Point(0, 1, 0);
		for(int i = 1; i <= ecOrder(E); i++) {
			Q = add(B, Q, E.a, E.p);
			System.out.println(i + "B = " + Q);
		}
	}

	public static void ECsoforderN(long p, long N) {
		try {
			PrintStream fout = new PrintStream(new FileOutputStream("data.txt"));
			for(int a = 0; a < p; a++)
				for(int b = 0; b < p; b++)
					if(isEC(a, b, p) && ecOrder(new EllipticCurve(a, b, p)) == N) {
						System.out.println(new EllipticCurve(a, b, p) + " has order " + N);
						//System.out.println(new EllipticCurve(a, b, p).forParsing());
						fout.print(a + " " + b + " " + p + "\n");
					}
		} catch(Exception e) {
			System.err.println(e);
		}
	}

	public static ArrayList<String> readData() {
		ArrayList<String> dat = new ArrayList<String>();
		try { //reads file
			BufferedReader in = new BufferedReader(new FileReader("data.txt"));
			for (int i = 0; in.ready(); i++)
				dat.add(in.readLine());
		} catch (Exception e) {
			System.err.println(e);
		}
		return dat;
	}

	public static long shanks(long a, long p) { /* compute sqrt(a) (mod p) */
		if((p % 4) == 3) /* easy for p = 3 (mod 4) */
			return modExp(a, (p+1)/4, p);
		if((p % 8) == 5 && modExp(a, (p-1)/4, p) == 1) /* sometimes easy for p = 5 (mod 8) */
			return modExp(a, (p+3)/8, p);
		long e = 0; long q = p-1;
		while(q%2==0 && q > 0) {
			e++;
			q >>= 1;
		}
		long n = 2;
		while(jacobi(n, p) != -1) {
			n = (long)(Math.random()*p);
		}
		long z = modExp(n, q, p);
		long y, r, x, b, t;
		y = z; r = e; x = modExp(a, (q-1)/2, p); b = (a*x*x)%p; x = (a*x)%p;
		while(b%p != 1) {
			int m = 1;
			m = smallestM(b, m, p);
			t = modExp(y, 1<<(r-m-1), p);
			y = (t*t)%p;
			r = m;
			x = (x*t)%p;
			b = (b*y)%p;
		}
		return x;
	}

	public static void main(String[] args) {
		long a, b, p, N;
		a = 1; b = 1; p = 17;
		N = 18;
		//ECsoforderN(17, 18); /* list all elliptic curves over F_p with order N */
		EllipticCurve E = new EllipticCurve(a, b, p);
		long order = ecOrder(E);
		//printECpts(E); /* list the points on E */
		Point B = new Point(0, 1, 1);
		printBmults(B, E); /* list the multiples of the generator B */
		/*ArrayList<String> data = readData();
		ArrayList<EllipticCurve> list = new ArrayList<EllipticCurve>();
		for(String s : data)
			list.add(parseEC(s));*/
	}
}