#!/usr/bin/env python -i

# 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):
	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)

	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):
	class RingElement(object):
		def __init__(self, a):
			self.a = a % n
			
		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')
			
		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')
			
		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')
			
		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')
			
		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')
			
		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')
		
		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')
		
		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')
		
		def __neg__(self): 
			return RingElement(-self.a)
		
		def __eq__(self, other): 
			return isinstance(other, RingElement) and self.a == other.a
		
		def __abs__(self): 
			return abs(self.a)
		
		def __str__(self): 
			return str(self.a)
		
		def __repr__(self): 
			return '%d (mod %d)' % (self.a, self.n)
			
		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
	
	RingElement.n = n
	RingElement.__name__ = 'Z/%d' % (n)
	return RingElement

class EllipticCurve(object):
	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
		
		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)
				
				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)))
			
			def __eq__(self, other):
				return isinstance(other, Point) and self.x == other.x and self.y == other.y
			
			def __add__(self, other):
				# TODO: You must implement this.
				
				# There are four cases:
				#if isinstance(other, Neutral):
				
				#elif (self.x == other.x) and (self.y == -other.y):
				
				#elif self == other:
				
				#else:
				
				
				#x = 
				#y = 
				#return Point(x, y)
				
				# Remove this line when you're done.
				return 0
			
			def __sub__(self, other):
				return self + (-other)
				
			def __neg__(self):
				return Point(self.x, -self.y)
				
			def __mul__(self, e):
				# TODO: You must implement this.
			
				# Make sure that the other operand is not a point.
				if isinstance(e, Point):
					raise Exception("Multiplication of two points is not defined.")
				
				# Recall that point multiplication corresponds to exponentiation, so
				# you can use the same idea as in lines 134-161.
						
				# Remove this line when you're done.
				return 0
				
			def __rmul__(self, e):
				return self * e
				
			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))
		self.Point = Point
		self.Neutral = Neutral()
		
		# Check if the curve is singular. If yes, raise exception.
		
	def __abs__(self):
		# Here goes point counting.
		# TODO: You must implement this.
		
		# There are several strategies you can try. Listing and counting
		# all the points is only feasible if the curve is small.
		
		# Better idea: Find some points on the curve and their orders.
		# Use this along with Hasse's bound, so that you can narrow the
		# possible candidates.
		
		# The ambitious might want to implement Shoof's algorithm.
		
		# Remove this line when you're done
		return 0
		
print("""
Usage: python -i ellipticcurve.py
   or: ./ellipticcurve.py (if you have set the executable flag)
	
The file provides the classes FinteRing and EllipticCurve. FiniteRing
is an implementation of computations modulo some n, and is a neat way
to not having to care about bringing the modulus around for other 
computations.

Syntax:

	ring = FiniteRing(n)
	element = ring(value)

Example:
	
	F71 = FinteRing(71)
	a = F71(70)
	b = F71(2)
	a + b  # Gives 1 (mod 71)
	a * b  # Gives 69 (mod 71)
	
Next, EllipticCurve should contain the algorithms used for elliptic
curves. You need to implement __add__, __mult__ and __abs__ yourself.

__add__ corresponds to the plus operation, allowing you to compute

	P + Q
	
__mult__ is multiplication,

	n * P, 		where n is an integers
	
__abs__ is |E|, so we interpret it as cardinality here.

	abs(E)

Syntax:

	curve = EllipticCurve(field, a, b)
	point = curve.Point(x, y)

where a, b, x and y should be field elements.

Example:
	E = EllipticCurve(F71, F71(1), F71(28))
	P = E.Point(F71(1), F71(32))
	O = E.Neutral 	# E.Neutral is the point at infinity
""")