"
Oliver Labs
RICAM (Linz, Austria)

www.OliverLabs.net
mail@OliverLabs.net

www.ricam.oeaw.ac.at
Oliver.Labs@oeaw.ac.at


This Singular code computes the objects used in the article
'A Septic with 99 Real Nodes'. 
So, it might be helpful if one wants to work
with these objects, e.g. with the conics through 6 singularities. 

The script also computes ideals describing the singularities. 
This can e.g. be used to check the reality assertion of the theorem.

Finally, we check using these explicit coordinates of the singularities
that the surface actually has 99 nodes and no other singularities.

This script does not follow the path of the paper, but it assumes from the beginning
that we already know the condition on alpha.


This version: November 2005, RICAM (Linz, Austria).
Original version: August 2004, University of Mainz (Mainz, Germany).
-----

";

LIB "all.lib";

LIB "olitools.lib";

option(redSB);

// To check how much time the computation takes:
int tStart = rtimer;

// If you want to see the code also in the output 
// then set echo = 1.
echo = 0;

// If you do not want to check the following things
// (these computations take some time), 
// then set nolongchecks = 1,
// else set nolongchecks = 0.
// - there are no other singularities (takes some seconds)
// - checks in 3-space over Q(alpha)
int nolongchecks = 0;

// If you do not have our Singular library defsnodal.lib,
// you have to set have_defsnodal = 0,
// else you can set have_defsnodal = 1.
int have_defsnodal = 0;
if(have_defsnodal == 1) {
    LIB "defsnodal.lib";
}

proc multdim(ideal I)
"
PURPOSE: Show some data related to I, e.g., its dim() and mult().
RETURN: nothing
ASSUME: I is a standard basis.
"
{
    degree(I);
    std_primdecGTZ(I,2); 
    "the first coordinate:";
    eliminate(I,prod_of_vars(I)/(ideal_of_vars(I)[1]));
}


///////////////////////////////
// the ring in which our 99-nodal septic is defined:
//
ring r = (0, alpha), (x,y,w,z), dp;
minpoly = 7*alpha^3 + 7*alpha + 1;

///////
// The parameters t, a(1), ..., a(5) 
// as indicated in the article:
//
number t = 1/(alpha^2+1);

// number a(4) = (alpha*(1+alpha^2)-1)*(1+alpha^2)^2;
// number a(3) = -( -256*a(4)^2*t^12-256*a(4)^2*t^11+192*a(4)^2*t^10+192*a(4)^2*t^9-512*a(4)*t^10-32*a(4)^2*t^8-512*a(4)*t^9-32*a(4)^2*t^7+384*a(4)*t^8-4*a(4)^2*t^6+384*a(4)*t^7-256*t^8-64*a(4)*t^6-256*t^7-64*a(4)*t^5-64*t^6-8*a(4)*t^4+192*t^5+416*t^4-32*t^3-4*a(4)*t-228*t^2+a(3)-3*a(4)+28 - a(3));
// number a(2) = -( a(3)*t-3*a(4)*t+a(2)-2*a(3)+3*a(4) - a(2));
// number a(1) = -( 8*a(4)*t^3+a(1)-(-( a(3)*t-3*a(4)*t+a(2)-2*a(3)+3*a(4) - a(2)))+a(3)-a(4)+8*t - a(1));
// number p_comp_a5_denom = leadcoef(t^2 * 
// (8*a(4)*t^5+a(3)*t^3-3*a(4)*t^3-a(3)*t^2+2*a(4)*t^2+8*t^3-1)^2);	
// number p_comp_a5_num = leadcoef(448*a(4)^2*t^12+96*a(3)*a(4)*t^10-288*a(4)^2*t^10-80*a(3)*a(4)*t^9+176*a(4)^2*t^9+896*a(4)*t^10+5*a(3)^2*t^8-30*a(3)*a(4)*t^8+45*a(4)^2*t^8-8*a(3)^2*t^7+42*a(3)*a(4)*t^7-54*a(4)^2*t^7+96*a(3)*t^8-288*a(4)*t^8+3*a(3)^2*t^6-14*a(3)*a(4)*t^6+16*a(4)^2*t^6-80*a(3)*t^7+112*a(4)*t^7+448*t^8-6*a(3)*t^5+18*a(4)*t^5+448*t^6+4*a(3)*t^4-10*a(4)*t^4-64*t^5-560*t^4+169*t^2-7);
// number a(5) = p_comp_a5_num;
// a(5) = a(5) / p_comp_a5_denom;

// the a(i) as given in the paper on page 4:
number a(1) = alpha^7+7*alpha^5-alpha^4+7*alpha^3-2*alpha^2-7*alpha-1;
number a(2) = (alpha^2+1)*(3*alpha^5+14*alpha^3-3*alpha^2+7*alpha-3);
number a(3) = (alpha^2+1)^2*(3*alpha^3+7*alpha-3);
number a(4) = (alpha*(1+alpha^2)-1)*(1+alpha^2)^2;
number a(5) = (1+alpha^2)/(alpha^2);

"";"The parameters:";

"t =", t;
"a_1 =", a(1);
"a_2 =", a(2);
"a_3 =", a(3);
"a_4 =", a(4);
"a_5 =", a(5);


///////
// The basic objects used for the construction in projective 3-space:
//
"";"";
"---------------------------------------";
"";

// the doubled cubic C:
poly dCxyzw = (x^2+y^2)*(z+w)+a(1)*z^3+a(2)*z^2*w+a(3)*z*w^2+a(4)*w^3;

// the 7 planes 
poly P = 21*x^4*y^2*z+x^7-64*z^7+21*x^2*y^4*z-7*y^6*x-21*x^5*y^2+35*x^3*y^4+112*x^2*z^5+112*z^5*y^2+7*z*y^6+7*x^6*z-56*z^3*x^4-56*z^3*y^4-112*z^3*x^2*y^2; 

// the 99-nodal septic:
poly S = P - (z+a(5)*w)*dCxyzw^2;
poly Saff = subst(S,w,1);


////////////////
// the objects in the projective y=0 plane:
//

// the plane curve S_y:
poly Sy = substitute(S,y,0);

// the doubled cubic dC:
poly C = substitute(dCxyzw,y,0);

// the three doubled lines L1L2L3:
poly L1L2L3 = -x^3-4*x^2*z+4*x*z^2+8*z^3;

// the line Sy1:
poly Sy1 = (z+w)-1/(1+alpha^2)*x;


////////
// switch to the affine chart w=1:
//

Sy = subst(Sy,w,1);
Sy1 = subst(Sy1,w,1);
L1L2L3 = subst(L1L2L3,w,1);
C = subst(C,w,1);


//////////
// Let us switch to the ring in which our affine plane curves live:
//
ring rxz = (0, alpha), (x,z), dp;
minpoly = 7*alpha^3 + 7*alpha + 1;

poly Sy = imap(r, Sy);
poly Sy1 = imap(r, Sy1);
poly L1L2L3 = imap(r, L1L2L3);
poly C = imap(r, C);

/////////////////////////////////
// compute all the objects used in the article:

// the sextic Sy6:
Sy1 = simplify(Sy1,1);
poly Sy6 = Sy / Sy1;

// all 9 generic singularities:
ideal I_gen_Sy = C, L1L2L3;
I_gen_Sy = std(I_gen_Sy);

// the 6 generic singularities of the sextic Sy6:
ideal I_gen_Sy6 = quotient(I_gen_Sy, Sy1);
I_gen_Sy6 = std(I_gen_Sy6);

// the 3 generic singularities on Sy1:
ideal I_gen_Sy1 = I_gen_Sy, Sy1;
I_gen_Sy1 = std(I_gen_Sy1);
"";"the 3 generic singularities on Sy1:";I_gen_Sy1;
multdim(I_gen_Sy1);

// the conic C_0 through these singularities:
poly C_0 = degreepart(I_gen_Sy6, 0, 2)[1];
"";"the conic C_0:", C_0;

// the two points of C_0 on the x=0 axes:
poly C_0_x = substitute(C_0, x, 0);
"";"the two points of C_0 on the x=0 axes:";
multdim(std(C_0_x));

// the two possible values of beta:
matrix cc = coef(C_0_x,z);
poly beta_discr = cc[2,2]^2 - 4*cc[2,1]*cc[2,3];
poly extquad = x^2-beta_discr;
list feq = factorize(extquad);
list sols_beta;
sols_beta = list(leadcoef(-(feq[1][2]-x)), 
		 leadcoef(-(feq[1][3]-x)));
"";"the two possible values of beta:";sols_beta;
string sBeta = string(sols_beta[1]);

// compute the two possible equations for C_1:
list nums = list((-cc[2,2] + sols_beta[1]),
		 (-cc[2,2] + sols_beta[2]));
poly denom = 2*cc[2,1];
list C1s = list((nums[1]+denom)*x^2 -4*denom*z*(z+1) + (nums[1]+denom)*4*z,
		   (nums[2]+denom)*x^2 -4*denom*z*(z+1) + (nums[2]+denom)*4*z);
C1s[1] = C1s[1] / leadcoef(C1s[1]);
C1s[2] = C1s[2] / leadcoef(C1s[2]);
//"C1s_a:", C1s[1];
//"C1s_b:", C1s[2];

// compute the three general singularities on C1s:
ideal I_gen_C1(1) = C1s[1], L1L2L3;
I_gen_C1(1) = quotient(I_gen_C1(1), z^3);
I_gen_C1(1) = std(I_gen_C1(1));
//dim(I_gen_C1(1)), mult(I_gen_C1(1));
ideal I_gen_C1(2) = C1s[2], L1L2L3;
I_gen_C1(2) = quotient(I_gen_C1(2), z^3);
I_gen_C1(2) = std(I_gen_C1(2));
//dim(I_gen_C1(2)), mult(I_gen_C1(2));
ideal I_gen_C1_check(1) = I_gen_C1(1), jacob(Sy6);
I_gen_C1_check(1) = std(I_gen_C1_check(1));
//dim(I_gen_C1_check(1)), mult(I_gen_C1_check(1));
ideal I_gen_C1_check(2) = I_gen_C1(2), jacob(Sy6);
I_gen_C1_check(2) = std(I_gen_C1_check(2));
//dim(I_gen_C1_check(2)), mult(I_gen_C1_check(2));
// check, which of the 2 values of beta is the one 
// on which the three special singularities are:
ideal I_gen_C1;
number sol_beta;
ideal C_1;
if(mult(I_gen_C1_check(1)) == 3) {
    I_gen_C1 = I_gen_C1(1);
    sol_beta = sols_beta[1];
    C_1 = C1s[1];
} else {
    I_gen_C1 = I_gen_C1(2);
    sol_beta = sols_beta[2];
    C_1 = C1s[2];
}

// C_1:
"";"the conic C_1:";C_1;

// the 3 generic singularities on C_1:
"";"the 3 generic singularities on C_1:";I_gen_C1;

// the 3 exceptional singularities of Sy6:
ideal I_E = C_1, jacob(Sy6);
I_E = quotient(I_E, I_gen_C1);
I_E = std(I_E);
"";"the 3 special singularities of Sy6";I_E;
multdim(I_E);

// the three generic singularities on C_2:
ideal I_gen_C2 = quotient(I_gen_Sy, Sy1);
I_gen_C2 = quotient(I_gen_C2, I_gen_C1);
I_gen_C2 = std(I_gen_C2);
"";"the three generic singularities on C_2:";I_gen_C2;
multdim(I_gen_C2);

// C_2:
ideal I_C_2 = intersect(I_gen_C2, I_E);
poly C_2 = degreepart(I_C_2, 0, 2)[1];
"";"the conic C_2:"; C_2;

// the two other singularities on Sy1, called O_{12}:
ideal O_12 = Sy6, Sy1;
O_12 = quotient(O_12, x); // quotient out P_x
O_12 = quotient(O_12, L1L2L3); // quotient out the generic singularities
O_12 = std(O_12);
"";"the two other singularities, O_12, on Sy1:";O_12;
multdim(O_12);

// the singularity on the rotation axes x=y=0, P_x:
ideal P_x = x, z+1;
"";"the singularity on the rotation axes x=y=0, P_x:"; P_x;


///////////////////////////
// Now we can check that the points are really 
// ordinary double points of the plane septic:
//
"";"";
"---------------------------------------";
"";

// The 15 singularities of the plane septic:
"";"";"the 15 singularities of the plane septic:";
"";"P_x:"; multdim(std(P_x));
"";"O_12:"; multdim(O_12);
"";"I_gen_Sy1:"; multdim(I_gen_Sy1);
"";"I_gen_C1:";multdim(I_gen_C1);
"";"I_gen_C2:";multdim(I_gen_C2);
"";"I_E:"; multdim(I_E);


proc num_A1_check(ideal S, ideal singpoints, int should_be_num)
{
    // Check that the ideal is reduced
    // and that these points are A_1 singularities of Sy.
    // We use the hessian of Sy and count their number.
    //
    // Return 1, if all this is fulfilled.
    //
    // Return 0, if not.
    // In this case, display the computed numbers.

    ideal numpts = singpoints;
    numpts = std(numpts);
    int d_numpts = dim(numpts); 
    int m_numpts = mult(numpts);
    numpts = radical(numpts);
    numpts = std(numpts);
    int d_numpts_rad = dim(numpts); 
    int m_numpts_rad = mult(numpts);

    ideal jS = jacob(S);
    ideal hessS = det(jacob(jS));
    ideal num_check;
    ideal A1_check;

    // check the total milnor number:
    num_check = jS, S, singpoints;
    num_check = std(num_check);
    // the following dimension should be 0:
    int d_num_check = dim(num_check);
    // the following number should be should_be_num:
    int m_num_check = mult(num_check);
    A1_check = num_check, hessS;
    A1_check = std(A1_check);
    // the following dimension should be -1:
    int d_A1_check = dim(A1_check);
    // the following number should be 0:
    int m_A1_check = mult(A1_check);

    if(d_num_check == 0 
       && m_num_check == should_be_num
       && d_A1_check == -1
       && m_A1_check == 0
       && d_numpts == 0
       && m_numpts == should_be_num
       && d_numpts_rad == 0
       && m_numpts_rad == should_be_num) {
	"";return("Okay!");
    } else {
	"";"---> Problem:";
	"";"The ideal:";
	"dim =",d_numpts;
	"mult =",m_numpts;
	"dim =",d_numpts_rad;
	"mult =",m_numpts_rad;
	"";"S:";
	"dim =",d_num_check;
	"mult =",m_num_check;
	"dim hess =",d_A1_check;
	"mult hess =",m_A1_check;
	"";return("PROBLEM!");
    }
}

"";"Check that each of the ideals is reduced";
"and check also that these are all ordinary double points of Sy:";
num_A1_check(Sy, P_x, 1);
num_A1_check(Sy, O_12, 2);
num_A1_check(Sy, I_gen_Sy1, 3);
num_A1_check(Sy, I_gen_C1, 3);
num_A1_check(Sy, I_gen_C2, 3);
num_A1_check(Sy, I_E, 3);

"";"";
"---------------------------------------";
"";

"";"Check that these are really 15 distinct points:";
ideal allsings = intersect(P_x, O_12, I_gen_Sy1, I_gen_C1, I_gen_C2, I_E);

ring rproj = (0,alpha),(x,w,z), dp; minpoly = 7*alpha^3+7*alpha+1;

ideal allsings = imap(rxz, allsings);
allsings = std(homog(allsings,w));
mult(allsings);

"";"Check that there are no other singularities on the proj. curve Sy:";
poly Sy = imap(rxz,Sy);
Sy = homog(Sy,w);
ideal noother_check = jacob(Sy);
noother_check = std(noother_check);
degree(noother_check);
setring rxz;

/////////////
// To verify the hypotheses of Lemma 1, 
// we do exactly the same checks for the 15
// points when viewing them as points in the (x,y,z)-space
// (with y-coordinate =0) by adding the equation y=0 to the ideals:

ring rxyz = (0, alpha), (x,y,z), dp;
minpoly = 7*alpha^3 + 7*alpha + 1;

poly Saff = imap(r, Saff);
ideal P_x = imap(rxz, P_x), y;
ideal O_12 = imap(rxz, O_12), y;
ideal I_gen_Sy1 = imap(rxz, I_gen_Sy1), y;
ideal I_gen_C1 = imap(rxz, I_gen_C1), y;
ideal I_gen_C2 = imap(rxz, I_gen_C2), y;
ideal I_E = imap(rxz, I_E), y;

/////
// First, check that (1:i:0:0) does not lie on the surface:
setring r;
poly Si = substitute(S, x,1,z,0,w,0);
ideal yi = y^2+1;
yi = std(yi);
"";"check that (1:i:0:0) does not lie on the surface,";
"i.e. the following is not zero:", reduce(Si,yi);
setring rxyz;

/////
// Check the condition on the nodes:
"";"Check that each of the ideals is reduced";
"and check also that these are all ordinary double points of S:";
num_A1_check(Saff, P_x, 1);
num_A1_check(Saff, O_12, 2);
num_A1_check(Saff, I_gen_Sy1, 3);
num_A1_check(Saff, I_gen_C1, 3);
num_A1_check(Saff, I_gen_C2, 3);
num_A1_check(Saff, I_E, 3);

"";

"Time needed for the whole computation:";
timeSinceHuman(tStart);
"";
$;