symplectic_basis.m

freeze;

import "misc.m": function_field;

function cup_product_matrix_split_places(basis, Q, g, r, W0 : prec := 15)

  d := Degree(Q);
  K := BaseRing(BaseRing(Q));

  // find the points at infinity:
  FF:=function_field(Q); // function field of curve over the rationals
  infplaces:=InfinitePlaces(FF);

  infplacesKinf := infplaces;
  Kinf := K;
  repeat 
    for Kinf_new in infplacesKinf do
      if not IsOne(Degree(Kinf_new)) then
        // field generated by points at infinity
        Kinf := Compositum(Kinf, NormalClosure(AbsoluteField(ResidueClassField(Kinf_new))));
      end if;
    end for;

    Kinfx:=RationalFunctionField(Kinf);
    Kinfxy:=PolynomialRing(Kinfx);
    FFKinf:=FunctionField(Kinfxy!Q); // function field of curve over Kinf

    infplacesKinf:=InfinitePlaces(FFKinf); // places at infinity all of degree 1, will be denoted by P
  until &and[IsOne(Degree(P)) : P in infplacesKinf];

  b0funKinf:=[]; // functions b^0 (finite integral basis)
  for i:=1 to d do
    b0i:=FFKinf!0;
    for j:=1 to d do
      b0i +:= Evaluate(W0[i,j],Kinfx.1)*FFKinf.1^(j-1);
    end for;
    b0funKinf[i]:=b0i;
  end for;

  L:=[];
  for i:=1 to 2*g do
    fun:=FFKinf!0;
    for j:=1 to d do
      fun +:= Evaluate(basis[i][j],Kinfx.1)*b0funKinf[j];
    end for;
    L[i]:=fun;
  end for;

  // compute the expansions omega_x

  omegax:=[]; // expansions of omega
  x_expansions_inf:=[];  // expansions of x

  for P in infplacesKinf do
    xfunx:=Expand(FFKinf!Kinfx.1,P : RelPrec:=prec+3); 
    //"Valuation of x", Valuation(xfunx); "L", L;
    dxdt:=Derivative(xfunx);
    zinv:=Expand(LeadingCoefficient(r)/(FFKinf!Evaluate(r,Kinfx.1)),P : RelPrec:=prec+3);
    omegaP:=[];
    for j:=1 to 2*g do
      omegaP[j]:=Expand(L[j],P : RelPrec:=prec+3)*dxdt*zinv;
    end for;
    Append(~omegax,omegaP);
    Append(~x_expansions_inf,xfunx);
  end for;

  C := ZeroMatrix(K, 2*g, 2*g);
  for i := 1 to 2*g do
    for j := 1 to 2*g do
      C[i,j] := &+[Coefficient(omegaP[i]*Integral(omegaP[j]), -1) : omegaP in omegax];
    end for;
  end for;
  return C, x_expansions_inf;
end function;


intrinsic CupProductMatrix(basis::SeqEnum[ModTupRngElt[RngUPol]], Q::RngUPolElt[RngUPol], g::RngIntElt, r::RngUPolElt, W0::AlgMatElt[FldFunRat] :
                           prec := 15, split := true)
  -> AlgMatElt, SeqEnum[RngSerLaurElt]
  {Given a basis of H^1_dR of a smooth projective curve such that the
  first g elements generate regular differentials, computes the cup product.}
  // split determines whether we work over an extension
  // where all places split. This is usually faster.

  //if Max([Degree(P) : P in infplaces]) le 1 then
  // Usually faster to work over extension
  if split then 
    return cup_product_matrix_split_places(basis, Q, g, r, W0 : prec := prec);
  end if;

  d := Degree(Q);
  K := BaseRing(BaseRing(Q));

  // find the points at infinity:
  FF:=function_field(Q); // function field of curve over K
  infplaces:=InfinitePlaces(FF);

  Kx := RationalFunctionField(K);
  Kxy := PolynomialRing(Kx);
  r0 := K!LeadingCoefficient(r);

  b0funK:=[]; // functions b^0 (finite integral basis)
  for i:=1 to d do
    b0i:=FF!0;
    for j:=1 to d do
      b0i +:= Evaluate(W0[i,j],Kx.1)*FF.1^(j-1);
    end for;
    b0funK[i]:=b0i;
  end for;

  L:=[];
  for i:=1 to 2*g do
    fun:=FF!0;
    for j:=1 to d do
      fun +:= Evaluate(basis[i][j],Kx.1)*b0funK[j];
    end for;
    L[i]:=fun;
  end for;

  omegax := [* *];
  x_expansions_inf := [* *];  // expansions of x

  for P in infplaces do
    xfunx:=Expand(Kx.1,P : RelPrec:=prec+3);
    Append(~x_expansions_inf, xfunx);
    dxdt:=Derivative(xfunx);
    
    zinv:=Expand(r0/(Evaluate(r,Kx.1)),P : RelPrec:=prec+3);

    omegaP:=[];
    for j:=1 to 2*g do
      omegaP[j]:=Expand(L[j],P : RelPrec:=prec+3)*dxdt*zinv;
    end for;
    Append(~omegax, omegaP);
  end for;

  C := ZeroMatrix(K, 2*g, 2*g);
  for i := 1 to 2*g do
    for j := 1 to 2*g do
      if K eq Rationals() then
        C[i,j] := &+[Trace(Coefficient(omegaP[i]*Integral(omegaP[j]), -1)) : omegaP in omegax];
      else
        C[i,j] := &+[Trace(Coefficient(omegaP[i]*Integral(omegaP[j]), -1), K) : omegaP in omegax];
      end if;
    end for;
  end for;
  return C, [elt : elt in x_expansions_inf];
end intrinsic;

// For coercing basis to sequence of vectors if need be
intrinsic CupProductMatrix(basis::SeqEnum[SeqEnum[RngUPol]], Q::RngUPolElt[RngUPol], g::RngIntElt, r::RngUPolElt, W0::AlgMatElt[FldFunRat] :
                           prec := 15, split := true)
  -> AlgMatElt, SeqEnum[RngSerLaurElt]
  {Given a basis of H^1_dR of a smooth projective curve such that the
  first g elements generate regular differentials, computes the cup product.}
  return CupProductMatrix([Vector(v) : v in basis], Q, g, r, W0 : prec := prec, split := split);
end intrinsic;


intrinsic SymplecticBasisH1(C::AlgMatElt) -> ModMatFldElt, ModTupFld
  {Given the cup product matrix of H^1_dR of a smooth projective curve, computes
  a symplectic basis with respect to the cup product.}
  // C  is the cup product matrix for some basis (omega_i)_{i=1,..,2g} of H^1 such that
  // omega_1,...,omega_g are holomorphic.
  // Compute linear combinations eta_1,...,eta_g of the omega_i such that
  // omega_1,...,omega_g,eta_1,...,eta_g form a symplectic basis of H^1

  assert C eq -Transpose(C);
  assert Nrows(C) eq Ncols(C);
  assert IsEven(Nrows(C));
  g := Nrows(C) div 2;
  assert IsZero(ExtractBlock(C, 1, 1, g, g));
  R := BaseRing(C);
  // The matrix A will contain the coefficients
  A := ZeroMatrix(R, g, 2*g);

  // First determine a_ij for g<j<2g+1
  C_small := ExtractBlock(C, 1, g+1, g, g);
  M := ZeroMatrix(R, g^2, g^2);

  for i := 0 to g-1 do
    InsertBlock(~M, C_small, i*g+1, i*g+1);
  end for;

  z := [1] cat [0: i in [1..g^2-1]];
  n := 1;
  for i := 1 to g-1 do
    n +:= g+1;
    z[n] := 1;
  end for;
  z := Vector(R, g^2, z);
  // M*v = z encodes [omega_i] cup [eta_j] = delta_ij
  v, L := Solution(Transpose(M), z);

  k := 0;
  for i := 1 to g do
    for j := g+1 to 2*g do
      k +:= 1;
      A[i,j] := v[k];
    end for;
  end for;
  
  function row(i,j)
    // This represents the row index of the equation [eta_i] cup [eta_j] =0 
    // in the matrix N below.
    return &+([0] cat [g-k : k in [1..i-1]]) + j-i;
  end function;

  // Now determine a_ij for 1<j<g+1
  N := ZeroMatrix(R, g*(g-1) div 2, g^2);
  m := 1;

  for i := 1 to g do
    for j := 1 to g do
      n := (i-1)*g + j; // nth column corresponds to a_ij
      // Find equations [eta_*] cup [eta_*] = 0 having nontrivial a_ij-coeff
      for q := i+1 to g do
        // [eta_i] cup [eta_q] = 0 has nontrivial a_ij-coeff
        N[row(i,q),n] +:= &+[A[q,l]*C[j,l] : l in [g+1..2*g]];
      end for;
      for r := 1 to i-1 do
        // [eta_r] cup [eta_i] = 0 has nontrivial a_ij-coeff
        N[row(r,i),n] +:= &+[A[r,k]*C[k,j] : k in [g+1..2*g]];
      end for;
    end for;
  end for;

  y := [R!0: i in [1..g*(g-1) div 2]];
  n := 0;
  for i := 1 to g do
    for j := i+1 to g do
      n +:= 1;
      // All of these entries of A have already been fixed above
      y[n] := -&+[&+[ A[i,k]*A[j,l]*C[k,l] : l in [g+1..2*g]] : k in [g+1..2*g]]; 
    end for;
  end for;
  y := Vector(R, #y, y);

  // N^T*x = y encodes [eta_i] cup [eta_j] = 0
  x, K := Solution(Transpose(N), y);
  // K can be used to create alternative bases, if desired

  k := 0;
  for i := 1 to g do
    for j := 1 to g do
      k +:= 1;
      A[i,j] := x[k];
    end for;
  end for;
  return A, K;
end intrinsic;