#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
#\\ SAGE Class for Elliptic Divisibility Sequences   v. 1.0            \\
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
#\\ Katherine E. Stange                                                \\
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
#\\ This is a class for use in SAGE which implements elliptic          \\
#\\ divisibility sequences.                                            \\
#\\                                                                    \\
#\\ See http://www.math.brown.edu/~stange/                             \\
#\\ or contact <stange at math dot brown dot edu> for info.            \\
#\\                                                                    \\
#\\ Throughout this script, it is assumed the sequences take values    \\
#\\ in a field (before division, checking if something is zero is all  \\
#\\ we do).  Much will work in general rings, but no promises.         \\
#\\                                                                    \\
#\\ Feel free to distribute this script.  If you alter it, add a       \\
#\\ description of the alteration.  Acknowledge my original version.   \\
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


class EDS:
	"""An elliptic divisibility sequence."""

	##########################################
	### The initialisation and representation
	##########################################

	def __init__(self, *indata):
		if len(indata) == 1:
			print "Creating an EDS with a list of first terms"
			self.initial = indata[0]
			if self.initial[0] != 1:
				print "The first term isn't 1: error"
			self.weierstrass = [0,0,0,0,0]
			a = self.initial[1]; b = self.initial[2]; c = self.initial[3]
			if a == 0 or b == 0 or c == 0:
				print "It's a degenerate sequence!  We can't handle it."
			self.weierstrass[0] = (c+a^5 - 2*a*b)/(a^2*b)
			self.weierstrass[1] = (a*b^2 + c + a^5 - a*b)/(a^3*b)
			self.weierstrass[2] = a; self.weierstrass[3] = 1; self.weierstrass[4] = 0
			self.pt = [0,0]
		elif len(indata) == 2:
			##print "Creating an EDS from a Weierstrass equation and point"
			if type(indata[0]) == list:
				##print "The Weierstrass equation is just a list of coefficients"
				a1 = indata[0][0]; a2 = indata[0][1]; a3 = indata[0][2]
				a4 = indata[0][3]; a6 = indata[0][4]
				b2 = a1^2 + 4*a2; b4 = 2*a4 + a1*a3; b6 = a3^2 + 4*a6
				b8 = a1^2*a6 + 4*a2*a6 - a1*a3*a4 + a2*a3^2 - a4^2
				pt = indata[1]
			else:
				##print "The Weierstrass equation is actually an EllipticCurve (I hope)"
				cur = indata[0]; pt = indata[1]; a1 = cur.a1(); a3 = cur.a3()
				a2 = cur.a2(); a4 = cur.a4(); a6 = cur.a6()
				b2 = cur.b2(); b4 = cur.b4(); b6 = cur.b6(); b8 = cur.b8()
			self.weierstrass = [a1, a2, a3, a4, a6]
			self.pt = pt
			self.initial = [1,0,0,0]
			self.initial[1] = 2*pt[1] + a1*pt[0] + a3
			self.initial[2] = 3*pt[0]^4 + b2*pt[0]^3 + 3*b4*pt[0]^2 + 3*b6*pt[0] \
				+ b8
			self.initial[3] = self.initial[1]*( 2*pt[0]^6 + b2*pt[0]^5 + 5*b4*pt[0]^4 \
				+ 10*b6*pt[0]^3 + 10*b8*pt[0]^2 + (b2*b8- b4*b6)*pt[0] \
				+ b4*b8 - b6^2)
		else:
			print "Not valid input"
		# Put in place various other data
		self.bs = [0,0,0,0]; self.cs = [0,0]
		self.bs[0] = self.weierstrass[0]^2 + 4*self.weierstrass[1]
		self.bs[1] = 2*self.weierstrass[3] + self.weierstrass[0]*self.weierstrass[2]
		self.bs[2] = self.weierstrass[2]^2 + 4*self.weierstrass[4]
		self.bs[3] = self.weierstrass[0]^2*self.weierstrass[4] + 4*self.weierstrass[1]*self.weierstrass[4] \
			- self.weierstrass[0]*self.weierstrass[2]*self.weierstrass[3] \
			+ self.weierstrass[1]*self.weierstrass[2]^2 - self.weierstrass[3]^2
		self.cs[0] = self.bs[0]^2 - 24*self.bs[1];
		self.cs[1] = -self.bs[0]^3 + 36*self.bs[0]*self.bs[1] - 216*self.bs[2]
		self.delta = -self.bs[0]^2*self.bs[3] - 8*self.bs[1]^3 - 27*self.bs[2]^2 \
			+ 9*self.bs[0]*self.bs[1]*self.bs[2]
		if self.delta == 0:
			self.j = "inf"
		else:
			self.j = self.cs[0]^3/self.delta
		self.termlist = [0, self.initial[0], self.initial[1], self.initial[2], self.initial[3]]

	def __repr__(self):
		return "An elliptic divisibility sequence with first terms %s, %s, %s, %s" % (self.term(1), self.term(2), \
			self.term(3), self.term(4))

	####################################
	### The following are user functions
	####################################

	def curve(self):
		### Returns the Weierstrass form [a1, a2, a3, a4, a6] of a cubic curve associated to the sequence
		return self.weierstrass

	def point(self):
		### Returns the point on self.curve() which is associated to the sequence
		return self.pt

	def delta(self):
		### Returns the discriminant of the sequence
		return self.delta

	def j(self):
		### Returns the j-invariant of the sequence
		return self.j

	def a(self,n):
		### Returns the n-th coefficient a_n of self.curve() (n = 1, 2, 3, 4, or 6)
		return self.weierstrass[n-1]

	def b(self,n):
		### Returns the b_n of self.curve() (n = 2, 4, 6, or 8)
		return self.bs[n/2-1]

	def c(self,n):
		### Returns the c_n of self.curve() (n = 4 or 6)
		return self.cs[n/2-2]

	def term(self,n):
		### Returns the n-th term of the sequence
		if n < 0:
			return -self.term(-n)
		elif n < 100:
			self.expand_termlist(n)
			return self.termlist[n]
		else:
			return self.shipsey_block(n)[5]

	def terms(self,n):
		### Returns the first n terms of the sequence in a list
		if n < 1:
			return "only takes positives"
		self.expand_termlist(n)
		if len(self.termlist) > n+1:
			temp_termlist = self.termlist[0:n]
		else:
			temp_termlist = self.termlist
		return temp_termlist

	def rank(self,n):
		### Returns the rank of the sequence, if it is not greater than n; otherwise an error message
		self.expand_termlist(n)
		for i in range(1,n+1):
			if self.termlist[i] == 0:
				return i
		return "Rank not found"

	def period(self,n):
		## Returns the period of the sequence, if it is not greater than n^2; otherwise an error message
		r = self.rank(n)
		for j in range(1,r+1):
			if (self.term(j*r+1) == 1) and (self.term(j*r+2) == self.term(2)) \
				and (self.term(j*r+3) == self.term(3)) and (self.term(j*r+4) == self.term(4)):
				return j*r
		return "Period not found"

	def subsequence(self,n): 
		### Returns the subsequence W_{nk}, k in Z, as an EDS object
		nterm = self.term(n)
		return EDS([1, self.term(2*n)/nterm^4, self.term(3*n)/nterm^9, self.term(4*n)/nterm^16])

	def equivalent(self, f):
		### Returns the sequence W_n*f^{n^2-1} as an EDS object
		return EDS([1, self.term(2)*f^3, self.term(3)*f^8, self.term(4)*f^15])

	##################################################
	### The following are meant as internal functions
	##################################################

	def expand_termlist(self, n):
		m = len(self.termlist)
		if m <= n:
			self.termlist = self.termlist+[0]*(n-m+1)
		for i in range(m, n+1):
			if self.termlist[i-4] != 0: 
			# In this case it is safe to divide by X[i-4]
				self.termlist[i] = (self.termlist[i-3]*self.termlist[i-1]*self.termlist[2]^2 \
					- self.termlist[i-2]^2*self.termlist[1]*self.termlist[3]) \
					/ ( self.termlist[i-4]*self.termlist[1]^2 )
			else: 
			# If X[i-4] == 0 then X[i-5] != 0 for non-degenerate sequences
				self.termlist[i] = (self.termlist[i-1]*self.termlist[i-4]*self.termlist[2] \
					*self.termlist[3]- self.termlist[i-2]*self.termlist[i-3] \
					*self.termlist[1]*self.termlist[4]) \
					/ (self.termlist[i-5]* self.termlist[2] * self.termlist[1] )

	def shipsey_block(self, l):
		# elliptic divisibility sequences satisfy an
		# antisymmetry property, so we calculate out in the
		# positive direction and then reverse the block and put
		# on minus signs
		if l < 0:
			X = [0]*6
			Xnew = self.shipsey_block(-l+4)
			for i in range(1,6):
				X[i] = -Xnew[6-i]
			return X

		# if we are requesting a block within termlist, it's pretty easy
		if l < len(self.termlist):
			return [self.term(l-4), self.term(l-3), self.term(l-2), self.term(l-1), self.term(l)]
	
		# Throughout, we use blocks (segments) of length 8 of the sequence.
		# A block is 'based on' k if it represents terms k-3 up through k+4 of the sequence
		X = [0]*9; Xnew = [0]*9; Xbase = [0]*5

		# to make sure final block contains len-4 up through len
		l = l-1

		# the initial terms of the sequence are stored in 'Xbase'
		Xbase[1]=self.term(1); Xbase[2]=self.term(2)
		Xbase[3]=self.term(3); Xbase[4]=self.term(4)

		# the inverse of the 2nd term is precomputed
		if self.term(2) == 0:
			invtwo = 0
		else:
			invtwo = 1/self.term(2)

		#\\ 'X' starts as a block based on 1 
		#	\\ (i.e. a block of length 8 whose 4th term is 
		#	\\ the first term of the sequence)
		X[1] = self.term(-2)
		X[2] = self.term(-1)
		X[3] = self.term(0)
		X[4] = self.term(1)
		X[5] = self.term(2)
		X[6] = self.term(3)
		X[7] = self.term(4)
		X[8] = self.term(5)
		Xnew = [0]*9

		# The length of the double-add chain is computed
		doubleaddchainlength = ceil(log(l+1)/log(2))
	
		# The chain will be stored in this vector
		doubleaddchain = [0]*(doubleaddchainlength+1)

		# compute the double-add-chain
		remainder = l; i = 0
		while remainder !=0:
			bit = remainder % 2
			doubleaddchain[doubleaddchainlength - i] = bit
			remainder = (remainder-bit)/2
			i = i+1
		
		# Loop through the double-add-chain doubling or adding
		for j in range(2,doubleaddchainlength+1):
			if doubleaddchain[j] == 0:
				# We double the block 
				# (get block based on 2*k from block
				# based on k)
				Xnew[1] = X[4]*X[2]^3 - X[1]*X[3]^3
				Xnew[3] = X[5]*X[3]^3 - X[2]*X[4]^3
				Xnew[5] = X[6]*X[4]^3 - X[3]*X[5]^3
				Xnew[7] = X[7]*X[5]^3 - X[4]*X[6]^3
				if invtwo != 0:
					Xnew[2] = X[3]*(X[5]*X[2]^2-X[1]*X[4]^2)*invtwo
					Xnew[4] = X[4]*(X[6]*X[3]^2-X[2]*X[5]^2)*invtwo
					Xnew[6] = X[5]*(X[7]*X[4]^2-X[3]*X[6]^2)*invtwo
					Xnew[8] = X[6]*(X[8]*X[5]^2-X[4]*X[7]^2)*invtwo
				else:
					Xnew[2] = 0; Xnew[4] = 0
					Xnew[6] = 0; Xnew[8] = 0
				X = Xnew[:]
			else:
				# We double-add the block 
				# (get block based on 2*k+1 from block 
				# based on k)
				Xnew[2] = X[5]*X[3]^3 - X[2]*X[4]^3
				Xnew[4] = X[6]*X[4]^3 - X[3]*X[5]^3
				Xnew[6] = X[7]*X[5]^3 - X[4]*X[6]^3
				Xnew[8] = X[8]*X[6]^3 - X[5]*X[7]^3
				if invtwo != 0:
					Xnew[1] = X[3]*(X[5]*X[2]^2 - X[1]*X[4]^2)*invtwo
					Xnew[3] = X[4]*(X[6]*X[3]^2 - X[2]*X[5]^2)*invtwo
					Xnew[5] = X[5]*(X[7]*X[4]^2 - X[3]*X[6]^2)*invtwo
					Xnew[7] = X[6]*(X[8]*X[5]^2 - X[4]*X[7]^2)*invtwo
				else:
					Xnew[1] = 0; Xnew[3] = 0
					Xnew[5] = 0; Xnew[7] = 0
				X = Xnew[:]
		return X