#!/usr/bin/env python -i

from math import sqrt

# This file might seem a bit complicated, and it far above the level
# one would expect from most students. So, feel free to add your own 
# code at the marked places, or just write your own file. How to do it:
#
#	1. Do all the symbolic computations on paper.
#	2. Tell the computer to do the same computations.
#
# This code has drawn much inspiration from 
# http://jeremykun.com/2014/02/24/elliptic-curves-as-python-objects/
# which is available under a CC-BY-NC license, see
# http://creativecommons.org/licenses/by-nc/3.0/deed.en_US
#
# Last modified: 2014-09-17 by Martin Strand

# We can do all kind of computations mod something using %. However, we
# would like to avoid dragging that information about everywhere, so we
# make the modulus an inherent property of the number instead.
#
# This class could be extended to handle finite field extensions, but
# that is outside the scope of this exercise.
#
# What happens here is that the function FiniteRing returns a
# description of a class RingElement, which we in turn can make 
# elements of.

class ToolBox(object):

	# Extended Euclidean algorithm, which returns the gcd and the coefficients.
	def egcd(self,a, b):
	    if a == 0:
	        return (b, 0, 1)
	    else:
	        g, y, x = self.egcd(b % a, a)
	        return (g, x - (b // a) * y, y)

	# Returns the inverse of an element a modulo m, given that it exists.
	def modinv(self, a, m):
	    g, x, y = self.egcd(a, m)
	    if g != 1:
	        raise Exception('modular inverse for '+str(a)+' modulo '+str(m)+
	        	' does not exist')
	    else:
	        return x % m

def FiniteRing(n):
# This is fancy. This function returns a class with the modulus hard-coded into it.

	# The class essentially implements the rings Z/nZ, where n is an integer.
	# It also allows multiplication from Z, so it is strictly speaking a
	# commutative Z-algebra.
	class RingElement(object):
		
		# All methods in the form __name__ are special methods recognised
		# by Python. Other methods shouldn't have underscores.
	
		# Sets the value a
		def __init__(self, a):
			self.a = a % n
			
		# Takes two objects of this class, and returns a new object holding
		# the sum, modulo n. Also allows addition with an ordinary int.
		def __add__(self, other):
			if isinstance(other, int):
				return RingElement(self.a + other)
			elif isinstance(other, RingElement):
				return RingElement(self.a + other.a)
			else:
				raise TypeError('Incompatible types')
			
		# As above, but providing commutativity with int addition since Python
		# uses the add method of the first class in x + y.
		def __radd__(self, other):
			if isinstance(other, int):
				return RingElement(self.a + other)
			elif isinstance(other, RingElement):
				return RingElement(self.a + other.a)
			else:
				raise TypeError('Incompatible types')
			
		# Subtraction. As above.
		def __sub__(self, other): 
			if isinstance(other, int):
				return RingElement(self.a - other)
			elif isinstance(other, RingElement):
				return RingElement(self.a - other.a)
			else:
				raise TypeError('Incompatible types')
			
		# Subtraction. As above.
		def __rsub__(self, other): 
			if isinstance(other, int):
				return RingElement(self.a - other)
			elif isinstance(other, RingElement):
				return RingElement(self.a - other.a)
			else:
				raise TypeError('Incompatible types')
			
		# Multiplication. As above.
		def __mul__(self, other): 
			if isinstance(other, int):
				return RingElement(self.a * other)
			elif isinstance(other, RingElement):
				return RingElement(self.a * other.a)
			else:
				raise TypeError('Incompatible types')
			
		# Multiplication. As above.
		def __rmul__(self, other): 
			if isinstance(other, int):
				return RingElement(self.a * other)
			elif isinstance(other, RingElement):
				return RingElement(self.a * other.a)
			else:
				raise TypeError('Incompatible types')
		
		# Division. It checks if gcd(divisor, modulus) = 1, and the
		# multiplies with the inverse of the divisor.
		# Only defined for ring elements.
		def __truediv__(self, other):
			if isinstance(other, RingElement):
				tb = ToolBox()
				try:
					f = self * RingElement(tb.modinv(other.a, n))
					return f
				except:
					raise Exception('Division not well-defined, ' + str(other) + ' is not invertible in ' + self.__name__)
			else:
				raise TypeError('Incompatible types')
		
		# As above. Technical thing.
		def __div__(self, other): 
			if isinstance(other, RingElement):
				tb = ToolBox()
				try:
					f = self * RingElement(tb.modinv(other.a, n))
					return f
				except:
					raise Exception('Division not well-defined, ' + str(other) + ' is not invertible in ' + self.__name__)
			else:
				raise TypeError('Incompatible types')
		
		# Negative element.
		def __neg__(self): 
			return RingElement(-self.a)
		
		# Checks equality. Requires same class membership, and same value.
		def __eq__(self, other): 
			return isinstance(other, RingElement) and self.a == other.a
		
		# Absolute value
		def __abs__(self): 
			return abs(self.a)
		
		# String representation of the value
		def __str__(self): 
			return str(self.a)
		
		# String representation of the object
		def __repr__(self): 
			return '%d (mod %d)' % (self.a, self.n)
			
		# Exponentiation. Implements square-and-multiply.
		def __pow__(self, e):
			# g is the base, e is the exponent, m is the modulus.
		
			# We want to handle all integers e
			if e < 0:
				tb = ToolBox()
				try:
					g = RingElement(1) / self
				except:
					raise Exception('Cannot raise to negative power, ' + repr(self) + ' is not invertible.')
				e = -e
			elif e == 0:
				return RingElement(1)
			else:
				g = self
			
			# Now assume that e > 0.

			# Get binary representation of e, and note that the first two
			# symbols are 0b when e is positive, so we cut them away.
			be = bin(e)[2:]
		
			a = RingElement(1)
			for b in be:
				a = (a * a)
				if b == '1':
					a = (a * g)
			return a
		# This ends the class definition.
	
	# This is in the function again. It sets the parameters of the class.
	RingElement.n = n
	RingElement.__name__ = 'Z/%d' % (n)
	return RingElement


class EllipticCurve(object):
# A new class, implementing elliptic curves.

	def __init__(self, k, a, b):
		# Create the curve. Create a point class.
		# Shorthand: k is a field
		self.k = k
		
		# Make sure the parameters are field elements, not just ints.
		if not isinstance(a, k):
			a = k(a)
		if not isinstance(b, k):
			b = k(b)
		
		# Check the discriminant. Multiplication is quicker than exponentiation.
		d = k(4) * a * a * a + k(27) * b * b
		if d == k(0):
			raise Exception("Curve is singular.")
		
		self.a = a
		self.b = b
		
		# The creator function will return a Point class, of which we can make
		# elements of later.
		class Point(object):
			def __init__(self, x, y):
				# Make sure the parameters are field elements, not just ints.
				if not isinstance(x, k):
					x = k(x)
				if not isinstance(y, k):
					y = k(y)
				
				# Check if the point is on the curve
				r = y * y - x * x * x - self.a * x - self.b
				if r == k(0):
					self.x = x
					self.y = y
				else:
					raise Exception("(%d, %d) is not on y^2 = x^3 + %dx + %d" % (abs(x), abs(y), abs(self.a), abs(self.b)))
			
			# Equality of points
			def __eq__(self, other):
				return isinstance(other, Point) and self.x == other.x and self.y == other.y
			
			# Addition of points.
			def __add__(self, other):
				# There are four cases:
				if isinstance(other, Neutral):
					return self
				elif (self.x == other.x) and (self.y == -other.y):
					return Neutral()
				elif self == other:
					l = (3 * self.x * self.x + self.a) / (2 * self.y)
				else:
					l = (other.y - self.y) / (other.x - self.x)
				
				x = l * l - self.x - other.x
				y = l * (self.x - x) - self.y
				return Point(x, y)
			
			# Subtraction of points. Special case of adding.
			def __sub__(self, other):
				return self + (-other)
				
			# Negation of point. Just reflection across the x-axis.
			def __neg__(self):
				return Point(self.x, -self.y)
				
			# Multiplication with an integer. Uses double-and-add
			# (square-and-multiply), see line 167
			def __mul__(self, e):
				if isinstance(e, Point):
					raise Exception("Multiplication of two points is not defined.")
				# g is the base, e is the exponent, m is the modulus.
		
				# We want to handle all integers e
				if e < 0:
					g = -self
					e = -e
				elif e == 0:
					return Neutral()
				else:
					g = self
			
				# Now assume that e > 0.

				# Get binary representation of e, and note that the first two
				# symbols are 0b when e is positive, so we cut them away.
				be = bin(e)[2:]
		
				a = Neutral()
				for b in be:
					a = (a + a)
					if b == '1':
						a = (a + g)
				return a
				
			# Technical detail to allow n * P, not just P * n
			def __rmul__(self, e):
				return self * e
				
			# String representation of the object.
			def __repr__(self):
				return "(%d, %d) on y^2 = x^3 + %dx + %d" % (abs(self.x), abs(self.y), abs(self.a), abs(self.b))
		
		class Neutral(Point):
			# This class describes the point at infinity.
			def __init__(self):
				return None
						
			def __eq__(self, other):
				return isinstance(other, Neutral)
		
			def __add__(self, other):
				return other
				
			def __neg__(self):
				return self
				
			def __repr__(self):
				return "Point at infinity"
		
		Point.a = a
		Point.b = b
		Point.__name__ = 'y^2 = x^3 + %dx + %d' % (abs(a), abs(b))
		
		# Sets a class
		self.Point = Point
		
		# Sets an object. Note the difference.
		self.Neutral = Neutral()
		
	def __abs__(self):
		# Here goes point counting.
		# TODO -- implement this.
		
		# The general case is hard, which we hope you were able to discover.
		# For the concrete case at hand, this code snippet would do:
		# q = 10000000019
		# lb = q + 1 - sqrt(q)  # Lower bound | sqrt is in the math library,
		# ub = q + 1 + sqrt(q)  # Upper bound | see line 3.
		# for i in range(int(lb), int(ub)+1):
		#	if i % 358573 == 0:
		#		print i
		# (Output: 9999883824)
		
		# See https://github.com/pdinges/python-schoof for an example 
		# implementation of Schoof's algorithm.
		return 1