// In the SINGULAR procedures below, one inputs coordinates of a point
// a in P^6 (in integers), and the final procedure will determine whether
// the divisor B(a) is smooth on E_1^(2) where, as in the paper, E_1 is 
// the elliptic curve which is the quotient of y^2=(x-1)(x^2+1) by the 
// subgroup <(1,0)> of order 2.  

int chr=0;

proc singV(int a0, int a1, int a2, int a3, int a4, int a5, int a6)
// Given point a in P^6 (specified by integer coordinates), outputs the
// number of singular points on \tilde B(a) in the fiber over \tilde 0
// for which such that Z1 is nonzero.
{
    	ring RV=chr,(t,V),dp;
		poly A=2;
    	poly psi0=A;
    	poly psi1=V6;
    	poly psi2=t*(V*(V4+A));
    	poly psi3=V*(V4-A);
    	poly psi4=V2;
    	poly psi5=V4;
    	poly psi6=t*V3;
    	poly comb=a0*psi0+a1*psi1+a2*psi2+a3*psi3+a4*psi4+a5*psi5+a6*psi6;
    	ideal sing=comb,t,diff(comb,t),diff(comb,V);
    	sing=std(sing);
    	"   Found", vdim(sing), "singularities with multiplicity";
}

proc singU(int a0, int a1, int a2, int a3, int a4, int a5, int a6)
// Given point a in P^6 (specified by integer coordinates), outputs the
// number of singular points on \tilde B(a) in the fiber over \tilde 0
// for which such that Z0' is nonzero.
{
    ring RU=chr,(t,U),dp;
	poly A=2;
	poly psi0=A*U6;
	poly psi1=1;
	poly psi2=t*(U*(1+A*U4));
	poly psi3=U*(1-A*U4);
	poly psi4=U4;
	poly psi5=U2;
	poly psi6=t*U3;
    	poly comb=a0*psi0+a1*psi1+a2*psi2+a3*psi3+a4*psi4+a5*psi5+a6*psi6;
    	ideal sing=comb,t,diff(comb,t),diff(comb,U);
    	sing=std(sing);
    	"   Found", vdim(sing), "singularities with multiplicity";
}

proc singAway(int a0, int a1, int a2, int a3, int a4, int a5, int a6)
// Given point a in P^6 (specified by integer coordinates), outputs the
// number of singular points of \tilde B(a) away from the fibers over \tilde 0
// and C, such that Z1 is nonzero.  (Note by symmetry under the involution that
// there are no singularities away from these two fibers if and only if there are
// no such singularities with Z1 nonzero.)
{
	ring RAway=chr,(x,y,t,Z),dp;
	poly curve=y2-x3+x2-x+1;
	poly alpha=1;
	poly A=2;
	poly inv=(x-alpha)*t-1;
	poly psi0=(x-alpha)*Z^6+A*t;
	poly psi1=(x-alpha)^2*Z^6+A^2*t2;
	poly psi2=(x-alpha)*Z^5+A*t*Z;
	poly psi3=y*Z^5-A*y*t2*Z;
	poly psi4=Z^4+Z^2;
	poly psi5=(x-alpha)*Z^4+A*t*Z^2;
	poly psi6=Z^3;
	poly comb=a0*psi0+a1*psi1+a2*psi2+a3*psi3+a4*psi4+a5*psi5+a6*psi6;
	ideal I=curve,comb,inv;
	ideal J=std(I);
	"   Equations give algebraic set of dimension ", dim(J);
	ideal sing=curve,comb,inv,minor(jacob(I),3);
	sing=std(sing);
	"   Singular locus has dimension", dim(sing);
	"   Found", vdim(sing), "singularities with multiplicity";
}

proc checkAll(int a0, int a1, int a2, int a3, int a4, int a5, int a6)
// Given point a in P^6 (specified by integer coordinates), performs all
// three procedures above.
{
	"";
	"Checking away from 0 and C...";
	singAway(a0,a1,a2,a3,a4,a5,a6);
	"";
	"Checking over 0 with Z1=1...";
	singV(a0,a1,a2,a3,a4,a5,a6);
	"";
	"Checking over 0 with Z0'=1...";
	singU(a0,a1,a2,a3,a4,a5,a6);
	"";
}