\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ 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 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);
}