\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ PARI/GP Script for Elliptic Divisibility Sequences   v. 2.0        \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ This script performs various manipulations related to 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.                                                        \\
\\                                                                    \\
\\ Feel free to distribute this script.                               \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ Set debug = 1 for some information on what is happening
debug = 0;








\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SECTION ONE:  TRANSLATING BETWEEN CURVES AND SEQUENCES   \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ curvetoeds(curve,point)
\\
\\ Given an elliptic curve 'curve' and a point 'point' on that curve,
\\ returns the first four terms of the associated elliptic divisibility
\\ sequences as a row vector.
\\
\\ The curve may be given either in initialized or non-initialized form.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

curvetoeds(curve, point) =

	local(bvals, initvals);

	bvals = vector(8);
	initvals = vector(4);

	\\ the variable 'bvals' stores the usual b_2, b_4, b_6 and b_8
	\\ defined for elliptic curves 
	\\ (see Silverman, Arithmetic of Elliptic Curves, 1986, p. 46)

	bvals[2] = curve[1]^2 + 4 * curve[2];
	bvals[4] = 2*curve[4] + curve[1]*curve[3];
	bvals[6] = curve[3]^2 + 4*curve[5];
	bvals[8] = curve[1]^2*curve[5] + 4*curve[2]*curve[5] 
		- curve[1]*curve[3]*curve[4] + curve[2]*curve[3]^2 
		- curve[4]^2;

	\\ the variable 'initvals' stores the first four terms of the
	\\ sequence; these formulae work for general Weierstrass eqns

	initvals[1] = 1;
	initvals[2] = 2*point[2] + curve[1]*point[1] + curve[3];
	initvals[3] = 3*point[1]^4 + bvals[2]*point[1]^3 
		+ 3*bvals[4]*point[1]^2 + 3*bvals[6]*point[1] 
		+ bvals[8];
	initvals[4] = initvals[2]*( 2*point[1]^6 + bvals[2]*point[1]^5 
		+ 5*bvals[4]*point[1]^4 + 10*bvals[6]*point[1]^3 
		+ 10*bvals[8]*point[1]^2 + (bvals[2]*bvals[8]
		- bvals[4]*bvals[6])*point[1] +bvals[4]*bvals[8] 
		- bvals[6]^2);

	return(initvals);

}






\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edstocurve(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns a curve such that that curve has the point (0,0)
\\ on it and the pair are associated to that sequence.  (I.e. eds
\\ is associated to edstocurve(eds) and point [0,0].)
\\
\\ The curve returned is not initialised.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ This uses formulae due to Christine Swart in her thesis, "Elliptic 
\\ Curves and related sequences", University of London, 2003.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

edstocurve(eds) =

	local(curve,a,b,c);
	curve = vector(5);
	a = eds[2];
	b = eds[3];
	c = eds[4];

	if(eds[1] != 1, print("First term not 1"); return(0); );

	curve[1] = (c+a^5 - 2*a*b)/(a^2*b);
	curve[2] = (a*b^2 + c + a^5 - a*b)/(a^3*b);
	curve[3] = a;
	curve[4] = 1;
	curve[5] = 0;

	return(curve);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ elllong(shortcurve)
\\
\\ Given a vector representing a cubic curve, calculates the standard
\\ values bi, ci, delta and j. (See Silverman, Arithmetic of Elliptic
\\ Curves, 1986, p. 46.)
\\
\\ The vector returned is of the form
\\
\\ [a1,a2,a3,a4,a6,b2,b4,b6,b8,c4,c6,delta,j]
\\ 
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

elllong(shortcurve) =

	local(a1,a2,a3,a4,a6,b2,b4,b6,b8,c4,c6,delta,j);

	a1 = shortcurve[1];
	a2 = shortcurve[2];
	a3 = shortcurve[3];
	a4 = shortcurve[4];
	a6 = shortcurve[5];
	
	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;

	c4 = b2^2 - 24*b4;
	c6 = -b2^3 + 36*b2*b4 - 216*b6;

	delta = -b2^2*b8 - 8*b4^3 - 27*b6^2 + 9*b2*b4*b6;

	if( delta == 0,
		j = "inf";
	,
		j = c4^3/delta;
	);

	return([a1,a2,a3,a4,a6,b2,b4,b6,b8,c4,c6,delta,j]);

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edstocurvefull(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns a curve such that that curve has the point (0,0)
\\ on it and the pair are associated to that sequence.  (I.e. eds
\\ is associated to edstocurve(eds) and point [0,0].)
\\
\\ The vector returned is of the form
\\
\\ [a1,a2,a3,a4,a6,b2,b4,b6,b8,c4,c6,delta,j]
\\
\\ where the constants are as denoted in Silverman, Arithmetic of 
\\ Elliptic Curves, 1986, p. 46
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ This uses formulae due to Christine Swart in her thesis, "Elliptic 
\\ Curves and related sequences", University of London, 2003.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

edstocurvefull(eds) =


	if(eds[1] != 1, print("First term not 1"); return(0); );

	return(elllong(edstocurve(eds)));

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsjinv(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the j-invariant of the associated curve.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
edsjinv(eds) = 

	if(eds[1] != 1, eds = eds/eds[1]; );

	return( edstocurvefull(eds)[13] );

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsdisc(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the discriminant of the associated curve.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
edsdisc(eds) = 

	if(eds[1] != 1, eds = eds/eds[1]; );
	
	return( edstocurve(eds)[12] );

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsc4(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the value c4 of the associated curve.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
edsc4(eds) = 

	if(eds[1] != 1, eds = eds/eds[1]; );
	
	return( edstocurve(eds)[10] );

}




\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsc6(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the value c6 of the associated curve.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
edsc6(eds) = 

	if(eds[1] != 1, eds = eds/eds[1]; );
	
	return( edstocurve(eds)[11] );

}







\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsgtwo(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the value g_2 for the associated curve.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ These use formulae due to M. Ward in his paper "Memoir on Elliptic
\\ Divisibility Sequences."
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edsgtwo(eds)=

	local(a,b,c);
	a = eds[2];
	b = eds[3];
	c = eds[4];

	if(eds[1] != 1, eds = eds/eds[1]; );

	return( 
		(a^20 + 4*a^15*c - 16*a^12*b^3 + 6*a^10*c^2 
		- 8*a^7*b^3*c + 4*a^5*c^3 + 16*a^4*b^6 + 8*a^2*b^3*c^2 
		+ c^4)/(12*a^8*b^4) 
	);

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsgthree(eds)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), with first
\\ term one, returns the value g_3 for the associated curve.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ If first term is not one, sequence is scaled to make it one.
\\
\\ These use formulae due to M. Ward in his paper "Memoir on Elliptic
\\ Divisibility Sequences."
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

edsgthree(eds)=

	local(a,b,c);
	a = eds[2];
	b = eds[3];
	c = eds[4];

	if(eds[1] != 1, eds = eds/eds[1]; );

	return( 
		- (a^30 + 6*a^25*c - 24*a^22*b^3 + 15*a^20*c^2
		- 60*a^17*b^3*c + 20*a^15*c^3 + 120*a^14*b^6 
		- 36*a^12*b^3*c^2 + 15*a^10*c^4 
		- 48*a^9*b^6*c + 12*a^7*b^3*c^3 + 64*a^6*b^9 
		+ 6*a^5*c^5 + 48*a^4*b^6*c^2 + 12*a^2*b^3*c^4 
		+ c^6)/(216*a^12*b^6) 
	);


}














\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SECTION TWO:  CALCULATING VALUES OF SEQUENCES            \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsget(eds,n)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' in the
\\ form of a vector of four terms (1st through 4th terms), and an integer
\\ n, returns the value of the sequence at index n.
\\
\\ Uses Shipsey's method (see edsblockships below).
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edsget(eds,n) = 

	if(debug, print("Calling edsget for W(", n, ")"); );
	
	if( n==0, return(0); );
	if( n < 5 && n > 0, return(eds[n]); );
	if( n < 0, return( -edsget(eds,-n) ); );
	return( edsblockships(eds,n)[5] );

}







\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsgen(eds,len)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' 
\\ (in the form of the first four terms as a vector), and a positive
\\  integer 'len', returns a vector of length 'len' containing the 
\\ first len terms of the sequence.
\\
\\ Uses the recurrence to generate terms one after another linearly.
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{ 

edsgen(eds, len) = 

	local(X, i);

	\\ Stupid to ask it to generate a shorter sequence than 4
	if( len < 4, 
		print("Length too short, returning eds."); 
		return(eds); 
	);

	X = vector(len);

	\\ The first four terms of the resulting vector are 
	\\ already at hand
	X[1]=eds[1];
	X[2]=eds[2];
	X[3]=eds[3];
	X[4]=eds[4];

	\\ For the rest of the terms, use the recurrence relation
	for(i=5,len,

		if( X[i-4] != 0, 
		\\ In this case it is safe to divide by X[i-4]

			X[i] = (X[i-3]*X[i-1]*X[2]^2 
				- X[i-2]^2*X[1]*X[3])
				/X[i-4]/X[1]/X[1];

		, 
		\\ If X[i-4] == 0 then X[i-5] != 0 for 
		\\ non-degenerate sequences
			
			X[i] = (X[i-1]*X[i-4]*X[2]*X[3] 
				- X[i-2]*X[i-3]*X[1]*X[4])
				/X[i-5]/X[2]/X[1];

		);

	);

	return(X);

}






\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsblocklin(eds, len)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' 
\\ (in the form of the first four terms as a vector), and an integer 
\\ 'len', returns a vector of length five containing the terms
\\ len-4 up through len of the sequence.
\\
\\ Uses the recurrence to generate terms one after another linearly.  
\\ Very slow runtime for large len (it is O(len)).
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ Accepts negative len.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ 

edsblocklin(eds, len) = 


	local(X,Xnew,Xbase,i);
	

	X = vector(5);
	Xnew = vector(5);
	Xbase = vector(5);

	\\ 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( len < 0,

		Xnew = edsblocklin(eds, -len+4);
		for(i=1,5,
			X[i] = -Xnew[6-i];
		);

		return(X);

	);

	\\ 'X' stores at each round the most recent block of 
	\\ five terms; start with the terms of eds, 
	\\ plus the zeroth term, which is 0
	X[1]=0;
	X[2]=eds[1];
	X[3]=eds[2];
	X[4]=eds[3];
	X[5]=eds[4];


	\\ if we are requesting a block very close to the origin, 
	\\ the calculation is easy
	if( len == 4, return(X); );
	if( len == 3, return([-X[2],0,X[2],X[3],X[4]]); );
	if( len == 2, return([-X[3],-X[2],0,X[2],X[3]]); );
	if( len == 1, return([-X[4],-X[3],-X[2],0,X[2]]); );
	if( len == 0, return([-X[5],-X[4],-X[3],-X[2],0]); );


	\\ 'Xbase' keeps a permanent record of these first five terms
	\\  for use in the recurrence
	Xbase = X;

	\\ Loop up to requested length, updating X each time to
	\\ represent shifting one term along the sequence
	for(i=1,len-4,

		X = edsblockincrement(X,Xbase);

	);

	return(X);

}







\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsblockships(eds, len)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' 
\\ (in the form of the first four terms as a vector), and an integer 
\\ 'len', returns a vector of length 4 containing the terms
\\ len-4 up through len of the sequence.
\\
\\ Uses Shipsey's double-and-add method to generate the final terms in
\\ O(log (len) ) time.  Very fast.
\\ Her method is described in her thesis (Rachel Shipsey, 
\\ Elliptic Divisibility Sequences, University of London, 2000)
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\
\\ Accepts negative len.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ 

edsblockships(eds, len) = 

	local(X, Xnew, Xbase, Xfinal, i, j, doubleaddchainlength,
		 doubleaddchain, scalefactor, remainder);


	\\ 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( len < 0,
		X = vector(5);
		Xnew = vector(5);

		Xnew = edsblockships(eds, -len+4);
		for(i=1,5,
			X[i] = -Xnew[6-i];
		);

		return(X);

	);

	\\ if we are requesting a block very close to the origin,
	\\ the calculation is easy
	if( len == 4, return([0,eds[1],eds[2],eds[3],eds[4]]); );
	if( len == 3, return([-eds[1],0,eds[1],eds[2],eds[3]]); );
	if( len == 2, return([-eds[2],-eds[1],0,eds[1],eds[2]]); );
	if( len == 1, return([-eds[3],-eds[2],-eds[1],0,eds[1]]); );
	if( len == 0, return([-eds[4],-eds[3],-eds[2],-eds[1],0]); );



	\\ 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 = vector(8);
	Xnew = vector(8);
	Xbase = vector(4);


	\\ if eds doesn't have first term one, we have to scale it
	\\ to one that does and then un-scale at the end
	scalefactor = eds[1];
	for(i=1,4,
		eds[i] = eds[i]/scalefactor;
	);

	\\to make sure final block contains len-4 up through len
	len = len-1;


	\\ the initial terms of the sequence are stored in 'Xbase'
	Xbase[1]=eds[1];
	Xbase[2]=eds[2];
	Xbase[3]=eds[3];
	Xbase[4]=eds[4];
	if( eds[2] == 0,
		invtwo = "inf";
	,
		\\ the inverse of the second term is precomputed
		invtwo = 1/eds[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] = -eds[2];
	X[2] = -eds[1];
	X[3] = 0;
	X[4] = eds[1];
	X[5] = eds[2];
	X[6] = eds[3];
	X[7] = eds[4];
	X[8] = eds[2]^3*eds[4] - eds[3]^2*eds[1]*eds[3];

	\\ The length of the double-add chain is computed
	doubleaddchainlength = ceil(log(len+1)/log(2));
	
	\\ The chain will be stored in this vector
	doubleaddchain = vector(doubleaddchainlength);

	\\ compute the double-add-chain
	remainder = len;
	i = 0;
	while(remainder !=0,
		bit = lift(Mod(remainder,2));
		doubleaddchain[doubleaddchainlength - i] = bit;
		remainder = (remainder-bit)/2;
		i = i+1;
	);

	\\ Show the double-add-chain if debug data is on
	if( debug == 1, 
		print("Double and add chain is: " doubleaddchain);
	 );

	\\ Loop through the double-add-chain doubling or adding
	for(j=2,doubleaddchainlength,

		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 != "inf",
				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;
			,
				Xnew[2] = 0;
				Xnew[4] = 0;
				Xnew[6] = 0;
				Xnew[8] = 0;
			);
			X = Xnew;
		,

			\\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 != "inf",
				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;
			,
				Xnew[1] = 0;
				Xnew[3] = 0;
				Xnew[5] = 0;
				Xnew[7] = 0;
			);
			X = Xnew;
		);

	);

	\\ Put the final data into a vector of length five and return
	Xfinal = vector(5);
	Xfinal[1] = X[1]*scalefactor;
	Xfinal[2] = X[2]*scalefactor;
	Xfinal[3] = X[3]*scalefactor;
	Xfinal[4] = X[4]*scalefactor;
	Xfinal[5] = X[5]*scalefactor;

	return(Xfinal);

}






\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsblockincrement(block,base)
\\
\\ Given an elliptic divisibility sequence 'base' in the form of a 
\\ vector of five terms (0th through 4th terms) or four terms 
\\ (1st through 4th), and a block
\\ of length 5 in that elliptic divisibility sequence or a Somos
\\ sequence translated from it,
\\ it calculates the block shifted one to the right.  The sequence
\\ must be non-degenerate.
\\
\\ Used as an internal function for edsblocklin and others.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edsblockincrement(block,base)=

	local(newblock);

	newblock = vector(5);

	\\ if base is a block of length four, make it length five
	if( length(base) == 4, 
		base = [0,base[1],base[2],base[3],base[4]]; 
	);

	\\ calculate newest term and store in 'newblock'
	if( block[2] != 0, \\ safe to divide by block[2]

		newblock[5] = (block[3]*block[5]*base[3]^2 
					- block[4]^2*base[2]*base[4])
					/block[2]/base[2]/base[2];

	, 
		\\ otherwise it is safe to divide by block[1] 
		\\(both cannot be zero in a non-degenerate sequence)

		newblock[5] = (block[5]*block[2]*base[3]*base[4] 
				- block[4]*block[3]*base[2]*base[5])
				/block[1]/base[3]/base[2];

	);

	\\ shift previous terms and store in 'newblock'
	newblock[4] = block[5];
	newblock[3] = block[4];
	newblock[2] = block[3];
	newblock[1] = block[2];

	return(newblock);


}










\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SECTION THREE:  TRANSFORMATIONS AND PROPERTIES OF SEQUENCES  \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsrankofapp(eds, len)
\\
\\ Given an elliptic divisibility sequence 'eds' 
\\ (in the form of the first four terms as a vector), and a positive
\\  integer 'len', returns the index of the first positive-indexed zero
\\  term before index len+1 or 0 if none is found up to that distance
\\
\\ If 'len' is omitted or is zero, will run until it finds a zero or 
\\ is interrupted.
\\
\\ Uses a linear method to calculate terms one after another until
\\ a zero is found.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ 

edsrankofapp(eds, len) = 

	local(X, Xnew, Xbase, i);

	X = vector(5);
	Xnew = vector(5);
	Xbase = vector(5);

	\\ 'X' stores at each round the most recent block of five terms
	\\ Start with the terms of eds, plus the zeroth term, 
	\\ which is 0.  Return index if a zero is found.
	X[1]=0;
	for(i=1,4,
		X[i+1]=eds[i];
		if( X[i+1] == 0, return(i) );
	);

	\\ 'Xbase' keeps a permanent record of these first five terms
	\\ for use in the recurrence
	Xbase = X;

	if( len > 1, 
		\\ loop up to requested length looking for a zero
		for(i=5,len,
			X = edsblockincrement(X,Xbase);
			if( X[5] == 0, return(i); );
		);
	,
		\\ loop up to requested length looking for a zero
		i = 4;
		while(X[5] != 0,
			X = edsblockincrement(X,Xbase);
			i = i+1;
		);
		return(i);

	);

	return(0);

}







\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsperiod(eds, len)
\\
\\ Given a non-degenerate elliptic divisibility sequence 'eds' 
\\ (in the form of the first four terms as a vector), and a positive
\\  integer 'len', returns the period of the sequence if it is less 
\\ than len+1, or 0 if the period is not found up to that distance.
\\
\\ Uses a linear method to calculate terms one after another until
\\ the period is found.
\\
\\ 'Non-degerate' means the first, second and third terms are non-zero.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ 

edsperiod(eds, len) = 

	local(X, Xnew, Xbase, i);

	X = vector(5);
	Xnew = vector(5);
	Xbase = vector(5);

	\\ 'X' stores at each round the most recent block of five terms
	\\ Start with the terms of eds, plus the zeroth term, 
	\\ which is 0
	X[1]=0;
	for(i=1,4,
		X[i+1]=eds[i];
	);

	\\ 'Xbase' keeps a permanent record of these first five terms 
	\\ for use in the recurrence
	Xbase = X;

	\\ loop up to requested length searching for a zero
	if( len > 0,
		for(i=1,len,
			X = edsblockincrement(X,Xbase);
			if( X == Xbase, return(i); );
		);
	,
		i = 1;
			X = edsblockincrement(X,Xbase);
		while(X != Xbase,
			X = edsblockincrement(X,Xbase);
			i = i+1;
		);
		return(i);

	);
	return(0);

}









\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edssubseq(eds, n)
\\
\\ Given an elliptic divisibility sequence 'eds' associated to P on E
\\ (in the form of the first four terms as a vector), and a positive
\\ integer 'n', the sequence associated to [n]P on E.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edssubseq(eds,n) = 

	local(neweds, nterm);

	nterm = edsget(eds,n);
	neweds = vector(4);

	neweds[1] = 1;
	neweds[2] = edsget(eds,2*n)/nterm^4;
	neweds[3] = edsget(eds,3*n)/nterm^9;
	neweds[4] = edsget(eds,4*n)/nterm^16;

	return(neweds);

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsmakeint(eds)
\\
\\ Given an elliptic divisibility sequence 'eds' associated to P on E
\\ of type "t_FRAC",
\\ return an eds equivalent to it but with integer values, and no 
\\ unnecessary gcd.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edsmakeint(eds) = 

	local(neweds, factors, factr,gg);

	neweds = eds;
	factr = 1;

	gg = gcd(numerator(neweds[2]),numerator(neweds[3]));
	while( gg != 1,
		neweds = edsequiv(neweds,g^(-1));		
		gg = gcd(numerator(neweds[2]),numerator(neweds[3]));
	);

	while( type(neweds[4]) != "t_INT",
		factors = factor(denominator(neweds[4]));
		factr = 1;
		for(i=1,matsize(factors)[1],
			factr = factr*factors[i,1];
		);
		neweds = edsequiv(neweds,factr);
	);


	return(neweds);

}






\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsequiv(eds, factr)
\\
\\ Given an elliptic divisibility sequence 'eds'
\\ return an eds equivalent to it by factor 'factor'^(n^2-1)
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edsequiv(eds, factr) = 

	local(neweds);

	neweds = vector(4);
	neweds[1] = eds[1];

	for(i=2,4,
		neweds[i] = eds[i]*factr^(i^2-1);
	);

	return(neweds);

}




\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edsperfper(eds, stop)
\\
\\ Given an elliptic divisibility sequence 'eds'
\\ return an eds equivalent to it which has rank of apparition equal
\\ to its period, if one is found before stop is reached.  (The
\\ algorithm loops through equivalences by i=1 to stop.  If stop=0,
\\ it will go until one is found (potentially never).
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\


{

edsperfper(eds,stop) =

	local(i,r);

	r = edsrankofapp(eds);
	if(debug, print("Rank of apparition is: ", r););

	i=1;
	while(edsget(edsequiv(eds,i),r+1) != 1 && i != stop,
		i = i+1;		
	);

	if( i == stop,
		return(0);
	);

	return(edsequiv(eds,i));

}