#==============================================================
# NUMERICAL COMPUTATION OF DETERMINANTAL REPRESENTATIONS
# OF PLANE CURVES USING THE HELTON-VINNIKOV FORMULA
# Maple worksheet illustrating the material in
# "COMPUTING LINEAR MATRIX REPRESENTATIONS OF
# HELTON-VINNIKOV CURVES"
# by Daniel Plaumann, Bernd Sturmfels, and Cynthia Vinzant
# arXiv ...
# =============================================================
restart:QygtSSV3aXRoRzYiNiNJJnBsb3RzRzYkJSpwcm90ZWN0ZWRHSShfc3lzbGliR0YlISIiLUYkNiNJKmFsZ2N1cnZlc0dGKEYrLUYkNiNJLkxpbmVhckFsZ2VicmFHRihGKw==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYrLUkjbWlHRiQ2JlEld2l0aEYnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1ElYm9sZEYnLyUrZm9udHdlaWdodEdGNC1JKG1mZW5jZWRHRiQ2JS1GIzYpLUYsNiZRKmFsZ2N1cnZlc0YnRi9GMkY1LyUnZmFtaWx5R1EwVGltZXN+TmV3flJvbWFuRicvJStmb3JlZ3JvdW5kR1EoWzAsMCwwXUYnLyUrYmFja2dyb3VuZEdRLlsyNTUsMjU1LDI1NV1GJy8lK2V4ZWN1dGFibGVHRjFGMkY1RjJGNS1JI21vR0YkNi5RIjpGJ0YyRjUvJSZmZW5jZUdGMS8lKnNlcGFyYXRvckdGMS8lKXN0cmV0Y2h5R0YxLyUqc3ltbWV0cmljR0YxLyUobGFyZ2VvcEdGMS8lLm1vdmFibGVsaW1pdHNHRjEvJSdhY2NlbnRHRjEvJSdsc3BhY2VHUSwwLjI3Nzc3NzhlbUYnLyUncnNwYWNlR0ZobkY/RkJGRUZIRjJGNQ==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYrLUkjbWlHRiQ2JlEld2l0aEYnLyUnaXRhbGljR1EmZmFsc2VGJy8lLG1hdGh2YXJpYW50R1ElYm9sZEYnLyUrZm9udHdlaWdodEdGNC1JKG1mZW5jZWRHRiQ2JS1GIzYpLUYsNiZRLkxpbmVhckFsZ2VicmFGJ0YvRjJGNS8lJ2ZhbWlseUdRMFRpbWVzfk5ld35Sb21hbkYnLyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GJy8lK2JhY2tncm91bmRHUS5bMjU1LDI1NSwyNTVdRicvJStleGVjdXRhYmxlR0YxRjJGNUYyRjUtSSNtb0dGJDYuUSI6RidGMkY1LyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGaG5GP0ZCRkVGSEYyRjU=precision:=20: Digits:=precision:#==============================================================
# INPUT: A polynomial f in two variables x and y
# OUTPUT: A determinantal representation of the form
# f = det(Id + x*D + y*R)
# with D diagonal and R symmetric.
# There exist such D and R with real entries if and only if
# f is a Helton-Vinnikov curve. The simplest method for
# finding such a curve as input, is to input it in the
# desired form using the matrices D0 and R0.
# This worksheet is still experimental and rather unstable.
# Polynomials with small integer coefficients seem to work best.
#===============================================================
D0:=DiagonalMatrix([-2,-1,1,2]):R0:=Matrix([[1,2,0,0],[-1,0,0,-1],[0,-1,0,0],[1,0,-1,2]]):
#We make R0 symmetric.
R0:=R0^++R0;
d:=Dimension(R0):f:=Determinant(IdentityMatrix(d)+x*R0+y*D0);
#Another sample curve:
#f:=(1/2)*x^4+(1/2)*y^4+(1/2)*x^2*y^2-(3/2)*x^2-(3/2)*y^2+1;QyQtSTBwbG90X3JlYWxfY3VydmVHNiI2JUkiZkdGJUkieEdGJUkieUdGJSIiIg==d:=degree(f):g:=genus(f,x,y):#Check whether the projective closure is smooth (=has the right genus)
evalb(g=(d-1)*(d-2)/2);#INTERSECTION POINTS with the line x=0:
Q0:=puiseux(f,x=0,y,0,t):Q:=[]: for rt in Q0 do
Q:=[op(Q),allvalues(rt)]
end do:#Get the points numerically, too:u:=[]:for rt in Q do u:=[op(u),evalf(subs(t=0,rhs(rt[2])))]end do:u;#Should have d points here:
evalb(nops(Q)=d);#BASE POINT for the Abel-Jacobi map. Choose to your liking.
P0List:=[]: P0puiseux:=puiseux(f,x=0,y,0,t):
for L in P0puiseux do
P0List:=[op(P0List),allvalues(L)]:
end do:P0:=P0List[1];#PERIOD MATRIX:
AB:=periodmatrix(f,x,y):A:=SubMatrix(AB,[1..g],[1..g]): B:=SubMatrix(AB,[1..g],[g+1..2*g]):Omega:=A^(-1).B:
#DIFFERENTIALS
w:=Vector(differentials(f,x,y)):w:=A^(-1).w:
#THETA-CHARACTERISTICS: You may pick four vectors with entries in 0,1/2. They must of length g.
# 4(a,b) must be EVEN and 4(e0,e1) ODD
a:=<0,0,0>: b:=<0,0,0>: e0:=<1/2,0,0>: e1:=<1/2,0,0>:
# We make sure the vectors are of the right length and change them in case they are not.
if (Dimension(a)!=g) or (Dimension(b)!=g) or (Dimension(e0)!=g) or (Dimension(e1)!=g) then
aold:=a: bold:=b: e0old:=e0: e1old:=e1:
a:=Vector(g): b:=Vector(g): e0:=Vector(g): e1:=Vector(g):
for i from 1 to Dimension(aold) do
a[i]:=aold[i]:
end do:
for i from 1 to Dimension(bold) do
b[i]:=bold[i]:
end do:
for i from 1 to Dimension(e0old) do
e0[i]:=e0old[i]:
end do:
for i from 1 to Dimension(e1old) do
e1[i]:=e0old[i]:
end do:
end if:#CHECK THAT (a,b) is an even theta-characteristic:
if not (2*a.(2*b) mod 2 =0) then
print("WARNING: (a,b) is not an even theta-characteristic! Things will go wrong!"):
end if: if not (2*e0.(2*e1) mod 2 = 1) then
print("WARNING: (e0,e1) is not an odd theta-characteristic! Things will go wrong!"):
end if:
#NOW DEFINE FUNCTIONS THAT GO INTO THE FORMULA
#THETA FUNCTION WITH CHARACTERISTIC:
ThetaChar := proc(P,a,b)
local z;
z:=Vector(P);
evalf(RiemannTheta(z+Omega.a+b,Omega) * exp(Pi*I*(2*a^+.(z+b)+a^+.Omega.a)));
end proc:
#DENOMINATOR OF THE PRIME FORM
wtheta := proc(Pnum)
local i,E,form;
form:=0;
for i from 1 to g do
E:=UnitVector(i,g);
form:=form+evalf(2*e0[i]*Pi*I*exp(Pi*I*(e0^+.Omega.e0+2*e0.e1))*RiemannTheta(Omega^+.e0+e1,Omega)+exp(Pi*I*(e0^+.Omega.e0+2*e0^+.e1))*RiemannTheta(Omega.e0+e1,Omega,[E]))*w[i];
end do;
form:=subs(x=Pnum[1],y=Pnum[2],form);
return form;
end proc:
#ABEL-JACOBI MAP
phi := proc(P)
Vector(AbelMap(f,x,y,P,P0,t,precision)):
end proc:
#DIAGONAL MATRIX D:
DiagEntries:=[]:
for uroot in u do
DiagEntries:=[op(DiagEntries),-1/uroot]
end do:
Diag:=DiagonalMatrix(DiagEntries);#COMPUTATION OF THE MATRIX R:
#DIAGONAL ENTRIES
R := Matrix(d):
for i from 1 to d do
R[i,i]:=subs(x=0,y=u[i],-diff(f,x)/diff(f,y))/u[i]:
end do:
#IMAGES OF INTERSECTION POINTS UNDER THE AJ MAP
xlist:=[]:
for i from 1 to d do
xlist:=[op(xlist),phi(Q[i])]:
end do:
#EVALUATION OF THE FORMULA FOR THE OFF-DIAGONAL ENTRIES
for i from 1 to d do
for j from 1 to i-1 do
x1:=xlist[i]; x2:=xlist[j];
R[i,j]:=(1/u[i]-1/u[j])/ThetaChar(ZeroVector(g),a,b)*ThetaChar(x2-x1,a,b)/ThetaChar(x2-x1,e0,e1)*sqrt(-subs(dx=1,wtheta([0,u[i]]))*u[i])*sqrt(-subs(dx=1,wtheta([0,u[j]]))*u[j]):
R[j,i]:=R[i,j]:
end do:
end do:
#HERE IS THE RESULT
f0:=Determinant(IdentityMatrix(d)+x*R+y*Diag);#AND THE DIFFERENCE WITH THE ORIGINAL POLYNOMIAL
f-f0;#IF YOU GET WRONG RESULTS, TRY RUNNING THE FOLLOWING LINES:
#FOR UNKNOWN REASONS, THE OFF-DIAGONAL ENTRIES SOMETIMES HAVE TO BE MULTIPLIED BY A COMMON FACTOR.
R1:=Matrix(d):
for i from 1 to d do
for j from 1 to i-1 do
R1[i,j]:=s*R[i,j]: R1[j,i]:=s*R[j,i]:
end do:
end do:
f2:=Determinant(IdentityMatrix(d)+x*R1+y*Diag):sfactors:=solve(subs(y=0,coeff(f,x^4))=subs(y=0,coeff(f2,x^4)),s);Rlist:=[]:
for s0 in [sfactors] do
R2:=subs(s=s0,R1):
Rlist:=[op(Rlist),R2]:
end do:#PRINT THE SCALINGS AND THE DIFFERENCES WITH THE ORIGINAL POLYNOMIAL f:
for i from 1 to nops(Rlist) do
f3:=Determinant(IdentityMatrix(d)+x*Rlist[i]+y*Diag):
print(sfactors[i]);
print(f-f3);
end do: