# Calibration from a torus picture
#
# This Maple script creates a random symmetric elliptic quartic curve Rd
# with fixed bitangents I*u+/-v=0
# and then uses Algorithm TorusCalibration
# to recover the equation of the elliptic absolute.
#
# run it by "read(calibtorus);"
# some wrappers for elimination functions
elimout := proc(polylist::list, vars::list)
description "Computes the elimination ideal of the ideal generated by polylist with respect to the variables vars.";
# Input: polylist, a list of polynomials;
# vars, a list of variables
# Output: a Groebner basis of the elimination ideal obtained eliminating
# the variables vars.
local allVars, minusVars, G, groebnerElim;
# computes the indeterminates in polylist
allVars := indets(polylist);
minusVars := [op(allVars minus {op(vars)})];
G := Groebner[Basis](polylist, lexdeg(vars, minusVars));
groebnerElim := remove(has, G, vars);
return groebnerElim;
end proc:
elimin := proc(polylist::list, vars::list)
description "Computes the elimination ideal of the ideal generated by polylist keeping only the variables vars.";
# Input: polylist, a list of polynomials;
# vars, a list of variables
# Output: a Groebner basis of the elimination ideal
# obtained eliminating all variables
# that appear in polylist, except for the ones in vars.
local allVars, minusVars, G, groebnerElim;
allVars := indets(polylist);
minusVars := [op(allVars minus {op(vars)})];
G := Groebner[Basis](polylist, lexdeg(minusVars, vars));
groebnerElim := remove(has, G, minusVars);
return groebnerElim;
end proc:
# Random integers and rational numbers
Randi := proc(n)
local arb;
# random integer
arb := rand(-n..n);
arb();
end proc:
Randr := proc(n)
# random rational number
local arb,r;
arb := rand(-n..n);
r := (arb()+n+1)/(arb()+n+1);
if rand() (-u:v:w)
Rd := u^4+(p0*v^2+p1*v*w+p2*w^2)*u^2+q0*v^4+q1*v^3*w+q2*v^2*w^2+q3*v*w^3+q4*w^4:
# force a singularity in the dual of the dual (a bitangent of Rd)
# we plug the line into the quartic and test if we get a perfect square
sbs := subs(u=1,factor(subs(v=I*u,Rd))):
FS := [op(subs(
b0=coeff(sbs,w,0),
b1=coeff(sbs,w,1),
b2=coeff(sbs,w,2),
b3=coeff(sbs,w,3),
b4=coeff(sbs,w,4),
SQr))]:
FF := [
op(map(T->coeff(expand(T),I,0),FS)),
op(map(T->coeff(expand(T),I,1),FS))
]:
F1 := elimout(FF,[]):
# (p0,..,q4) is in the zero set of F1 iff both lines v+I*u=0 and v-I*u=0 are bitangents
# Force a random singularity in Rd itself
# we may assume w=0 because one can show there is a projective transformation
# fixing symmetry and bitangent and bringing the w to 0
s1 := Randr(siz):
F2 := map(T->subs(u=s1,v=1,w=0,diff(Rd,T)),[u,v,w]):
# select a random example with these two constraints
# we observed that the variables p2,q1,q2 can be chosen generically
Fi := elimout([op(F1),op(F2),p2-Randr(siz),q1-Randr(siz),q2-Randr(siz)],[]):
Rd := subs(solve(Fi),Rd);
# create the torus image by dualizing the dual
R := elimout([x-diff(Rd,u),y-diff(Rd,v),z-diff(Rd,w),Rd],[u,v,w])[1]:
# R has 8 nodes, 6 outside the symmetry line
# one complex pair is (1:+/-I:0), used in the algorithm
# we test if there are more complex pairs
# the y-coordinates of the nodes to check appear as a factor of the discriminant
# this factor has degree 4 and multiplicity 2
# we test if it has real zeroes
DisF := factors(discrim(R,z))[2]:
onodes := subs(x=1,select(T->(T[2]=2 and degree(T[1])=4),DisF)[1][1]):
if nops([fsolve(onodes)]) = 4 then "no other complex nodes"; else "more complex nodes"; fi;
# compute symmetry of Rd
ssu := {u=a11*u+a12*v+a13*w,v=a21*u+a22*v+a23*w,w=a31*u+a32*v+a33*w}:
Eq := [coeffs(expand(subs(ssu,Rd)-Rd),[u,v,w])]:
Eq := [op(Eq),coeffs(expand(charpoly(matrix([[a11,a12,a13],[a21,a22,a23],[a31,a32,a33]]),t)-(t+1)^2*(t-1)),t)]:
SYS := elimout([op(elimout([op(Eq)],[])),a33*nz-1],[nz]):
syd := subs(solve(elimout([op(SYS)],[])),ssu):
if nops([syd])=1 then "only one symmetry"; else "more symmeties"; fi;
# only one symmetry!
# candidate for dual elliptic absolute - symmetric
Cand := a1*u^2+a2*v^2+a3*v*w+a4*w^2:
# condition imposed by bitangent of Rd (=tangent for Cand)
E1 := discrim(subs(u=1,subs(v=I*u,Cand)),w):
# tangential intersection test
# since we are taking resultant with respect to u,
# the resultant is always a perfect square of a quartic
sbs := subs(w=1,v=t,sqrt(resultant(Rd,Cand,u),symbolic)):
# we test if this quartic is the perfect square of a quadric
FS := [op(subs(
b0=coeff(sbs,t,0),
b1=coeff(sbs,t,1),
b2=coeff(sbs,t,2),
b3=coeff(sbs,t,3),
b4=coeff(sbs,t,4),
SQr)),a1*nz-1]:
FF := elimout(FS,[nz]):
# (a1,a2,a3,a4) is in the zero set of FF iff Cand and and Rd have double intersection
# inequalities: Cand should not pass through the singular points of Rd (it suffices to test one)
# and Cand should have nonzero hessian
ineq := subs(u=s1,v=1,w=0,Cand)*det(hessian(Cand,[u,v,w])):
# solve system of equations and inequalities
G := elimout([op(FF),E1,ineq*nz-1,a1-1],[nz]):
if map(degree,G) = [1,1,1,1] then "system has a single solution"; else "system has more solutions"; fi;
Ad := subs(solve(G),Cand):
# dualize the dual elliptic absolute
A := elimout([x-diff(Ad,u),y-diff(Ad,v),z-diff(Ad,w),Ad],[u,v,w])[1];
if coeff(A,x,2)*det(hessian(A,[x,y,z])) > 0 then "positive definite"; else "not positive definite"; fi;