proc isTrue(int c) { if(c==0) { return("FALSE"); } else { return("TRUE"); } }

ring r = (17,alpha), (x,y,z,w), dp; minpoly = 7*alpha^3 + 7*alpha + 1;
number a(1) = -12/7*alpha^2 - 384/49*alpha - 8/7;
number a(2) = -32/7*alpha^2  + 24/49*alpha - 4;
number a(3) =    -4*alpha^2  + 24/49*alpha - 4;
number a(4) =  -8/7*alpha^2   + 8/49*alpha - 8/7;
number a(5) =    49*alpha^2      - 7*alpha + 50;
poly P = x*(x^6-3*7*x^4*y^2+5*7*x^2*y^4-7*y^6) 
    + 7*z*((x^2+y^2)^3-2^3*z^2*(x^2+y^2)^2+2^4*z^4*(x^2+y^2)) - 2^6*z^7;
poly U = (z+a(5)*w)*(a(1)*z^3+a(2)*z^2*w+a(3)*z*w^2+a(4)*w^3+(z+w)*(x^2+y^2))^2;
poly S = P-U;

"1. Check that the point (1:i:0:0) is not on S:";
poly Si = subst(subst(subst(S, x,1), z,0), w,0); 
ideal yi = y^2+1; yi = std(yi); 
isTrue(reduce(Si, yi) != 0), ": S(1:i:0:0) =", reduce(Si, yi); "";
"2. Check that there is exactly one node on x=y=0 axes:";
ideal jSxy = x, y, jacob(S); jSxy = std(jSxy);
isTrue(mult(jSxy)==1), ": milnor =", mult(jSxy), ", dim =", (dim(jSxy)-1); "";
"3. Check that milnor(S_y) = 15 (takes some minutes):";
ideal jSy = y, jacob(S); jSy = std(jSy);
isTrue(mult(jSy)==15), ": milnor =", mult(jSy), ", dim =", (dim(jSy)-1); "";
"4. Switch to the affine chart w=1."; "";
S = subst(S,w,1);
ring raff = (17,alpha), (x,y,z), dp; minpoly = 7*alpha^3 + 7*alpha + 1;
poly S = imap(r,S);
"5. Check that all the 15 sing. are in the affine chart w=1 (takes some minutes):";
ideal jSyaff = y, S, jacob(S); jSyaff = std(jSyaff);
isTrue(mult(jSyaff)==15), ": milnor =", mult(jSyaff), ", dim =", dim(jSyaff); "";
"6. Check that all the 15 singularities are nodes (takes approx. an hour):";
poly hessian = det(jacob(jacob(S)));
ideal nonnodes = y, S, jacob(S), hessian; nonnodes = std(nonnodes);
isTrue(mult(nonnodes)==0), ": milnor =", mult(nonnodes), ", dim =", dim(nonnodes); "";
$;