fruit.py

N=4

def gcd(a, b):
	while b != 0:
		(a, b) = (b, a % b)
	return a

class Q(object):
	def __init__(self, p, q=1):
		d = gcd(p, q)
		self.p = p // d
		self.q = q // d

	def __add__(self, other):
		return Q(self.p * other.q + other.p * self.q, self.q * other.q)

	def __sub__(self, other):
		return Q(self.p * other.q - other.p * self.q, self.q * other.q)

	def __neg__(self):
		return Q(-self.p, self.q)

	def __mul__(self, other):
		return Q(self.p * other.p, self.q * other.q)

	def __truediv__(self, other):
		return Q(self.p * other.q, self.q * other.p)

	def __repr__(self):
		return "%d / %d" % (self.p, self.q)


class P(object):
	def __init__(self, ec, x, y):
		self.ec = ec
		self.x = x
		self.y = y

	# From the paper: recover a, b, and c from x and y
	def abc(self):
		aa = Q(8*(N+3)) - self.x + self.y
		bb = Q(8*(N+3)) - self.x - self.y
		cc = Q(-8*(N+3)) - Q(2*(N+2))*self.x
		q = aa.q * bb.q * cc.q
		a = (aa * Q(q)).p
		b = (bb * Q(q)).p
		c = (cc * Q(q)).p
		d = gcd(gcd(a, b), c)
		return (a // d, b // d, c // d)

	def check(self):
		(a, b, c) = self.abc()
		(a, b, c) = (Q(a), Q(b), Q(c))
		return a/(b+c) + b/(a+c) + c/(a+b)

	# Elliptic curve algebraic operations taken from 
	# https://www.math.brown.edu/johsilve/Presentations/WyomingEllipticCurve.pdf

	# Implements P1 + P2 by drawing a line through 2 points and finding where it intersects the curve
	def sec(self, other):
		m = (other.y - self.y) / (other.x - self.x)
		b = self.y - self.x * m
		x3 = m*m - Q(ec.A) - self.x - other.x
		y3 = m * x3 + b
		return P(self.ec, x3, -y3)

	# Implements P + P by drawing a tangent line at P and finding where it intersects the curve
	def tan(self):
		m = (self.x * self.x * Q(3) + self.x * Q(ec.A * 2) + Q(ec.B)) / (self.y * Q(2))
		b = self.y - self.x * m
		x3 = m*m - Q(ec.A) - self.x - self.x
		y3 = m * x3 + b
		return P(self.ec, x3, -y3)

	def __repr__(self):
		return "(%s, %s)" % (self.x, self.y)


class EC(object):
	def __init__(self, A, B):
		self.A = A
		self.B = B

	def pt(self, x, y):
		return P(self, Q(x), Q(y))

# From the paper: Cubic model of the elliptic curve
ec = EC(4*N*N + 12*N - 3, 32*(N+3))
# Fom the paper: the point G
g = ec.pt(-4, 28)

# some other fun points to play around with
o = ec.pt(0, 0)
g2 = g.sec(o)
g3 = g2.tan()

# From the paper: 9G is the first positive solution
g8 = g.tan().tan().tan()
g9 = g8.sec(g)
(a, b, c) = g9.abc()
g9.check()