\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ PARI/GP Script for Rank Two Elliptic Nets            v. 1.1        \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ This script performs various manipulations related to elliptic     \\
\\ nets.  It requires the script for elliptic divisibility sequences  \\
\\ edstools.gp version 2.0.                                           \\
\\                                                                    \\
\\ Version 1.1 corrects a minor error in the function                 \\
\\ twonettocurvechar2().                                              \\   
\\                                                                    \\
\\ 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 nets take values         \\
\\ in a field.  Sometimes this field is required to have              \\
\\ characteristic not equal to two.  Many things will work            \\
\\ for general rings, but no guarantees there or anywhere.            \\
\\ Zero divisors in particular are a big problem.                     \\
\\                                                                    \\
\\ The functions in this script are restricted to rank 2.             \\
\\                                                                    \\
\\ A rank two elliptic net is represented as a vector of four         \\
\\ entries:                                                           \\
\\                                                                    \\
\\ [ W(2,0), W(0,2), W(2,1), W(1,2) ]                                 \\
\\                                                                    \\
\\ and is called 'non-degenerate' if none of the following occurs     \\
\\ 1) W(2,1) = W(1,2)                                                 \\
\\ 2) W(2,0) = W(2,1) = 0                                             \\
\\ 3) W(0,2) = W(1,2) = 0                                             \\
\\ Degenerate curves are generally not allowed anywhere.              \\
\\                                                                    \\
\\ Feel free to distribute this script.  If you alter it, add a       \\
\\ description of the alteration.  Acknowledge my original version.   \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\r edstools

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





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





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ curvetotwonet(curve,pointa,pointb)
\\
\\ Given an elliptic curve 'curve' and a points 'pointa' and 'pointb'
\\ on that curve,
\\ returns the terms (2,0), (0,2), (2,1), (1,2) of the net in that
\\ order as a vector of length four.
\\
\\ Will return '0' in the case that pointa, pointb, pointa+pointb or 
\\ pointa-pointb is the zero point on the curve.
\\
\\ The curve may be given either in initialised or non-initialised form.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

curvetotwonet(curve, pointa, pointb) = 

	local(initvals);

	\\ The vector initvals will hold the terms
	\\ (2,0), (0,2), (2,1), (1,2)	
	initvals = vector(4); 


	\\ It makes no sense to form a net with P, Q, P+Q or P-Q trivial
	if( pointa == [0] || pointb == [0],
		if(debug, 
			print("curvetotwonet does not accept zero
				 points"); 
		);
		return(0);
	);

	if( pointa[1] == pointb[1],
		if(debug, 
			print("curvetotwonet does not accept points
				 which are equal or inverses"); 
		);
		return(0);
	);


	\\ Formulae for net polynomials gives the initial values
	initvals[1] = 2*pointa[2] + curve[1]*pointa[1] + curve[3];
	initvals[2] = 2*pointb[2] + curve[1]*pointb[1] + curve[3];
	initvals[3] = 2*pointa[1] + pointb[1] - ((pointb[2] 
		- pointa[2])/(pointb[1] - pointa[1]))^2 
		- curve[1]*((pointb[2] - pointa[2])/(pointb[1] 
		- pointa[1])) + curve[2];
	initvals[4] = 2*pointb[1] + pointa[1] - ((pointb[2] 
		- pointa[2])/(pointb[1] - pointa[1]))^2 
		- curve[1]*((pointb[2] - pointa[2])/(pointb[1] 
		- pointa[1])) + curve[2];

	\\ Error trap: the resulting net is degenerate, which
	\\ should only happen
	\\ if P+Q or P-Q is [0] (trapped earlier)
	if(initvals[3] == initvals[4],
		if(debug, print("The resulting net is degenerate.
			  This should have been caught earlier."); );
		print("ERROR, PLEASE REPORT THIS BUG 4301982357
			 to stange@math.brown.edu");
	);

	return(initvals);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonettocurve(net)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of length four (terms (2,0), (0,2), (2,1), (1,2)),
\\ returns a vector of length two whose first component is a length
\\ five vector representing a curve, and whose second component is
\\ a value 'x' such that that curve has the point (x,0) and (-x,0)
\\ on it and the triple are associated to that net.  (I.e. net
\\ is associated to curve twonettocurve(net)[1] and points [x,0]
\\ and [-x,0].)
\\
\\ The curve returned is not initialised.
\\
\\ The field must be of characteristic not equal to two.
\\
\\ This uses formulae from my thesis "Elliptic Nets and Elliptic 
\\ Curves."
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonettocurve(net) =

	local(curve,n20,n02,n21,n12,xcoord, returnvector);

	\\ set up return vector (curve and x coordinate)
	returnvector = vector(2);
	curve = vector(5);

	\\ Return zero if the net is degenerate.
	if( istwonetdegen(net),
		if(debug, print("Degenerate nets taste bad!"); );
		return(0);
	);

	\\ gives names to net vals for ease of formulae
	n20 = net[1];
	n02 = net[2];
	n21 = net[3];
	n12 = net[4];


	\\ Coefficients of curve
	curve[1] = (n20 - n02)/(n21 - n12);
	curve[2] = (n21 + n12)/2;
	curve[3] = (n20 + n02)/2;
	curve[4] = -(n21 - n12)^2/4;
	curve[5] = -(n21 - n12)^2*(n21 + n12)/8;

	\\ X coordinate of first point (and negative of second)
	xcoord = (n21 - n12)/2;

	\\ Create vector to return (curve and coordinate)
	returnvector = [curve, xcoord];

	return(returnvector);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonettocurvechar2(net)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of length four (terms (2,0), (0,2), (2,1), (1,2)),
\\ returns a vector of length two whose first component is a length
\\ five vector representing a curve, and whose second component is
\\ a value 'x' such that that curve has the point (0,0) and (x,0)
\\ on it and the triple are associated to that net.  (I.e. net
\\ is associated to curve twonettocurve(net)[1] and points [0,0]
\\ and [x,0].)
\\
\\ The curve returned is not initialised.
\\
\\ The field may be of characteristic equal to two (or may not).
\\
\\ This uses formulae from my thesis "Elliptic Nets and Elliptic 
\\ Curves."
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonettocurvechar2(net) =

	local(curve,n20,n02,n21,n12,xcoord, returnvector);

	\\ Set up return vector (curve and x coordinate)
	returnvector = vector(2);
	curve = vector(5);

	\\ Catch degenerate nets (return zero)
	if( istwonetdegen(net),
		if(debug, print("Degenerate nets taste bad!"); );
		return(0);
	);

	\\ Give names to variables for ease of formulae
	n20 = net[1];
	n02 = net[2];
	n21 = net[3];
	n12 = net[4];

	\\ Coefficients of curve
	curve[1] = (n20 - n02)/(n21 - n12);
	curve[2] = 2*n21 - n12;
	curve[3] = (n20);
	curve[4] = (n21 - n12)*n21;
	curve[5] = 0;

	\\ X coordinate of second point
	xcoord = (n12-n21);

	\\ create vector of curve and coordinate
	returnvector = [curve, xcoord];

	return(returnvector);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetjinv(net)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms,
\\ returns the j-invariant of the associated curve.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
twonetjinv(net) = 

	return( elllong(twonettocurvechar2(net)[1])[13] );

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetdisc(net)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms,
\\ returns the discriminant of the associated curve.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
twonetdisc(net) = 

	return( elllong(twonettocurvechar2(net)[1])[12] );

}





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SECTION TWO:  RELATING NETS AND SEQUENCES                \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonettoeds(net,a,b)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms,
\\ returns the normalised elliptic divisibility sequence
\\ associated to the direction (a,b).  That is, it returns
\\ [1, W(2a,2b)*W(a,b)^(-4), W(3a,3b)*W(a,b)^(-9), 
\\ W(4a,4b)*W(a,b)^(-16)]
\\
\\ Returns '0' on degenerate nets.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonettoeds(net,a,b) =

	local( returnvec, twonetab );

	returnvec = vector(4);

	\\ Return zero if the net is degenerate.
	if( istwonetdegen(net),
		if(debug, print("Degenerate nets taste bad!"); );
		return(0);
	);

	\\ Get the value of the net W(a,b)
	twonetab = twonetget(net,a,b);

	\\ Set up the eds
	returnvec[1] = 1;
	returnvec[2] = twonetget(net,2*a,2*b)*twonetab^(-4);
	returnvec[3] = twonetget(net,3*a,3*b)*twonetab^(-9);
	returnvec[4] = twonetget(net,4*a,4*b)*twonetab^(-16);

	\\ It's possible that the resulting eds is
	\\ degenerate, so notify if debug=1
	if( returnvec[2] == 0 || returnvec[3] == 0,
		if(debug, print("The resulting elliptic
			 divisibility sequence is degenerate."); );
	);

	return(returnvec);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ edstotwonet(eds1,eds2)
\\
\\ Given two elliptic divisibility sequences associated to the same
\\ curve, returns a net (as a 4-vector) of the curve and both points.
\\ 
\\ If the two eds are not from the same curve, returns zero.
\\
\\ Elliptic divisibility sequences must be normalised (first term 1).
\\ If they are not normalised, they will be scaled.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

edstotwonet(eds1,eds2) = 

	local(curvepoint1, curvepoint2, returnnet, coordchange,
		fullcurve1, fullcurve2);

	\\ get the curves and points of the sequences
	curvepoint1 = [edstocurve(eds1),[0,0]];
	curvepoint2 = [edstocurve(eds2),[0,0]];

	\\ check that the curves are the same up to unihomothetic 
	\\ change of variables
	\\ (actually this checks up to u = plus/minus 1)
	fullcurve1 = edstocurvefull(eds1);
	fullcurve2 = edstocurvefull(eds2);
	if( fullcurve1[13] != fullcurve2[13] ||
		fullcurve1[12] != fullcurve2[12] ||
		fullcurve1[11] != fullcurve2[11] ||
		fullcurve1[10] != fullcurve2[10],
		if( debug, 
			print( "Elliptic divisibility sequences are
				 not from the same curve." ); 
		);
		return(0);
	);

	\\ get coordinate change required to go from second curve
	\\ to first
	\\ since u = plus/minus 1, this should always be possible
	coordchange = getellcoordchange(curvepoint2[1],
			 curvepoint1[1]);
	if( coordchange == 0,
		if(debug, 
			print("Coordchange failed.");
		);
		return(0);
	);

	\\ change the second curve to go to the first
	curvepoint2 = [ ellchangecurve(curvepoint2[1],coordchange),
		 ellchangepoint(curvepoint2[2],coordchange) ];

	\\ get the net associated to the curve from eds1, and the
	\\ point from eds1 and the changed point from eds2
	returnnet = curvetotwonet(curvepoint1[1], curvepoint1[2],
		 curvepoint2[2]);

	return(returnnet);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ getellcoordchange(curve1,curve2)
\\
\\ Given two elliptic curves of the same j-invariant, calculate the
\\ change of variables required to go from first to second, if
\\ possible.  The change of variables may lie over an extension field,
\\ in which case getellcoordchange may fail and you can try defining
\\ the same curves over the extension field and trying again.
\\ 
\\ This algorithm will return 0 if it fails.  It may fail for many
\\ reasons, the most common being that the curves are not isomorphic
\\ or the isomorphism lies over an extension.  It will also fail
\\ if the j-invariant is zero.
\\
\\ Requires curves to be isomorphic and to be of the same pari type.  
\\ Mixed t_INT and t_FRAC is ok.
\\
\\ If you're having trouble with a curve of seemingly mixed type like
\\ [1,x,x^2,...], you can make 1 of type t_POL by using 1+0*x instead
\\ for example.  Pari doesn't square-root polys well, though, so you'll
\\ have better luck making it t_SER.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

getellcoordchange(curve1,curve2) =

	local(u,r,s,t, long, short, numer, denom, numerroot, denomroot,
		 uroot, u2, u4, i, k);
	long = vector(2);
	short = vector(2);
	uroot = matrix(2,2);


	\\ initialise curves for extended data
	\\ recall that this works for singular curves
	long[1] = elllong(curve1);
	long[2] = elllong(curve2);


	\\ also store a convenient short version of curves
	short[1] = vector(5);
	short[2] = vector(5);
	for(i=1,5,
		short[1][i] = curve1[i];
		short[2][i] = curve2[i];
	);


	\\ if the j-invariants are different there's no hope anyway
	if (long[1][13] != long[2][13],
		if(debug, 
			print("different j-invariants in");
			print(" getellcoordchange"); 
		);
		return(0);
	);

	\\ if the j-invariant is zero, this is an annoying and
	\\ difficult case which is not yet implemented.
	if ( long[1][13] == 0,
		if(debug,
			print("The j-invariant is zero.");
			print("  This case not implemented.");
		);
		return(0);
	);


	\\ in the case that only a unihomothetic change of variables
	\\ is needed, this routine is easier (no issues of
	\\ type/roots)
	\\ this is the basic algorithm used below also
	if( long[1][13] == long[2][13] &&
		long[1][12] == long[2][12] &&
		long[1][11] == long[2][11] &&
		long[1][10] == long[2][10],

		for(i=1,2,
			u = (-1)^i;
			s = 1/2*(curve2[1]*u - curve1[1]);
			r = 1/3*(curve2[2]*u^2 - curve1[2] 
				+ s*curve1[1] + s^2);
			t = 1/2*(curve2[3]*u^3 - curve1[3] 
				- r*curve1[1]);
	
			if(debug, print("trying: ", [u,r,s,t] ); );
			if( ellchangecurve(short[1],[u,r,s,t]) 
				== short[2],
				return([u,r,s,t]);
			);
		);
		\\ if we reached this point, u was plus/minus 1, but
		\\ somehow neither of the changes of coords worked
		if(debug,
			print("coordchange failed for u=\pm 1");
		);
		return(0);

	);

	\\ If execution gets here, u was not just plus/minus 1


	\\ There are issues with type comparison, so if it's mixed
	\\ type, return error.
	\\ t_INT and t_FRAC mixed are ok

	\\ obtain first type
	typebase = type(short[1][1]);

	\\ consider integers as if they are fractions
	if( typebase == "t_INT", typebase = "t_FRAC" );

	\\ loop through other types, considering integers as
	\\ as fractions, and watch for mismatch with first
	for(i=1,5, for(k=1,2,
		typecheck = type(short[k][i]);
		if( typecheck == "t_INT", typecheck = "t_FRAC" );
		if( typecheck != typebase,
			if(debug, print("input curves must be
				 all same types"); );
			return(0);
		);
	); );


	\\ fourth power of u is the ration of the c4's
	u4 = long[1][10]/long[2][10];


	\\ make sure taking square root of u^4 is possible
	if( !issquare(u4) && (type(u4) == "t_INTMOD" 
		|| type(u4) == "t_POLMOD"),

		\\ in this case u^4 is not a square
		\\ and the isomorphism of the curves
		\\ is defined over an extension field
		\\ try again in an extension field
		if(debug, 
			print("fourth power of u isn't a square!!");
			print("u^4 = ", u4 );
			print("try an extension field.");
		);

		return(0);
	);


	\\ get a square root.  
	u2 = sqrt(u4);

	\\ if it's a rational and square,
	\\ make sure to get it as type rational
	if( (type(u4) == "t_FRAC" || type(u4) == "t_INT" )
		&& u4 > 0,

		\\ break it up as numerator and denominator and
		\\ take the integer roots of each
		numer = numerator(u4);
		denom = denominator(u4);
		numerroot = sqrtint(numer);
		denomroot = sqrtint(denom);

		\\ set u2 as the root, if this worked
		\\ if this didn't work, give a message
		\\ (didn't work means wasn't a rational
		\\ square, that's all)
		if( (numerroot/denomroot)^2 == u4,
			u2 = numerroot/denomroot;
		,
			if(debug,
				print("Square root as rational");
				print("	didn't work in");
				print(" getellcoordchange.");
			);
		);
	);


	\\ loop through both square roots of u4
	for(i = 1,2,

		\\ do u2 and negative u2 in turn
		u2 = (-1) * u2;
	
		\\ check if it's a square, and if it is,
		\\ put the two roots in uroot
		\\ if it is not, just put uroot=1 as placeholder
		if( !issquare(u2) && (type(u2) == "t_INTMOD" 
			|| type(u2) == "t_POLMOD"),

			\\ u2 is not a square
			\\ and we are working with a modulus
			if(debug, 
				print("Quadratic non-residue");
				print(" in getellcoordchange.");
				print("  May need extension");
				print(" field?");
			); 
			uroot[1,i]=1;
			uroot[2,i]=1;
		,

			\\ we are working in t_FRAC, t_COMPLEX etc.
			\\ so try to take a square root
			\\ depending on type this may produce
			\\ a pari error for weird types
			u = sqrt(u2);

			\\ if the type was rational (and positive)
			\\ try to do the root as a rational if
			\\ possible.

			\\ if this fails, it means the root is over
			\\ an extension of the rationals
			if( (type(u2) == "t_FRAC" 
				|| type(u2) == "t_INT" ) && u2 > 0,

				numer = numerator(u2);
				denom = denominator(u2);
				numerroot = sqrtint(numer);
				denomroot = sqrtint(denom);
				if( (numerroot/denomroot)^2 == u2,
					u = numerroot/denomroot;
				,
					if(debug,
					    print("Square root");
					    print(" as rational");
					    print(" didn't work");
					    print(" in getell");
					    print("coordchange.");
					);
				);
			);

			\\ store the roots in uroot
			uroot[1,i] = u;
			uroot[2,i] = -u;
		);

	);

	\\ at this point, if all has gone well, we've stored four
	\\ roots of u^4 in uroot and we can test them all
	\\ as possibilities.

	\\ if all has not gone well, some roots were missed
	\\ and it is possible the change of variables
	\\ requires working over an extension field.

	for(i=1,2,for(k=1,2,

		\\ select the root for testing
		u = uroot[i,k];

		\\ setup change of variables
		s = 1/2*(curve2[1]*u - curve1[1]);
		r = 1/3*(curve2[2]*u^2 - curve1[2] + s*curve1[1] 
			+ s^2);
		t = 1/2*(curve2[3]*u^3 - curve1[3] - r*curve1[1]);
	
		\\ report if debug is on
		if(debug, print("trying: ", [u,r,s,t] ); );

		\\ if one of them works return it
		if( ellchangecurve(short[1],[u,r,s,t]) == short[2],
			return([u,r,s,t]);
		);
	););

	\\ if not change of coordinates was found, report this
	\\ and return 0
	if(debug, print("No change of coordinates found."); );

	return(0);

}





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





\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ istwonetdegen(net)
\\
\\ Returns 1 if net is degenerate, otherwise 0.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

istwonetdegen(net) = 


	\\ check each of the possible degenerate cases
	if( net[3] == net[4],
		if(debug, print("The net has P-Q = 0"); );
		return(1);
	);
	if( net[2] == 0 && net[4] == 0,
		if(debug, print("The net has P=0"); );
		return(1);
	);
	if( net[1] == 0 && net[3] == 0,
		if(debug, print("The net has Q=0"); );
		return(1);
	);
	return(0);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ istwonetsing(net)
\\
\\ Returns 1 if net is singular, otherwise 0.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

istwonetsing(net) = 

	\\ check if net is singular by looking at discriminant
	\\ of associated curve
	if( elllong(twonettocurvechar2(net)[1])[12] == 0,
		return(1);
	,
		return(0);
	);

}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetbasischange(net,a,b,c,d)
\\
\\ Given a non-degenerate elliptic net 'net' associated to E,P,Q in the
\\ form of a vector of four terms, and integers a,b,c,d,
\\ returns a net associated to E and aP + bQ, cP + dQ
\\
\\ Does this directly via formulas.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetbasischange(nett,a,b,c,d) = 

	local(newnetter,nettab, nettcd, nettacbd, nettsing);

	\\ the new (post basis change) net will be stored here
	newnetter = vector(4);

	\\ store some useful values of the old net
	nettab = twonetget(nett,a,b);
	nettcd = twonetget(nett,c,d);
	nettacbd = twonetget(nett,a+c,b+d);
	nettsing = twonetget(nett,a+c,-b-d);

	\\ don't allow basis change that involves zero terms
	if( nettab == 0 || nettcd == 0 || nettacbd == 0 
		|| nettsing == 0,
		if(debug, 
			print("illegal basis change (or resulting
				 net is degenerate)");
		);
		return(0);
	);

	\\ compute terms of new net
	newnetter[1] = twonetget(nett,2*a,2*b)/nettab^4;
	newnetter[2] = twonetget(nett,2*c,2*d)/nettcd^4;
	newnetter[3] = twonetget(nett,2*a+c,2*b+d)/nettab^2
			/nettacbd^2*nettcd;
	newnetter[4] = twonetget(nett,a+2*c,b+2*d)*nettab
			/nettacbd^2/nettcd^2;

	\\ error trap
	if( istwonetdegen(newnetter),
		if(debug, 
			print("Watch out!  New net is degenerate.
			(This should have been caught earlier.)");
		);
		print("ERROR, PLEASE REPORT THIS BUG 98273243611
			 to stange@math.brown.edu");
	);

	return(newnetter);


}


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetbasischangeviacurve(net,a,b,c,d)
\\
\\ Given a non-degenerate elliptic net 'net' associated to E,P,Q in the
\\ form of a vector of four terms, and integers a,b,c,d,
\\ returns a net associated to E and aP + bQ, cP + dQ
\\
\\ Does this by translating to curve and then back to net.
\\
\\ This produces the same output as twonetbasischange, but is often
\\ faster since it doens't require computing all sorts of elements
\\ of elliptic nets.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetbasischangeviacurve(nett,a,b,c,d) = 

	local(newnetter, curve, point1, point2, curvepoint);

	\\ vector to store new (post basis change) net
	newnetter = vector(4);

	\\ get curve and point associated to the net
	curvepoint = twonettocurve(nett);
	curve = curvepoint[1];
	point1 = [curvepoint[2],0];
	point2 = [-curvepoint[2],0];

	\\ error trap:  points should be on the curve
	if( !ellisoncurve(curve,point1) || !ellisoncurve(curve,point2),
		if(debug, 
			print("points not on curve"); 
		);
		print("ERROR, PLEASE REPORT THIS BUG 944466112 to
			stange@math.brown.edu");
		return(0);
	);

	\\ make a net from the curve and the new basis points
	newnetter = curvetotwonet(curve, elladd(curve,
		 ellpow(curve, point1, a), ellpow(curve, point2, b) ),
		 elladd(curve, ellpow(curve, point1, c),
		 ellpow(curve, point2, d) ) );

	\\ error trap:  making the new net failed
	\\ maybe it was degenerate
	if( newnetter == 0, 
		if(debug, print("curvetotwonet failed."); );
		return(0); 
	);

	\\ error trap:  curvetotwonet should not produce a degenerate
	\\ net
	if( istwonetdegen(newnetter),
		if(debug, 
			print("Watch out!  New net is degenerate.
		  	(This should have been caught earlier.)"); 
		);
		print("ERROR, PLEASE REPORT THIS BUG 98276112 to
			stange@math.brown.edu");
	);

	return(newnetter);


}




\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ SECTION FOUR:  CALCULATING VALUES OF NETS                \\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\




\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetget(net,x,y)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms, and integers x and y
\\ returns the value of the net at index (x,y).
\\
\\ Uses a log n algorithm based on shipsey's thesis.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetget(nett,x,y) =

	local(X,k,newnett, ar, cr);

	\\ This function is recursive, so this keeps track
	\\ of the recursion if debug is on.
	\\ debug = 1 in general with twonetget will produce
	\\ waaaaaayyy too much data.
	if(debug, print("Calling twonetget on W(", x, "," ,y, ")"); );

	\\ Degenerate nets are not allowed.  It is possible (but
	\\ rarely seen?) that this algorithm will recursively call
	\\ itself on a degenerate net somewhere in the process.
	\\ I would appreciate reports on this.
	if( istwonetdegen(nett),
		if(debug, 
			print("Degenerate nets taste like belly
			 button lint!"); 
		);
		\\ can't return 0 since that's a value!
		return("failed");
	);


	\\ sometimes there's REALLY no work to do
	\\ just return one of these simple values near the origin
	if( [x,y] == [2,0], return( nett[1] ); );
	if( [x,y] == [2,1], return( nett[3] ); );
	if( [x,y] == [1,1], return( 1 ); );
	if( [x,y] == [1,0], return( 1 ); );
	if( [x,y] == [0,0], return( 0 ); );	
	if( [x,y] == [0,-2], return( -nett[2] ); );
	if( [x,y] == [-1,-2], return( -nett[4] ); );
	if( [x,y] == [-1,-1], return( -1 ); );
	if( [x,y] == [0,-1], return( -1 ); );


	\\ if the answer is really small, just do it directly
	\\ by calling twonetarray to build the small values
	\\ this is more efficient than the recursive algorithm
	\\ for small values
	if( abs(x) < 5 && abs(y) < 5, 
		return( twonetarray(nett,5)[6+x,6+y] ); 
	);


	\\ make sure x >= y
	if( y > x, 
		return( twonetget( 
			[nett[2],nett[1],nett[4],nett[3]] ,y,x) );
	);


	\\ make sure x >= 0
	if( x < 0, 
		return( -twonetget( nett, -x, -y ) ); 
	);


	\\ if the requested coordinate is of the form
	\\ (positive,negative), use a basis change to change
	\\ it to (pos,pos)
	if( y < 0, 

		\\ change basis to P, -Q to do the calculation
		\\ with two positive values

		newnett = twonetbasischange(nett,1,0,0,-1);

		return( twonetget(newnett,x,-y)*(-1)^(x*y) );

	);


	\\ if y is 0 or 1, we'll need to get the x-axis of
	\\ the net to start things off
	if( y == 0 || y == 1,

		X = twonetarray(nett,5);
		k = 5+1;

	);


	\\ if y=0, just need EDS associated to first point
	if( y == 0,

		return(	edsget([X[k+1,k],X[k+2,k],
				X[k+3,k],X[k+4,k]],x) );

	);

	\\ if y=1, just need to calculate out Shipsey block
	\\ along the x-axis
	\\ See my Tate Pairing via Elliptic Nets paper
	\\ for calculating Shipsey blocks.
	if( y == 1,

		return( shipsey_block([X[k-2,k],X[k-1,k],X[k,k],
			X[k+1,k],X[k+2,k],X[k+3,k],X[k+4,k],
			X[k+5,k];X[k,k+1],X[k+1,k+1],X[k+2,k+1],
			0,0,0,0,0], nett, x ) );


	);

	\\ at this point, y >= 2, and x >= y
	if(debug, print("got to y >=2 case"); );

	if( twonetget(nett,0,y) == 0, 

		\\ if W(0,y) = 0, translate by it to get a
		\\ simpler thing to return
		if(debug, 
			print("We have W(0,y) = 0 for y=", y); 
		);

		if( twonetget(nett,0,2) != 0 
			&& twonetget(nett,2,0) != 0, 
			\\ then no division by zero in following
			ar = (twonetget(nett,2,y)
				/twonetget(nett,2,0)
				/twonetget(nett,1,y));
			cr = (twonetget(nett,1,y+1)
				*twonetget(nett,0,2)
				*twonetget(nett,0,y+1)
				/twonetget(nett,0,y+2)/ar);
			return( twonetget(nett, x, 0)*ar^x*cr);	
	
		, \\ one of those is zero

			if( twonetget(nett,2,0) != 0, 
				\\ W(0,2) is zero, so W(1,2) and
				\\ W(2,2) are not, then y is even
				ar = twonetget(nett,2,2)
					/twonetget(nett,2,0)
					/twonetget(nett,1,2);
				cr = twonetget(nett,1,2)/ar;
				return( twonetget(nett,x,0)
					*ar^(y/2)*cr );

			, 
				\\ in this case W(2,0) = 0 
				\\ so W(3,0) and W(3,1) not zero
				br = twonetget(nett,3,1)
					/twonetget(nett,3,0);
				ar = -twonetget(nett,2,1)/br;
				cr = -ar;

				if( Mod(x,2) == 0, 
					\\ then W(x,y) = 0
					return( 0 );					
				,
					return( twonetget(nett,1,y)
						*ar^((x-1)/2)
						*br^(y*(x-1)/2)
						*cr^(((x-1)/2)^2) );
				);

			);
		);

	);

	if( twonetget(nett,1,y) == 0, 

		\\ if W(1,y) = 0, translate by it to get
		\\ a simpler thing to return
		if(debug, 
			print("We have W(1,y) = 0 for y=", y);
		);

		if( twonetget(nett,0,2) != 0 
			&& twonetget(nett,2,0) != 0, 
			\\ then no division by zero in
			\\ the following
			ar = (twonetget(nett,3,y)
				/twonetget(nett,2,0)
				/twonetget(nett,2,y));
			cr = (twonetget(nett,2,y+1)
				*twonetget(nett,0,2)
				*twonetget(nett,1,y+1)
				/twonetget(nett,1,y+2)/ar);

			return( twonetget(nett, x-1, 0)
				*ar^(x-1)*cr);
		
		, 
			\\ one of those is zero
			if( twonetget(nett,2,0) != 0,
				\\ W(0,2) is zero, so W(1,2)
				\\ and W(2,2) are not, y odd
				ar = twonetget(nett,2,2)
					/twonetget(nett,2,0)
					/twonetget(nett,1,2);
				cr = twonetget(nett,1,2)/ar;
				return( twonetget(nett,x,1)
					*ar^((y-1)/2)*cr );

			, 
				\\ in this case W(2,0) = 0 
				\\ so W(3,0) and W(3,1) not 0
				br = twonetget(nett,3,1)
					/twonetget(nett,3,0);
				ar = -twonetget(nett,2,1)/br;
				cr = -ar;

				if( Mod(x,2) == 0, 
					return( twonetget(nett,
						0,y)*ar^(x/2)
						*br^(y*x/2)
						*cr^((x/2)^2) );					
				,  \\ then W(x,y) = 0
					return( 0 );
				);

			);
		);

	);


	if( twonetget(nett,-1,y) == 0, 
		\\ if W(1,-y) = 0, translate by it to get a 
		\\ simpler thing to return
		if(debug, 
			print("We have W(-1,y) = 0 for y=", y);
		);

		if( twonetget(nett,0,2) != 0 
			&& twonetget(nett,2,0) != 0, 
			\\ then no division by zero in following

			ar = (twonetget(nett,1,y)
				/twonetget(nett,2,0)
				/twonetget(nett,0,y));
			cr = (twonetget(nett,0,y+1)
				*twonetget(nett,0,2)
				*twonetget(nett,-1,y+1)
				/twonetget(nett,-1,y+2)/ar);
			return( twonetget(nett, x+1, 0)
				*ar^(x+1)*cr);	
	
		, 
			\\ one of those is zero
			if( twonetget(nett,2,0) != 0,
				\\ W(0,2) is zero, so W(1,2) and 
				\\ W(2,2) are not, y odd
				ar = twonetget(nett,2,2)
					/twonetget(nett,2,0)
					/twonetget(nett,1,2);
				cr = twonetget(nett,1,2)/ar;
				return( twonetget(nett,x,1)
					*ar^((y-1)/2)*cr );

			, 
				\\ in this case W(2,0) = 0 
				\\ so W(3,0) and W(3,1) not zero
				br = twonetget(nett,3,1)
					/twonetget(nett,3,0);
				ar = -twonetget(nett,2,1)/br;
				cr = -ar;

				if( Mod(x,2) == 0, 
					return( twonetget(nett,
						0,y)*ar^(x/2)
						*br^(y*x/2)
						*cr^((x/2)^2) );					
				,  
					\\ then W(x,y) = 0
					return( 0 );
				);

			);
		);

	);


	\\ otherwise change basis to P, yQ and use Shipsey block

	newnett = twonetbasischange(nett,1,0,0,y);
	if( newnett == 0, 
		\\ we just had a failed basis change since 
		\\ W(1,y) or W(0,y) = 0 or W(1,-y) = 0
		if(debug, 
			print("failed basis change that
				 should not have happened"); 
		);
	);

	return( twonetget(newnett,x,1)*(twonetget(nett,1,y)^x)
		/twonetget(nett,0,y)^(x-1) );


}



{

shipsey_block( theblock, nett, val ) = 
\\ given: theblock, the block centred on 1, in the net 'nett'
\\ cannot accept nett[1] = 0 or net such that
\\ twonetget(nett, 2, -1) = 0
\\ returns: the central value of block centred on val in the
\\ net 'nett' i.e. W(val,1)

	local(initial_data);
	initial_data = vector(4);

	if( debug, print("Calling Shipsey Block"); );

	if( nett[1] == 0 || twonetget(nett,2,-1),
		if(debug, 
			print("Shipsey Block Request with
				 zero W(2,0) or W(2,-1)"); 
		);
	);

	if( nett[1] == 0, 
		\\ in case W(2,0) = 0
		initial_data[1] = "inf";
	,
		initial_data[1] = 1/nett[1];
	);
	initial_data[2] = 1;
	\\ the following division should never be zero 
	\\ in non-degenerage net
	initial_data[3] = 1/twonetget(nett,-1,1); 

	if( twonetget(nett,2,-1) == 0, 
		\\ in case W(2,-1) = 0, we cannot
		\\ call net_loop
		if(debug, 
			print("Avoiding a W(2,-1) issue"); 
		);
		return( -twonetget(nett,val+2,0)
				*twonetget(nett,-3,1)
				/(twonetget(nett,-1,1)
				*twonetget(nett,2,0))^(val+3) );
	,
		initial_data[4] = 1/twonetget(nett,2,-1);
	);

	return( net_loop( theblock, initial_data, val )[2,2] );

}	


{

double_or_add(V, initial_data, add) =
\\ double_or_add
\\ Given a block V centred at k, and the initial data relevant to 
\\ the elliptic net, returns either a block centred at 2k or 2k+1
\\ depending on whether variable "add" is 0 or 1 respectively.


	local(doubleV, m, i, j);


	\\ Create the output block
	doubleV = matrix(2,8);

	
	\\ initial_data contains the precomputed inverses
	inverse_20 = initial_data[1]; \\ inverse of W(2,0) 
	inverse_11 = initial_data[2]; \\ inverse of W(1,1)
	inverse_n1 = initial_data[3]; \\ inverse of W(-1,1)
	inverse_2n = initial_data[4]; \\ inverse of W(2,-1) 



	\\ Fill out first vector of output block

	for(j=-1,2,

		i = j;
		m = 4; \\ index to middle of block

		if( inverse_20 == "inf", 
			\\ in case it's a net where 2P = 0,
			\\ even terms are all 0
			doubleV[1,m+2*i-add] = 0;
		,
			doubleV[1,m + 2*i - add] = ((V[1,m+i])
				*(V[1,m+i+2])*(V[1,m+i-1])^2 
				- (V[1,m+i])*(V[1,m+i-2])
				*(V[1,m+i+1])^2)*inverse_20;
		);

		\\ when we hit j=-1, if add=1,
		\\ calculate W(2k+5,0) instead of W(2k-3,0)
		if( i == -1 && add == 1, i = 3; );

		doubleV[1,m + 2*i - 1 - add] = (V[1,m+i+1])
			*(V[1,m+i-1])^3 - (V[1,m+i-2])
			*(V[1,m+i])^3;
	);


	\\Fill out second vector of output block

	m2 = 2; \\ index to middle of second vector of block
	m1 = 4; \\ index to middle of first vector of block

	if( add == 0,
		doubleV[2,1] = ( V[2,m2+1]*V[2,m2-1]
			*V[1,m1-1]^2 - V[1,m1]*V[1,m1-2]
			*V[2,m2]^2 )*inverse_11;
	);

	doubleV[2,2-add] = ( V[2,m2-1]*V[2,m2+1]*V[1,m1]^2
		 - V[1,m1-1]*V[1,m1+1]*V[2,m2]^2);
	doubleV[2,3-add] = ( V[2,m2+1]*V[2,m2-1]*V[1,m1+1]^2 
		- V[1,m1]*V[1,m1+2]*V[2,m2]^2)*inverse_n1;

	if( add == 1,
		if( inverse_2n == "inf",
			if(debug, 
				print("Oh dear, W(2,-1) = 0
				 in shipsey"); 
			);
			doubleV[2,3] = "inf";
		,
			doubleV[2,3] = ( V[1,m1+1]*V[1,m1+3]
				*V[2,m2]^2 - V[2,m2-1]*V[2,m2+1]
				*V[1,m1+2]^2 )*inverse_2n;
		);
	);


	return(doubleV);

}



{

net_loop(V,initial_data, m) =
\\ net_loop
\\ Given a starting block V centred at 1, initial data relevant to 
\\ the elliptic net, and an integer m, returns the block centred at m

	local(currentV, m_size, add, i, j);


	\\ determine the number of steps in the double-and-add loop
	m_size = ceil(log(m+1)/log(2));


	\\ the variable storing the current block
	currentV = V;


	\\ ignore the first "1" in the binary expansion of m
	m = m - 2^(m_size-1);


	\\ step through the digits in the binary expansion of m
	for(j=1,m_size-1,
		i = m_size - j; \\kludgy version of "down to"

		\\ determine if this is a double step or a
		\\ double-and-add step
		\\ based on current digit of m; set "add" accordingly
		if( m - 2^(i-1) >= 0,
			add = 1;
			m = m - 2^(i-1);
		,
			add = 0;
		);


		\\ call the double or double-and-add function to
		\\ update the current block
		currentV = double_or_add(currentV, initial_data, add);


		\\ print information about the current step
		if( debug,
			print("The digit of $m$ for this step is "
				 add, ".");
			\\print("The resulting block from this
			\\	 step is ", currentV);
		);

	);


	return(currentV);

}



\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetgetslow(eds,a,b)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms, and integers a and b
\\ returns the value of the net at index (a,b).
\\
\\ Horribly inefficient; generates the whole array to large width.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetgetslow(nety,x,y) =

	local(k,width);

	if( istwonetdegen(nety),
		if(debug, 
			print("Degenerate nets taste like toe jam!"); 
		);
		return("failed");
	);

	width = 2*max(abs(x),abs(y))+8;
	k = width + 1;

	return(twonetarray(nety,width)[k+x,k+y]);

}



\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetarray(net,width)
\\
\\ Given a non-degenerate elliptic net 'net' in the
\\ form of a vector of four terms, and integer width > 3,
\\ returns a large matrix (2*width+1 times 2*width+1)
\\ whose entries represent the entries of the elliptic net
\\
\\ If the returned matrix is called X, then 
\\ W(x,y) is stored in entry X[width+1+x,width+1+y].
\\ In particular, the array contains all W(x,y) with |x|,|y| <= width.
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetarray(startnet,width) =

	local( i,j,k,l,w,X,m );

	if(debug, print("Calling twonetarray"); );

	if( istwonetdegen(startnet),
		if(debug, 
			print("Degenerate nets taste like earwax!"); 
		);
		return(0);
	);

	\\ check to be sure width is an integer
	if( type(width) != "t_INT",
		if(debug, print("Width not an integer!"););
		return(0);
	);

	\\ check to be sure width is not too small
	if( width < 4, 
		print("width too small, increasing to width=4"); 
		width = 4; 
	);

	\\ set up the matrix
	w = 2*width+1;
	X = matrix(w,w);

	\\ Useful variable 'k' so that W(x,y) = X[k+x,k+y]
	k = width+1;

	\\ Error-catching background filler for matrix 
	\\ (if in the end you see 'c' there's been an error,
	\\ as this should all be overwritten)
	for(l=1,w,for(m=1,w,
		X[l,m] = c;
	););


	\\ Initial values of the net.
	X[k+0,k+0] = 0;
	X[k+0,k+1] = 1;
	X[k+1,k+0] = 1;
	X[k+1,k+1] = 1;
	X[k+2,k+0] = startnet[1];
	X[k+2,k+1] = startnet[3];
	X[k+0,k+2] = startnet[2];
	X[k+1,k+2] = startnet[4];



	\\ Starter recurrences fill out the area near the origin.
	\\ No division by zero possible here if net is nondegenerate.

	X[k+1,k-1] = (-X[k+2,k+1]*X[k,k+1]^3 + X[k+1,k+2]*X[k+1,k]^3)
		/(X[k+1,k+1]^3);

	if(X[k+1,k-1] == 0, 
		if(debug,print("Oops, P+Q=0");); 
		return(0);
	);

	X[k+2,k-1] = (X[k+2,k]*X[k+1,k]^2*X[k+0,k+2] - X[k+2,k+1]
		*X[k+1,k-1]^2*X[k,k+1])/(X[k,k+1]*X[k+1,k+1]^2);

	X[k+1,k-2] = (-X[k+2,k-1]*X[k,k+1]^3 + X[k+1,k+1]*X[k+1,k-1]^3)
		/(X[k+1,k]^3);

	X[k+3,k] = ( -X[k+1,k+1]*X[k+1,k-1]*X[k+2,k]^2 + X[k+2,k+1]
		*X[k+2,k-1]*X[k+1,k]^2) / (X[k+1,k]*X[k,k+1]^2);

	X[k,k+3] = ( X[k+1,k+1]*X[k+1,k-1]*X[k,k+2]^2 - X[k+1,k+2]
		*X[k+1,k-2]*X[k,k+1]^2) / ( X[k,k+1]*X[k+1,k]^2);

	X[k+2,k+2] = ( - X[k+2,k]*X[k+1,k+1]*X[k+1,k+2]*X[k,k+1] 
		+ X[k,k+2]*X[k+2,k+1]*X[k+1,k]*X[k+1,k+1])/(X[k,k+1]
		*X[k+1,k]*X[k+1,k-1]);

	X[k+2,k-2] = ( X[k+2,k]*X[k,k+1]*X[k+1,k-2]*X[k+1,k-1] 
		- X[k+2,k-1]*X[k,k+2]*X[k+1,k-1]*X[k+1,k])
		/(-X[k,k+1]*X[k+1,k]*X[k+1,k+1]);

	X[k+3,k+1] = ( - X[k+1,k+1]^2*X[k+2,k-1]*X[k+2,k+1] 
		+ X[k+2,k]^2*X[k+1,k]*X[k+1,k+2])/(-X[k+1,k-1]
		*X[k,k+1]);

	X[k+1,k+3] = ( X[k+1,k+1]^2*X[k+1,k-2]*X[k+1,k+2] + X[k,k+2]
		*X[k,k+2]*X[k,k+1]*X[k+2,k+1])/(X[k+1,k-1]*X[k+1,k]^2);

	X[k+2,k+3] = ( X[k+1,k]*X[k+1,k+2]*X[k,k+2]*X[k+2,k+2] 
		- X[k+1,k+3]*X[k+1,k+1]*X[k+2,k+1]*X[k,k+1] ) 
		/ ( X[k,k+1]*X[k+1,k-1]*X[k+1,k+1] );

	X[k+3,k+2] = ( X[k,k+1]*X[k+2,k+1]*X[k+2,k]*X[k+2,k+2] 
		- X[k+3,k+1]*X[k+1,k+1]*X[k+1,k+2]*X[k+1,k] ) 
		/ ( -X[k+1,k]*X[k+1,k-1]*X[k+1,k+1] );

	X[k+3,k+3] = ( - X[k+3,k+2]*X[k+1,k+1]*X[k+1,k+2]*X[k+1,k+1] 
		+ X[k+2,k+2]*X[k,k+1]*X[k+2,k+2]*X[k+2,k+1])
		/(-X[k+1,k]^2*X[k+1,k+1]);

	X[k+3,k-1] = ( X[k+2,k]^3*X[k,k+2] + X[k+1,k-1]^3*X[k+3,k+1] ) 
		/ (X[k+1,k+1]^3);

	X[k+1,k-3] = -( X[k,k+2]^3*X[k+2,k] - X[k+1,k-1]^3*X[k+1,k+3] ) 
		/ (X[k+1,k+1]^3);

	if( X[k+2,k] != 0, 

		\\ in this case 2P != 0 on the curve

		X[k,k+4] = ( -X[k+1,k+3]*X[k+1,k+1]*X[k,k+2]*X[k+2,k-2] 
			-X[k+1,k-1]*X[k+1,k-3]*X[k,k+2]*X[k+2,k+2] )
			/ (-X[k+2,k]*X[k+1,k-1]*X[k+1,k+1] );

	,  
		\\ now 2P = 0.  If also 2P - Q = 0,
		\\ then Q = 0 -- a degen. net.  So 
		\\ division by X[k+2,k-1] ok.

		if( X[k+2,k-1] == 0, 
			if(debug, print("error, Q=0"); );
			return(0); 
		);

		X[k,k+4] = ( -X[k+2,k-2]*X[k+1,k+2]*X[k,k+3]*X[k+1,k+1] 
			- X[k+2,k+2]*X[k+1,k-2]*X[k,k+1]*X[k+1,k-3] ) 
			/ ( -X[k+1,k]*X[k+1,k+1]*X[k+2,k-1] );

	);

	if( X[k,k+2] != 0, 
		\\ in this case 2Q != 0

		X[k+4,k] = ( X[k+3,k+1]*X[k+1,k+1]*X[k+2,k]*X[k+2,k-2] 
			- X[k+1,k-1]*X[k+3,k-1]*X[k+2,k]*X[k+2,k+2] ) 
			/ (X[k,k+2]*X[k+1,k-1]*X[k+1,k+1] );

	,  
		\\ now 2Q = 0.  If also P - 2Q = 0, then
		\\ P = 0 -- a degen. net.  So division by X[k+1,k-2] ok.

		if( X[k+1,k-2] == 0, 
			if(debug, print("error, P=0"); );
			return(0); 
		);
		X[k+4,k] = ( X[k+2,k-2]*X[k+2,k+1]*X[k+3,k]*X[k+1,k+1] 
			- X[k+2,k+2]*X[k+2,k-1]*X[k+1,k]*X[k+3,k-1] ) 
			/ ( X[k,k+1]*X[k+1,k+1]*X[k+1,k-2] );

	);

	X[k+4,k-1] = ( -X[k+3,k]*X[k+1,k-2]*X[k+2,k]^2 + X[k+1,k-1]
		*X[k+3,k+1]*X[k+2,k-1]^2 ) / (X[k,k+1]*X[k+1,k+1]^2);
	X[k+1,k-4] = ( -X[k,k+3]*X[k+2,k-1]*X[k,k+2]^2 + X[k+1,k-1]
		*X[k+1,k+3]*X[k+1,k-2]^2 ) / (X[k+1,k]*X[k+1,k+1]^2);

	\\if(debug, print("done starter recurrences"); );


	\\ Fill out axes in positive direction with EDS recurrence
	\\ The terms W(0,i) and W(i,0) for i=1,4 are already done
	\\ by starter recurrences

	for(i=5,width,

		\\ y-axis

		if( X[k,k+i-4] != 0, 

			\\ safe to divide by W(0,i-4)
			X[k,k+i]=(X[k,k+i-1]*X[k,k+i-3]*X[k,k+2]^2
				-X[k,k+3]*X[k,k+1]*X[k,k+i-2]^2)
				/(X[k,k+i-4]*X[k,k+1]^2);
		, 
			\\ if W(0,i-4) = 0 then W(0,i-5) != 0, so safe
			\\ to divide by W(0,i-5)

			\\ further, if W(0,2) = 0 = W(0,i-4), then net
			\\ nondegen => i is even, so W(0,i)=0 too

			if(X[k,k+2] == 0,
				if(debug && Mod(i,2) != 0, 
					print("Degenerate net error
						W(0,odd)=W(0,2)=0"); 
				); 
				X[k,k+i] = 0; 

			,

				if(i==5, 
					if(debug, 
						print("W(0,1)=0 error"); 
					); 
				);

				X[k,k+i]=(X[k,k+i-1]*X[k,k+i-4]*X[k,k+2]
					*X[k,k+3]-X[k,k+4]*X[k,k+1]
					*X[k,k+i-2]*X[k,k+i-3])
					/(X[k,k+i-5]
					*X[k,k+1]*X[k,k+2]);
			);
		);

		\\ x-axis
		if( X[k+i-4,k] != 0, 
			\\ safe to divide by W(i-4,0)
			X[k+i,k]=(X[k+i-1,k]*X[k+i-3,k]*X[k+2,k]^2
				-X[k+3,k]*X[k+1,k]*X[k+i-2,k]^2)
				/(X[k+i-4,k]*X[k+1,k]^2);
		,  
			\\ if W(i-4,0) = 0 then W(i-5,0) != 0,
			\\ so safe to divide by W(i-5,0)

			\\ further, if W(2,0) = 0 = W(i-4,0),
			\\ then net nondegen => i is even, so W(i,0)=0
			\\ too
			if(X[k+2,k] == 0,
				if(debug && Mod(i,2) != 0, 
					print("Degenerate net error 
						W(odd,0)=W(2,0)=0"); 
				); 
				X[k+i,k] = 0; 

			,

				if(i==5, 
					if(debug, 
						print("W(1,0)=0 error"); 
					); 
				);
				X[k+i,k]=(X[k+i-1,k]*X[k+i-4,k]
					*X[k+2,k]*X[k+3,k]
					-X[k+4,k]*X[k+1,k]
					*X[k+i-2,k]*X[k+i-3,k])
					/(X[k+i-5,k]
					*X[k+1,k]*X[k+2,k]);
			);
		);

	);


	\\if(debug, print("done axes"); );

	\\Fill out the terms 0-4 of W(-1,*).
	for(j=0,4,
		X[k-1,k+j] = -X[k+1,k-j];
	);

	
	\\Using translated sequence recurrences, fill out positive
	\\ rows and columns (full first quadrant)

	for(i=1,4,
		X = dnetgen_helper_row(X,width,i,-1);
		X = dnetgen_helper_col(X,width,i,-1);
	);

	for(i=5,width,
		X = dnetgen_helper_row(X,width,i,0);
		X = dnetgen_helper_col(X,width,i,0);
	);

	\\if(debug, print("done translated recs positive"); );

	\\ Fill out columns in negative direction

	for(j=1,width,
	for(m=1,width,
	i = -m;	

		\\ error-check
		if( debug && X[k+j,k+i+4] == 0 && X[k+j,k+i+5] == 0, 
			print("Yikes degenerate: W(", j, ",", i+4,
				 ")=0=W(", j, ",", i+5, ")" ); 
		);

		\\ fill out column j downwards
		if( X[k+j,k+i+4] != 0, 
			\\ okay to divide by W(j,i+4)
			X[k+j,k+i]=(X[k+j,k+i+1]*X[k+j,k+i+3]
				*X[k,k+2]^2-X[k,k+3]*X[k,k+1]
				*X[k+j,k+i+2]^2)
				/(X[k+j,k+i+4]*X[k,k+1]^2);
		, 
			\\ W(j,i+4)=0 so okay to divide by W(j,i+5)

			\\ further, if W(0,2) = 0 = W(j,i+4), then net
			\\ nondegen => i is even, so W(j,i)=0 too

			if(X[k,k+2] == 0,
				X[k+j,k+i] = 0; 
			,
				X[k+j,k+i]=(X[k+j,k+i+1]*X[k+j,k+i+4]
					*X[k,k+2]*X[k,k+3]-X[k,k+4]
					*X[k,k+1]*X[k+j,k+i+2]
					*X[k+j,k+i+3])
					/(X[k+j,k+i+5]*X[k,k+1]
					*X[k,k+2]);
			);
		);

	);
	);


	\\Fill out the negatives.
	for(i=-width,-1,for(j=-width,width,
		X[k+i,k+j]=-X[k-i,k-j];
	););

	for(i=-width,-1,
		X[k,k+i]=-X[k,k-i];
	);

	if(debug, print("returning twonetarray"); );

	return(X);

}


{
dnetgen_helper_row(X, width, row, startblock)=
\\ internal function that fills out a row in the
\\ positive direction of a 2d array

	local(j,i);

	j=row;
	for(i=startblock+5,width,

		if( X[k+i-4,k+j] != 0, 
			\\ okay to divide by W(i-4,j)
			X[k+i,k+j]=(X[k+i-1,k+j]*X[k+i-3,k+j]
				*X[k+2,k]^2-X[k+3,k]*X[k+1,k]
				*X[k+i-2,k+j]^2)/(X[k+i-4,k+j]
				*X[k+1,k]^2);

		, 

			\\ W(i-4,j) = 0 so okay to divide by W(i-3,j)

			\\ further, if W(2,0) = 0 = W(i-4,j), then
			\\ net nondegen => i is even, so W(i,j)=0 too
			if(X[k+2,k] == 0,
				X[k+i,k+j] = 0; 
			,
				X[k+i,k+j]=(X[k+i-1,k+j]
					*X[k+i-4,k+j]*X[k+2,k]
					*X[k+3,k]-X[k+4,k]*X[k+1,k]
					*X[k+i-2,k+j]*X[k+i-3,k+j])
					/(X[k+i-5,k+j]*X[k+1,k]
					*X[k+2,k]);
			);
		);
	);

	return(X);

}


{
dnetgen_helper_col(X, width, col, startblock)=
\\ internal function that fills out a column in the positive
\\ direction of a 2d array

	local(j,i);

	j=col;
	for(i=startblock+5,width,
		if( X[k+j,k+i-4] != 0, 
			\\ okay to divide by W(j,i-4)
			X[k+j,k+i]=(X[k+j,k+i-1]*X[k+j,k+i-3]
				*X[k,k+2]^2-X[k,k+3]*X[k,k+1]
				*X[k+j,k+i-2]^2)/(X[k+j,k+i-4]
				*X[k,k+1]^2);

		, 

			\\ W(j,i-4) = 0 so okay to divide by W(j,i-3)

			\\ further, if W(0,2) = 0 = W(j,i-4), then
			\\ net nondegen => i is even, so W(j,i)=0 too
			if(X[k,k+2] == 0,
				X[k+j,k+i] = 0; 
			,
				X[k+j,k+i]=(X[k+j,k+i-1]*X[k+j,k+i-4]
					*X[k,k+2]*X[k,k+3]-X[k,k+4]
					*X[k,k+1]*X[k+j,k+i-2]
					*X[k+j,k+i-3])/(X[k+j,k+i-5]
					*X[k,k+1]*X[k,k+2]);
			);
		);
	);

	return(X);

}



\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ twonetarrayprettify(array,typeofpretty)
\\
\\ Given an array produced by twonetarray, prettifies it depending
\\ on the value of typeofpretty:
\\
\\ cartesian:  makes x and y increase right and up from center
\\ posquad: shows positive quadrant only, cartesian style
\\ latex: outputs it in latex format
\\ maple: outputs it as matrix for maple
\\ mathematica:  outputs it as matrix for mathematica
\\ 
\\ For best results, write latex'd, maple'd or mathematica'd matrices
\\ to file, where the escape characters don't show (copy and paste
\\ from screen causes problems).  Use write command.
\\
\\ In general, these apply to any square matrix, and one should apply
\\ them, in the order desired, one at a time (for example, to make a
\\ nice latex table, do cartesian first, then latex.
\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{

twonetarrayprettify(array,typeofpretty)=

	local(newmat, newsize);

	\\ all cases except posquad use newsize
	newsize = length(array);

	if( typeofpretty == "cartesian",
		newmat = array;
		array = mattranspose(array);
		for(i=1,newsize,for(j=1,newsize,
			newmat[i,j] = array[newsize+1-i,j];
		););
		return(newmat);

	);
	if( typeofpretty == "posquad",
		newsize = (length(array)-1)/2;
		newmat = matrix(newsize, newsize);
		for(i=1,newsize,for(j=1,newsize,
			newmat[i,j] = array[newsize+i,newsize+j]
		););
		newmat = twonetarrayprettify(newmat,"cartesian");
		return(newmat);
	);
	if( typeofpretty == "latex",
		newmat = "";
		newmat = concat(newmat,"\\begin{matrix}");
		for(i=1,newsize-1,

			newmat = concat(newmat, array[i,1]);
			for(j=2,newsize,
				newmat = concat(newmat," & ");
				newmat = concat(newmat,array[i,j]);
			);
			newmat = concat(newmat," \\\\");

		);

		newmat = concat(newmat, array[newsize,1]);
		for(j=2,newsize,
			newmat = concat(newmat," & ");
			newmat = concat(newmat,array[newsize,j]);
		);
		
		newmat = concat(newmat,"\\end{matrix}");

		return(newmat);

	);
	if( typeofpretty == "maple",
		newmat = "";
		newmat = concat(newmat,"Matrix(");
		newmat = concat(newmat,newsize);
		newmat = concat(newmat,",");
		newmat = concat(newmat,newsize);
		newmat = concat(newmat,",[[");
		for(i=1,newsize-1,

			newmat = concat(newmat, array[i,1]);
			for(j=2,newsize,
				newmat = concat(newmat,",");
				newmat = concat(newmat,array[i,j]);
			);
			newmat = concat(newmat,"],[");

		);

		newmat = concat(newmat, array[newsize,1]);
		for(j=2,newsize,
			newmat = concat(newmat,",");
			newmat = concat(newmat,array[newsize,j]);
		);
		
		newmat = concat(newmat,"]])");

		return(newmat);


	);
	if( typeofpretty == "mathematica",
		newmat = "";
		newmat = concat(newmat,"{{");
		for(i=1,newsize-1,

			newmat = concat(newmat, array[i,1]);
			for(j=2,newsize,
				newmat = concat(newmat,",");
				newmat = concat(newmat,array[i,j]);
			);
			newmat = concat(newmat,"},{");

		);

		newmat = concat(newmat, array[newsize,1]);
		for(j=2,newsize,
			newmat = concat(newmat,",");
			newmat = concat(newmat,array[newsize,j]);
		);
		
		newmat = concat(newmat,"}}");

		return(newmat);

	);

	if( debug,
		print("not a valid type of prettification");
	);
	return(0);

}