LIB "enumpoints.lib";
// The set of nodes is a 0-dimensional variety over a number field.
ring R2 = (0,q),(x,y,z,w),dp;
minpoly = q16-q12+q8-q4+1;
poly s = 2q13-q9-q7+q5-q3-q;
ideal Togliatti = 64*(x-w)*(x^4 -4*x^3*w -10*x^2*y^2 -4*x^2*w^2 +16*x*w^3 -20*x*y^2*w+5*y^4 +16*w^4 -20*y^2*w^2) -5*s*(2*z -s*w)*(4*(x^2+y^2-z^2) +(1+3*(5-s^2))*w^2)^2;
matrix Jac = jacob(Togliatti);
ideal I2 = Togliatti+Jac;
list L=projPoints(std(I2));
L[1];
==> [1]:
==>    1
==> [2]:
==>    0
==> [3]:
==>    (1/2q13-1/4q9-1/4q7+1/4q5-1/4q3-1/4q)
==> [4]:
==>    (-q12+q8+3/2)
size(L);
==> 31