
 //////////////////////////////////////////////////////
 //                                                  //
 //   kleinmatrix.m                                  //
 //                                                  //
 //   Computations accompanying the article          //
 //   "Genus one curves defined by Pfaffians"        //
 //                                                  //
 //   Author: T.A. Fisher                            // 
 //   Version: September 2004                        //
 //                                                  //
 //////////////////////////////////////////////////////

function PolyCoeffs(poly) 
   R:=Parent(poly);
   d:=TotalDegree(poly);
   mons:=MonomialsOfDegree(R,d);
   return [MonomialCoefficient(poly,mon):mon in mons];
end function;

function MatDelete(mat,d)
   n:=Nrows(mat);
   for i in [d..n-1] do
      mat:=SwapRows(mat,i,i+1);
      mat:=SwapColumns(mat,i,i+1);
   end for;
   return Submatrix(mat,1,1,n-1,n-1);
end function;

function Pfaffian(mat)
   n:=Nrows(mat);
   if n eq 0 then return 1; end if;
   if n eq 1 then return 0; end if;
   if n ge 2 then 
     mat1:=MatDelete(mat,1);
     return &+[(-1)^i*mat[1,i]*Pfaffian(MatDelete(mat1,i-1)):i in [2..n]];
   end if;
end function;

function LeadingSubmatrices(r,mat)
   n:=Nrows(mat);
   ss:=[SetToSequence(s): s in Subsets({1..n},r)];
   return [Matrix(r,r,[mat[s[i],s[j]]: i,j in [1..r]]): s in ss];
end function;

function PfIdeal(r,mat)
   R:=Parent(mat[1,1]);
   submats:=LeadingSubmatrices(r,mat);
   return ideal<R|[Pfaffian(m): m in submats]>;
end function;

function PfSubmax(mat)
   n:=Nrows(mat);
   return [(-1)^i*Pfaffian(MatDelete(mat,i)):i in [1..n]];
end function;

beta:=func<r,n|Binomial(n-r,r)+Binomial(n-r-1,r-1)>;

function IsHomogeneousBasis(ff)
   if #ff eq 0 then return true; end if;
   R:=Parent(ff[1]);
   r:=TotalDegree(ff[1]);
   for f in ff do
     if not IsHomogeneous(f) then return false; end if; 
     if TotalDegree(f) ne r then return false; end if;  
   end for;
   V:=VectorSpace(BaseRing(R),#MonomialsOfDegree(R,r));
   FF:=[V!PolyCoeffs(f): f in ff];
   return IsIndependent(FF);
end function;

function NextSecantVariety1(ff);
   R:=Parent(ff[1]);
   K:=BaseRing(R);
   n:=Rank(R);
   r:=TotalDegree(ff[1]);
   mons1:=MonomialsOfDegree(R,r);
   mons2:=MonomialsOfDegree(R,r+1);
   V1:=VectorSpace(K,#mons1);
   V2:=VectorSpace(K,#mons2);  
   FF:=[V1!PolyCoeffs(f):f in ff];
   W,map1:=quo<V1|FF>;
   map1:=Hom(V1,W)!map1;
   U:=V2;
   for i in [1..n] do
     map2:=Hom(V2,V1)![Degree(mon,R.i) eq 0 select V1!0 else
       Degree(mon,R.i)*V1.Position(mons1,mon/R.i): mon in mons2];
     U:=U meet Kernel(map2*map1);
   end for;
   GG:=[ElementToSequence(u): u in Basis(U)];
   return [&+[G[i]*mons2[i]: i in [1..#mons2]]: G in GG];
end function;

function NextSecantVariety2(ff);
   R:=Parent(ff[1]);
   K:=BaseRing(R);
   n:=Rank(R);
   r:=TotalDegree(ff[1]);
   mons0:=MonomialsOfDegree(R,r-1);
   mons1:=MonomialsOfDegree(R,r);
   mons2:=MonomialsOfDegree(R,r+1);
   V0:=VectorSpace(K,#mons0);
   V1:=VectorSpace(K,#mons1);
   V2:=VectorSpace(K,#mons2);
   VV:=VectorSpace(K,#ff*n);
   FF:=[V1!PolyCoeffs(f):f in ff];
   W,map1:=quo<V1|FF>;
   map1:=Hom(V1,W)!map1;
   dd:=[[V0![MonomialCoefficient(Derivative(f,R.j),mon): mon in mons0]
     : j in [1..n]]: f in ff];
   homs:=[Hom(V0,V1)![V1.Position(mons1,mon*R.i): mon in mons0]: i in [1..n]];
   U:=VV;
   targetdim:=(n-2*r)*beta(r,n);
   for k in [1..n] do
     map2:=Hom(VV,V1)![homs[i](dd[j][k]): i in [1..n],j in [1..#ff]];
     U:=U meet Kernel(map2*map1);
     dimU:=Dimension(U);
     if dimU eq targetdim then break; end if;
   end for;
   GG:=[ElementToSequence(u): u in Basis(U)];
   polys:=[R.i*f: i in [1..n],f in ff];
   solns:=[&+[G[i]*polys[i]: i in [1..#ff*n]]: G in GG];
   U:=sub<V2|[V2![MonomialCoefficient(q,mon):mon in mons2]: q in solns]>;
   GG:=Basis(U); 
   return [&+[G[i]*mons2[i]: i in [1..#mons2]]: G in GG]; 
end function;

function PolynomialContent(poly)
   if poly eq 0 then return 0; end if;
   cc:=Coefficients(poly);
   d:=LCM([Denominator(c):c in cc]);
   return GCD([IntegerRing()!(c*d): c in cc])/d;
end function;

function MatrixContent(mat)
   cc:=Eltseq(mat);
   d:=LCM([Denominator(c):c in cc]);
   return GCD([IntegerRing()!(c*d): c in cc])/d;
end function;

function PolyMatrixContent(mat)
   m:=Nrows(mat);n:=Ncols(mat);
   cc:=[PolynomialContent(mat[i,j]): i in [1..m], j in [1..n]];
   d:=LCM([Denominator(c):c in cc]);
   return GCD([IntegerRing()!(c*d): c in cc])/d;
end function;

function BEmatrix(ff,t)
   R:=Parent(ff[1]);
   K:=BaseRing(R);
   n:=#ff;
   zero:=ZeroMatrix(R,n,n);
   mons0:=MonomialsOfDegree(R,t);
   mons1:=MonomialsOfDegree(R,t+1);
   mat:=ZeroMatrix(BaseRing(R),n^2,#mons1);
   for i in [1..n] do
     for j in [1..#mons0] do
       pos:=Position(mons1,mons0[j]*R.i);
       for k in [1..#ff] do
         mat[#ff*(i-1)+k,pos]:=MonomialCoefficient(ff[k],mons0[j]);
       end for;
     end for;  
   end for;
   ns:=NullspaceMatrix(mat);
   if Nrows(ns) ne n then return zero; end if;
   Psi:=Matrix(R,n,n,[&+[ns[i,n*(j-1)+k]*R.j:j in [1..n]]:i,k in [1..n]]);
   mat:=ZeroMatrix(BaseRing(R),n^2,n^3);  
   for i in [1..n] do
     for j in [1..n] do
       for ii in [[i,j],[j,i]] do
       for r in [1..n] do
       for s in [1..n] do
         mat[n*(ii[1]-1)+r,n^2*(i-1)+n*(j-1)+s]+:=
           MonomialCoefficient(Psi[r,ii[2]],R.s);
       end for; end for; end for; 
     end for;  
     if Rank(mat) ge n^2-1 then break; end if; 
   end for;  
   ns:=NullspaceMatrix(mat);
   if Nrows(ns) ne 1 then return zero; end if;
   B:=Matrix(n,n,[ns[1,n*(i-1)+j]:i,j in [1..n]]);
   Phi:=B*Psi;
   return Phi;
end function;

function kmat(ff,n)
   assert IsOdd(n) and (n ge 5);
   r:=IntegerRing()!((n-3)/2);
   gg:=ff;     
   for rho in [1..r] do
     assert IsHomogeneousBasis(gg);
     assert TotalDegree(gg[1]) eq rho+1;
     assert #gg eq beta(rho+1,n);  
     printf "I(Sec^%o C) generated by %o %o-ics\n",rho,#gg,rho+1;
     if rho eq r then break; end if;
     if n*beta(rho+1,n) gt Binomial(n+rho+1,rho+2) then 
       gg:=NextSecantVariety1(gg);
     else
       gg:=NextSecantVariety2(gg);
     end if;
   end for;
   Phi:=BEmatrix(gg,r+1);
   pp:=PfSubmax(Phi);
   c:=pp[1]/gg[1];
   if pp ne [c*g:g in gg] then c:=0; end if;
   return Phi,c;
end function;

function KleinMatrix(C)
   I:=Ideal(C);
   n:=Degree(C);
   ff:=MinimalBasis(I);
   return kmat(ff,n);
end function;

function IsKleinMatrix(C,Phi);
   R:=CoordinateRing(AmbientSpace(C));
   n:=Nrows(Phi);
   Phi1:=Matrix(R,n,n,[Phi[i,j]: i,j in [1..n]]);
   flag0:=(Transpose(Phi1) eq -Phi1);
   flag1:=(PfIdeal(4,Phi1) subset Ideal(C)); 
   flag2:=IsHomogeneousBasis(PfSubmax(Phi1));
   return flag0 and flag1 and flag2;
end function;

function MatrixCongruence(Phi1,Phi2)
   P1:=PfSubmax(Phi1);
   P2:=PfSubmax(Phi2);
   n:=Nrows(Phi1);
   r:=IntegerRing()!((n-3)/2);
   R:=Parent(Phi1[1,1]);
   V:=VectorSpace(BaseRing(R),Binomial(n+r,r+1));
   W:=VectorSpace(BaseRing(R),n);
   phi1:=hom<W->V|[V!PolyCoeffs(p): p in P1]>;
   phi2:=hom<W->V|[V!PolyCoeffs(p): p in P2]>;
   phi:=phi2*phi1^(-1);
   B:=Matrix(n,n,[ElementToSequence(phi(W.i)): i in [1..n]]);
   Psi:=Transpose(B)*Phi2*B;
   for i in [1..n] do
   for j in [1..n] do
     if Psi[i,j] ne 0 then
       c:=Phi1[i,j]/Psi[i,j]; 
       break; 
     end if;
   end for; end for; 
   assert Phi1 eq c*Psi;
   return B,c;
end function;

////////////////////////////////////////////////////////////////////

 ///////////////////////////////////////////////////
 //                                               //
 //   Routines for testing and giving examples    //
 //                                               //
 ///////////////////////////////////////////////////

// Same as magma's "AntisymmetricMatrix", except we specify 
// the entries _above_ the diagonal.

function AlternatingMatrix(n,ll)
   mat:=ZeroMatrix(Parent(ll[1]),n,n);
   ctr:=0;
   for i in [1..n] do
   for j in [i+1..n] do
     ctr+:=1;
     mat[i,j]:=ll[ctr];
   end for;
   end for;
   return mat-Transpose(mat);
end function;

//  Section 8. EXAMPLES

//  8.1. Elliptic curves

function Phi5(E,PP)
   x1,x2,x3,x4,x5:=Explode([PP.i:i in [1..5]]);
   a1,a2,a3,a4,a6:=Explode(aInvariants(E));
   ll:=[a2*x4+a6*x1,x5,x4-a1*x3+a4*x1,-x3,-x4,-x3-a3*x1,x2,-x2,0,-x1];
   return AlternatingMatrix(5,ll);
end function;

function Phi7(E,PP)
   x1,x2,x3,x4,x5,x6,x7:=Explode([PP.i:i in [1..7]]);
   a1,a2,a3,a4,a6:=Explode(aInvariants(E));
   ll:=[x7,0,-p,x3,x5,0,x6,-q,-x5,0,x4,x5,-x2,-x4,0,r,s,-x3,x3,x1,x2]
     where p is x6-a1*x5+a2*x4-a3*x3+a4*x2+a6*x1
     where q is a3*x5-a4*x4-a6*x2
     where r is x4-a1*x3+a2*x2
     where s is a3*x3-a4*x2-a6*x1;  
   return AlternatingMatrix(7,ll);
end function;

function Weierstrass5(E)
   k:=BaseField(E);
   P2<x,y,z>:=AmbientSpace(E);
   P4<x1,x2,x3,x4,x5>:=ProjectiveSpace(k,4);
   phi:=map<E->P4|[z^2,x*z,y*z,x^2,x*y]>;
   C:=Image(phi);
   Phi1,c:=KleinMatrix(C);
   assert IsKleinMatrix(C,Phi1);
   Phi2:=Phi5(E,P4);
   assert IsKleinMatrix(C,Phi2);
   B,c:=MatrixCongruence(Phi1,Phi2);
   return Phi1,Phi2,B,c;
end function;

function Weierstrass7(E)
   k:=BaseField(E);
   P2<x,y,z>:=AmbientSpace(E);
   P6<x1,x2,x3,x4,x5,x6,x7>:=ProjectiveSpace(k,6);
   phi:=map<E->P6|[z^3,x*z^2,y*z^2,x^2*z,x*y*z,x^3,x^2*y]>;
   C:=Image(phi);
   Phi1,c:=KleinMatrix(C);
   assert IsKleinMatrix(C,Phi1);
   Phi2:=Phi7(E,P6);
   assert IsKleinMatrix(C,Phi2);
   B,c:=MatrixCongruence(Phi1,Phi2);
   return Phi1,Phi2,B,c;
end function;

//  8.2. Plane cubics

function DecomposeCubic(U)
   R<x,y,z>:=Parent(U);
   mons:=MonomialsOfDegree(R,2);
   U0:=U;
   aa:=[R!0: i in [1..3]];
   for i in [1..3] do
     aa[i]:=&+[MonomialCoefficient(U0,mon*R.i)*mon: mon in mons];
     U0:=U0-aa[i]*R.i;
   end for;
   assert U0 eq 0;
   assert U eq &+[aa[i]*R.i: i in [1..3]];
   return aa;
end function;

function Phi9(C3,PP)
   U:=Equation(C3);
   R3<x,y,z>:=Parent(U);
   aa:=DecomposeCubic(U);
   M:=AlternatingMatrix(4,[z,-y,aa[1],x,aa[2],aa[3]]);
   P:=Matrix(R3,4,10,[<i,3*(j-1)+i,R3.j>: i,j in [1..3]] cat [<4,10,1>]);
   Phi0:=Transpose(P)*M*P;
   mons:=MonomialsOfDegree(R3,3);
   cc:=[MonomialCoefficient(U,mon): mon in mons];
   assert exists(t){i: i in [1..10] | cc[i] ne 0};
   subst:=[cc[t]*PP.i :i in [1..t] cat [t..9]];
   subst[t]:=-&+[cc[i]*subst[i]: i in [1..10]| i ne t]/cc[t];
   Phi1:=Matrix(10,10,[&+[MonomialCoefficient(Phi0[i,j],mons[k])*subst[k]
     : k in [1..10]] : i,j in [1..10]]);
   return Matrix(9,9,[Phi1[i,j]: i,j in [2..10]]);
end function;

function Cubic9(C3);
   k:=BaseField(C3);
   P2<x,y,z>:=AmbientSpace(C3);
   U:=Equation(C3);
   R:=CoordinateRing(P2);
   P8<x1,x2,x3,x4,x5,x6,x7,x8,x9>:=ProjectiveSpace(k,8);
   mons:=MonomialsOfDegree(R,3);
   cc:=[MonomialCoefficient(U,mon): mon in mons];
   assert exists(t){i: i in [1..10] | cc[i] ne 0};
   phi:=map<C3->P8|[mons[i]: i in [1..10] | i ne t]>;
   C9:=Image(phi);
   Phi1,c:=KleinMatrix(C9);
   assert IsKleinMatrix(C9,Phi1);
   Phi2:=Phi9(C3,P8);
   assert IsKleinMatrix(C9,Phi2);
   B,c:=MatrixCongruence(Phi1,Phi2);
   return Phi1,Phi2,B,c;
end function;

//  8.3. Curves with a diagonal action of mu_n

function T5(la,tt,PP)
   assert la eq &*[t: t in tt];
   x:=&cat[[PP.i : i in [1..5]]: j in [1..2]];
   t:=tt cat tt;
   eqns:=[t[i]*x[i]^2+x[i+1]*x[i+4]-t[i+2]*t[i+3]*x[i+2]*x[i+3]
	   : i in [1..5]];
   return Scheme(PP,eqns);
end function;

function T7(la,tt,PP)
   assert la^4*(la-1) eq &*[t: t in tt];
   x:=&cat[[PP.i : i in [1..7]]: j in [1..2]];
   t:=tt cat tt;
   eqns1:=[t[i]*x[i]^2+x[i+1]*x[i+6]
      -(1/la^2)*t[i+2]*t[i+3]*t[i+4]*t[i+5]*x[i+2]*x[i+5] : i in [1..7]];
   eqns2:=[t[i]*x[i]^2+la*x[i+1]*x[i+6]
      -(1/la^3)*t[i+2]*t[i+3]^2*t[i+4]^2*t[i+5]*x[i+3]*x[i+4] : i in [1..7]];
   return Scheme(PP,eqns1 cat eqns2);
end function;

function PhiMu5(la,tt,PP)
   assert la eq &*[t: t in tt];
   x0,x1,x2,x3,x4:=Explode([PP.i:i in [1..5]]);
   t0,t1,t2,t3,t4:=Explode(tt);
   ll:=[t1*x1,x2,-x3,-t4*x4,t3*x3,x4,-x0,t0*x0,x1,t2*x2];
   return AlternatingMatrix(5,ll);
end function;

function PhiMu7(la,tt,PP)
   assert la^4*(la-1) eq &*[t: t in tt];
   x0,x1,x2,x3,x4,x5,x6:=Explode([PP.i:i in [1..7]]);
   t0,t1,t2,t3,t4,t5,t6:=Explode(tt);
   ll:=[t0*t1*t2*x1,-la*t2*x2,-la^2*x3,la^2*x4,la*t5*x5,-t5*t6*t0*x6,
	  t2*t3*t4*x3,-la*t4*x4,-la^2*x5,la^2*x6,la*t0*x0,
	  t4*t5*t6*x5,-la*t6*x6,-la^2*x0,la^2*x1,
	  t6*t0*t1*x0,-la*t1*x1,-la^2*x2,
	  t1*t2*t3*x2,-la*t3*x3,
	  t3*t4*t5*x4];
   return AlternatingMatrix(7,ll);
end function;

function Diagonal5(la,tt,k);
   assert la eq &*[t: t in tt];
   PP<x0,x1,x2,x3,x4>:=ProjectiveSpace(k,4);  
   C:=T5(la,tt,PP);
   Phi1,c:=KleinMatrix(C);
   assert IsKleinMatrix(C,Phi1);
   Phi2:=PhiMu5(la,tt,PP);
   assert IsKleinMatrix(C,Phi2);
   B,c:=MatrixCongruence(Phi1,Phi2);
   return Phi1,Phi2,B,c;
end function;

function Diagonal7(la,tt,k);
   assert la^4*(la-1) eq &*[t: t in tt];
   PP<x0,x1,x2,x3,x4,x5,x6>:=ProjectiveSpace(k,6);  
   C:=T7(la,tt,PP);
   Phi1,c:=KleinMatrix(C);
   assert IsKleinMatrix(C,Phi1);
   Phi2:=PhiMu7(la,tt,PP);
   assert IsKleinMatrix(C,Phi2);
   B,c:=MatrixCongruence(Phi1,Phi2);
   return Phi1,Phi2,B,c;
end function;

//  8.4. Numerical examples 

function Example275B3(P4)
   x1,x2,x3,x4,x5:=Explode([P4.i:i in [1..5]]);
   eqns:=[
     11*x1^2 - x2^2 + 5*x3^2 + x4^2 - 5*x5^2,
     2*x1*x4 + x2^2 + 10*x2*x3 - x2*x4 - 5*x2*x5 + 5*x3^2 - 5*x3*x4 - 5*x3*x5,
     2*x1*x5 + x2^2 + 2*x2*x3 + x2*x4 + x2*x5 + 5*x3^2 + x3*x4 + 5*x3*x5,
     11*x1*x2 + 55*x1*x3 + 2*x2*x4 - 10*x3*x5 + 2*x4^2 + 10*x5^2,
     11*x1*x2 + 11*x1*x3 + 2*x2*x5 - 2*x3*x4 - 4*x4*x5
   ];
  ll:=[-2*x5,2*x4,-x2-x3,x2+5*x3,-11*x1,x3+x5,x2+x4,-x2+x4,-5*x3+5*x5,-x1];
  return Scheme(P4,eqns),AlternatingMatrix(5,ll);
end function;

function Example37A1(P6)
   x1,x2,x3,x4,x5,x6,x7:=Explode([P6.i:i in [1..7]]);
   E<x,y,z>:=EllipticCurve([0,0,1,-1,0]);
   phi:=map<E->P6|[x*z^3,x^2*z^2,x*y*z^2,x^3*z,x^2*y*z,x^4,(y+z)*z^3]>;
   ll:=[x1,-x2,0,x2,x4,-x1+x3,x7,x2,-2*x7,x1-x3,-x2,-x4,x1+x3,-2*x2,x1+x7,
        -x4,-x6,x2-x5,3*x2+x5,-2*x1+x4-2*x7,x1-x2-3*x3];
   return E,Image(phi),AlternatingMatrix(7,ll);
end function;
   
function DivisorExample(D)
   n:=Degree(D);
   phi:=DivisorMap(D);
   C:=phi(Curve(D));
   if n eq 5 then 
     PP<x1,x2,x3,x4,x5>:=AmbientSpace(C); end if;
   if n eq 7 then 
     PP<x1,x2,x3,x4,x5,x6,x7>:=AmbientSpace(C); end if;
   if n eq 9 then 
     PP<x1,x2,x3,x4,x5,x6,x7,x8,x9>:=AmbientSpace(C); end if;
   if n eq 11 then 
     PP<x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11>:=AmbientSpace(C); end if;
   Phi,c:=KleinMatrix(C);
   assert IsKleinMatrix(C,Phi);
   return Phi,c;
end function;

//  Evidence for Conjecture 6.8 
  
function SecantExperiment(Phi);
   R:=Parent(Phi[1,1]);
   pp:=PfSubmax(Phi);
   n:=#pp;
   J:=Matrix(n,n,[Derivative(pp[i],R.j):i,j in [1..n]]);
   print "Computing determinant of J = (d p_i/d x_j) ....";
   ss:=Determinant(J);
   print "Factoring ....";
   ff:=Factorization(ss);
   if #ff eq 1 and ff[1,2] eq 2 then 
      print "We find that det(J) is a square!"; 
   end if;
   return ss,ff;
end function;

//  Sample code running the above examples
//  (with lots of testing/redundancy).

/*

E<x,y,z>:=EllipticCurve([1,1,1,-43,89]);

printf "\nEmbedding the %o by |5.0| we obtain ",E;
printf "a curve in P^4 with Klein matrix\n";
Phi1,Phi2,B,c:=Weierstrass5(E);
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

printf "\nEmbedding the %o by |7.0| we obtain ",E;
printf "a curve in P^6 with Klein matrix\n";
Phi1,Phi2,B,c:=Weierstrass7(E);
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

printf "\nEmbedding the %o by |9.0| we obtain ",E;
printf "a curve in P^8 with Klein matrix\n";
Phi1,Phi2,B,c:=Cubic9(E);
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

P2<x,y,z>:=ProjectiveSpace(Rationals(),2);
C:=Curve(P2,3*x^3+2*x^2*y-5*x^2*z-6*x*y^2+9*x*y*z-9*y^3+2*y*z^2+z^3);
printf "\nEmbedding the %o by |3.H| we obtain ",C;
printf "a curve in P^8 with Klein matrix\n";
Phi1,Phi2,B,c:=Cubic9(C);
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

tt:=[1,1,2,3,5];la:=30;
printf "\nT[%o;%o] has Klein matrix\n",la,tt;
Phi1,Phi2,B,c:=Diagonal5(la,tt,Rationals());
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

tt:=[1,1,1,1,2^4,3^4,5];la:=6;
printf "\nT[%o;%o] has Klein matrix\n",la,tt;
Phi1,Phi2,B,c:=Diagonal7(la,tt,Rationals());
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;

P4<x1,x2,x3,x4,x5>:=ProjectiveSpace(Rationals(),4);
C,Phi2:=Example275B3(P4);
printf "\nThe %o\nhas Klein matrix\n",C;
Phi1,c:=KleinMatrix(C);
assert IsKleinMatrix(C,Phi1);
assert IsKleinMatrix(C,Phi2);
B,c:=MatrixCongruence(Phi1,Phi2);
print Phi2;

P6<x1,x2,x3,x4,x5,x6,x7>:=ProjectiveSpace(Rationals(),6);
E,C,Phi2:=Example37A1(P6);
printf "\nEmbedding the %o by |6.0+P| where P=(0,0) we obtain ",E;
printf "a curve in P^6 with Klein matrix\n";
Phi1,c:=KleinMatrix(C);
assert IsKleinMatrix(C,Phi1);
assert IsKleinMatrix(C,Phi2);
B,c:=MatrixCongruence(Phi1,Phi2);
print Phi2;

E<x,y,z>:=EllipticCurve([0,0,1,-1,0]);
P2:=AmbientSpace(E);
C:=Curve(P2,Equation(E));
O:=C![0,1,0];
P:=C![0,0,1];
for n in [5..11 by 2] do
  D:=(n-1)*Divisor(O)+Divisor(P);
  printf "\nEmbedding the %o by |%o.0+P| where P=(0,0) we obtain ",E,n-1;
  printf "a curve in P^%o with Klein matrix\n",n-1;
  Phi,c:=DivisorExample(D);
  assert c ne 0;
  print Phi;
end for;

E<x,y,z>:=EllipticCurve([1,1,1,-5,2]);
printf "\nEmbedding the %o by |7.0| we obtain ",E;
printf "a curve in P^6 with Klein matrix\n";
Phi1,Phi2,B,c:=Weierstrass7(E);
assert Phi1 eq c*Transpose(B)*Phi2*B;
print Phi2;
ss,ff:=SecantExperiment(Phi2);

*/

////////////////////////////////////////////////////////////////////////

/*

Magma V2.11-2     Tue Sep  7 2004 16:29:27    [Seed = 810511170]
Type ? for help.  Type <Ctrl>-D to quit.
> load "kleinmatrix.m";

Embedding the Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 43*x + 89 
over Rational Field by |5.0| we obtain a curve in P^4 with Klein matrix
I(Sec^1 C) generated by 5 2-ics
[0   89*x1 + x4   x5   -43*x1 - x3 + x4   -x3]
[-89*x1 - x4   0   -x4   -x1 - x3   x2]
[-x5   x4   0   -x2   0]
[43*x1 + x3 - x4   x1 + x3   x2   0   -x1]
[x3   -x2   0   x1   0]

Embedding the Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 43*x + 89 
over Rational Field by |7.0| we obtain a curve in P^6 with Klein matrix
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
[0   x7   0   -89*x1 + 43*x2 + x3 - x4 + x5 - x6   x3   x5   0]
[-x7   0   x6   89*x2 - 43*x4 - x5   -x5   0   x4]
[0   -x6   0   x5   -x2   -x4   0]
[89*x1 - 43*x2 - x3 + x4 - x5 + x6   -89*x2 + 43*x4 + x5   -x5   0   x2 - x3 + 
    x4   -89*x1 + 43*x2 + x3   -x3]
[-x3   x5   x2   -x2 + x3 - x4   0   x3   x1]
[-x5   0   x4   89*x1 - 43*x2 - x3   -x3   0   x2]
[0   -x4   0   x3   -x1   -x2   0]

Embedding the Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 43*x + 89 
over Rational Field by |9.0| we obtain a curve in P^8 with Klein matrix
I(Sec^1 C) generated by 27 2-ics
I(Sec^2 C) generated by 30 3-ics
I(Sec^3 C) generated by 9 4-ics
[0   x2 - x4 - 43*x5 - x7 - x8 + 89*x9   x4   0   -x1   x5   0   -x2   -x4 - x5]
[-x2 + x4 + 43*x5 + x7 + x8 - 89*x9   0   -x3   x1   0   -x4   x2   0   89*x5]
[-x4   x3   0   -x7   x6   0   -x8   x7   x1 + x4 - x7 - 43*x8]
[0   -x1   x7   0   -x3   x8   0   -x4   -x7 - x8]
[x1   0   -x6   x3   0   -x7   x4   0   89*x8]
[-x5   x4   0   -x8   x7   0   -x9   x8   x2 + x5 - x8 - 43*x9]
[0   -x2   x8   0   -x4   x9   0   -x5   -x8 - x9]
[x2   0   -x7   x4   0   -x8   x5   0   89*x9]
[x4 + x5   -89*x5   -x1 - x4 + x7 + 43*x8   x7 + x8   -89*x8   -x2 - x5 + x8 + 
    43*x9   x8 + x9   -89*x9   0]

Embedding the Curve over Rational Field defined by
3*x^3 + 2*x^2*y - 5*x^2*z - 6*x*y^2 + 9*x*y*z - 9*y^3 + 2*y*z^2 + z^3 by |3.H| 
we obtain a curve in P^8 with Klein matrix
I(Sec^1 C) generated by 27 2-ics
I(Sec^2 C) generated by 30 3-ics
I(Sec^3 C) generated by 9 4-ics
[0   -2*x1 + 5*x2 + 6*x3 - 9*x4 + 9*x6 - 2*x8 - x9   -3*x4   0   3*x1   -3*x5   
    0   3*x2   -27*x3 + 6*x5]
[2*x1 - 5*x2 - 6*x3 + 9*x4 - 9*x6 + 2*x8 + x9   0   3*x3   -3*x1   0   3*x4   
    -3*x2   0   3*x5]
[3*x4   -3*x3   0   3*x7   -3*x6   0   3*x8   -3*x7   9*x1 + 6*x3 - 15*x4 - 
    18*x6 + 27*x7]
[0   3*x1   -3*x7   0   3*x3   -3*x8   0   3*x4   -27*x6 + 6*x8]
[-3*x1   0   3*x6   -3*x3   0   3*x7   -3*x4   0   3*x8]
[3*x5   -3*x4   0   3*x8   -3*x7   0   3*x9   -3*x8   9*x2 + 6*x4 - 15*x5 - 
    18*x7 + 27*x8]
[0   3*x2   -3*x8   0   3*x4   -3*x9   0   3*x5   -27*x7 + 6*x9]
[-3*x2   0   3*x7   -3*x4   0   3*x8   -3*x5   0   3*x9]
[27*x3 - 6*x5   -3*x5   -9*x1 - 6*x3 + 15*x4 + 18*x6 - 27*x7   27*x6 - 6*x8   
    -3*x8   -9*x2 - 6*x4 + 15*x5 + 18*x7 - 27*x8   27*x7 - 6*x9   -3*x9   0]

T[30;[ 1, 1, 2, 3, 5 ]] has Klein matrix
I(Sec^1 C) generated by 5 2-ics
[    0    x1    x2   -x3 -5*x4]
[  -x1     0  3*x3    x4   -x0]
[  -x2 -3*x3     0    x0    x1]
[   x3   -x4   -x0     0  2*x2]
[ 5*x4    x0   -x1 -2*x2     0]

T[6;[ 1, 1, 1, 1, 16, 81, 5 ]] has Klein matrix
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
[       0       x1    -6*x2   -36*x3    36*x4   486*x5  -405*x6]
[     -x1        0    16*x3   -96*x4   -36*x5    36*x6     6*x0]
[    6*x2   -16*x3        0  6480*x5   -30*x6   -36*x0    36*x1]
[   36*x3    96*x4 -6480*x5        0     5*x0    -6*x1   -36*x2]
[  -36*x4    36*x5    30*x6    -5*x0        0       x2    -6*x3]
[ -486*x5   -36*x6    36*x0     6*x1      -x2        0  1296*x4]
[  405*x6    -6*x0   -36*x1    36*x2     6*x3 -1296*x4        0]

The Scheme over Rational Field defined by
11*x1^2 - x2^2 + 5*x3^2 + x4^2 - 5*x5^2,
2*x1*x4 + x2^2 + 10*x2*x3 - x2*x4 - 5*x2*x5 + 5*x3^2 - 5*x3*x4 - 5*x3*x5,
2*x1*x5 + x2^2 + 2*x2*x3 + x2*x4 + x2*x5 + 5*x3^2 + x3*x4 + 5*x3*x5,
11*x1*x2 + 55*x1*x3 + 2*x2*x4 - 10*x3*x5 + 2*x4^2 + 10*x5^2,
11*x1*x2 + 11*x1*x3 + 2*x2*x5 - 2*x3*x4 - 4*x4*x5
has Klein matrix
I(Sec^1 C) generated by 5 2-ics
[           0        -2*x5         2*x4     -x2 - x3    x2 + 5*x3]
[        2*x5            0       -11*x1      x3 + x5      x2 + x4]
[       -2*x4        11*x1            0     -x2 + x4 -5*x3 + 5*x5]
[     x2 + x3     -x3 - x5      x2 - x4            0          -x1]
[  -x2 - 5*x3     -x2 - x4  5*x3 - 5*x5           x1            0]

Embedding the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field by
|6.0+P| where P=(0,0) we obtain a curve in P^6 with Klein matrix
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
[0   x1   -x2   0   x2   x4   -x1 + x3]
[-x1   0   x7   x2   -2*x7   x1 - x3   -x2]
[x2   -x7   0   -x4   x1 + x3   -2*x2   x1 + x7]
[0   -x2   x4   0   -x4   -x6   x2 - x5]
[-x2   2*x7   -x1 - x3   x4   0   3*x2 + x5   -2*x1 + x4 - 2*x7]
[-x4   -x1 + x3   2*x2   x6   -3*x2 - x5   0   x1 - x2 - 3*x3]
[x1 - x3   x2   -x1 - x7   -x2 + x5   2*x1 - x4 + 2*x7   -x1 + x2 + 3*x3   0]

Embedding the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field by
|4.0+P| where P=(0,0) we obtain a curve in P^4 with Klein matrix
I(Sec^1 C) generated by 5 2-ics
[       0 -x2 - x5       x4  x1 + x5      -x2]
[ x2 + x5        0 -x1 + x5      -x3       x4]
[     -x4  x1 - x5        0       x4      -x5]
[-x1 - x5       x3      -x4        0        0]
[      x2      -x4       x5        0        0]

Embedding the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field by
|6.0+P| where P=(0,0) we obtain a curve in P^6 with Klein matrix
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
[0   x7   -x6   0   x6   x5   x2 - x7]
[-x7   0   x3   x6   -2*x3   -x2 + x7   -x6]
[x6   -x3   0   -x5   x2 + x7   -2*x6   x3 + x7]
[0   -x6   x5   0   -x5   -x4   -x1 + x6]
[-x6   2*x3   -x2 - x7   x5   0   x1 + 3*x6   -2*x3 + x5 - 2*x7]
[-x5   x2 - x7   2*x6   x4   -x1 - 3*x6   0   -3*x2 - x6 + x7]
[-x2 + x7   x6   -x3 - x7   x1 - x6   2*x3 - x5 + 2*x7   3*x2 + x6 - x7   0]

Embedding the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field by
|8.0+P| where P=(0,0) we obtain a curve in P^8 with Klein matrix
I(Sec^1 C) generated by 27 2-ics
I(Sec^2 C) generated by 30 3-ics
I(Sec^3 C) generated by 9 4-ics
[0   0   -x9   x8   -x4   0   x8   -x7   x3 + x9]
[0   0   x8   -x7   x3 + x9   0   -x7   x6   -x2 - x8]
[x9   -x8   0   -x3 + x9   x8   x7   -2*x9   x2 + x8   -x7 - 2*x9]
[-x8   x7   x3 - x9   0   -x4 - x9   -x6   -x2 + 3*x8   -2*x7   x3 + 3*x8 + x9]
[x4   -x3 - x9   -x8   x4 + x9   0   x2 + x8   -2*x4 + x7   x3 - x8 + x9   
    -2*x4]
[0   0   -x7   x6   -x2 - x8   0   x6   -x5   x1 + x7]
[-x8   x7   2*x9   x2 - 3*x8   2*x4 - x7   -x6   0   -x1 + x7   -2*x3 + x6 + 
    2*x8 - 2*x9]
[x7   -x6   -x2 - x8   2*x7   -x3 + x8 - x9   x5   x1 - x7   0   x2 - 3*x7 + x8]
[-x3 - x9   x2 + x8   x7 + 2*x9   -x3 - 3*x8 - x9   2*x4   -x1 - x7   2*x3 - x6 
    - 2*x8 + 2*x9   -x2 + 3*x7 - x8   0]

Embedding the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field by
|10.0+P| where P=(0,0) we obtain a curve in P^10 with Klein matrix
I(Sec^1 C) generated by 44 2-ics
I(Sec^2 C) generated by 77 3-ics
I(Sec^3 C) generated by 55 4-ics
I(Sec^4 C) generated by 11 5-ics
[0   -x11   x10   0   -x10   -x9   -x4 + x11   0   -3*x9   x8   x3 - x10]
[x11   0   -x5   -x10   2*x5   x4 - x11   x10   x9   2*x4 + 4*x11   -x3 - 5*x10 
    2*x5 - x9]
[-x10   x5   0   x9   -x4 - x11   2*x10   -x5 - x11   -x8   x3 - x10   4*x9   
    -x4 + x10 - x11]
[0   x10   -x9   0   x9   x8   x3 - x10   0   3*x8   -x7   -x2 + x9]
[x10   -2*x5   x4 + x11   -x9   0   -x3 - 3*x10   2*x5 - x9 + 2*x11   x8   -4*x3
    - 2*x10   x2 - 3*x9   x8 - 2*x10]
[x9   -x4 + x11   -2*x10   -x8   x3 + 3*x10   0   3*x4 + x10 - x11   x7   -x2 + 
    7*x9   -6*x8   -x3 - x9 + 3*x10]
[x4 - x11   -x10   x5 + x11   -x3 + x10   -2*x5 + x9 - 2*x11   -3*x4 - x10 + x11
    0   x2 - x9   -x8 - 2*x10 - 4*x11   -3*x3 + x9 + 5*x10   -2*x5 + 2*x9 - 
    2*x11]
[0   -x9   x8   0   -x8   -x7   -x2 + x9   0   -3*x7   x6   x1 - x8]
[3*x9   -2*x4 - 4*x11   -x3 + x10   -3*x8   4*x3 + 2*x10   x2 - 7*x9   x8 + 
    2*x10 + 4*x11   3*x7   0   -x1 - 11*x8   6*x3 - x7 - 2*x9 + 2*x10]
[-x8   x3 + 5*x10   -4*x9   x7   -x2 + 3*x9   6*x8   3*x3 - x9 - 5*x10   -x6   
    x1 + 11*x8   0   -5*x2 + x8 + 3*x9]
[-x3 + x10   -2*x5 + x9   x4 - x10 + x11   x2 - x9   -x8 + 2*x10   x3 + x9 - 
    3*x10   2*x5 - 2*x9 + 2*x11   -x1 + x8   -6*x3 + x7 + 2*x9 - 2*x10   5*x2 - 
    x8 - 3*x9   0]

Embedding the Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 5*x + 2 over
Rational Field by |7.0| we obtain a curve in P^6 with Klein matrix
I(Sec^1 C) generated by 14 2-ics
I(Sec^2 C) generated by 7 3-ics
[0   x7   0   -2*x1 + 5*x2 + x3 - x4 + x5 - x6   x3   x5   0]
[-x7   0   x6   2*x2 - 5*x4 - x5   -x5   0   x4]
[0   -x6   0   x5   -x2   -x4   0]
[2*x1 - 5*x2 - x3 + x4 - x5 + x6   -2*x2 + 5*x4 + x5   -x5   0   x2 - x3 + x4   
    -2*x1 + 5*x2 + x3   -x3]
[-x3   x5   x2   -x2 + x3 - x4   0   x3   x1]
[-x5   0   x4   2*x1 - 5*x2 - x3   -x3   0   x2]
[0   -x4   0   x3   -x1   -x2   0]
Computing determinant of J = (d p_i/d x_j) ....
Factoring ....
We find that det(J) is a square!

Total time: 105.569 seconds, Total memory usage: 88.28MB

*/

