hodge.m

freeze;

//////////////////////////////////////////////
// Functions for computing Hodge structures //
//////////////////////////////////////////////

import "auxpolys.m": auxpolys;
import "misc.m": function_field;

intrinsic HodgeData(Q::RngUPolElt[RngUPol], g::RngIntElt, W0::AlgMatElt[FldFunRat], basis::SeqEnum[ModTupRngElt[RngUPol]],
                    Z::AlgMatElt, bpt::PlcFunElt : r := auxpolys(Q), prec := 5)
  -> ModTupFldElt, ModTupFldElt, ModTupRngElt[RngUPol], RngIntElt
  {Compute the 1-form eta, as a vector of coefficients
  w.r.t. basis[i] for i=2g+1,...,2g+k-1 where k is the 
  number of points lying over x=infinity.}
  // prec is the relative precision to be used for the expansions
  // To do: We don't want to work over the splitting field!

  // Magma only allows 6 arguments in signature so we make r an optional argument
  require ISA(Type(r), RngUPolElt): "r must be of type RngUPolElt.";

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

  // find the points at infinity:

  Kx := RationalFunctionField(K);
  Kxy := PolynomialRing(Kx);

  FF := function_field(Q); // function field of curve over K
  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+#infplacesKinf-1) 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, Omega_x, b^0_x

  omegax:=[]; // expansions of omega
  Omegax:=[]; // expansions of Omega
  b0funx:=[]; // expansions of b^0
  xfunx:=[];  // expansions of x

  for i:=1 to #infplacesKinf do

    P:=infplacesKinf[i];
    xfunx[i]:=Expand(FFKinf!Kinfx.1,P : RelPrec:=prec+3); 
    dxdt:=Derivative(xfunx[i]);
    
    zinv:=Expand(LeadingCoefficient(r)/(FFKinf!Evaluate(r,Kinfx.1)),P : RelPrec:=prec+3);
    
    omegaP:=[];
    for j:=1 to 2*g+#infplacesKinf-1 do
      omegaP[j]:=Expand(L[j],P : RelPrec:=prec+3)*dxdt*zinv;
    end for;
    Append(~omegax,omegaP);
    
    OmegaP:=[];
    for j:=1 to 2*g do
      OmegaP[j]:=Integral(omegaP[j]); 
    end for;
    Append(~Omegax,OmegaP);

    b0funP:=[];
    for j:=1 to d do
      b0funP[j]:=Expand(b0funKinf[j],P : RelPrec:=prec+3);
    end for;
    Append(~b0funx,b0funP);

  end for;

  // compute expansions of Omega*Z*omega at all points at infinity

  omegaZOmega:=[];
  for i:=1 to #infplacesKinf do
    omegaZOmegaP:=0;
    for j:=1 to 2*g do
      for k:=1 to 2*g do
        omegaZOmegaP +:= omegax[i][j]*Z[j,k]*Omegax[i][k];
      end for;
    end for;
    Append(~omegaZOmega, omegaZOmegaP);
  end for;

  // set up the linear system eta*A=v satisfied by eta
  
  v:=[];
  A:=ZeroMatrix(Kinf,#infplacesKinf-1,#infplacesKinf);
  for i:=1 to #infplacesKinf do
    v[i]:=-Coefficient(omegaZOmega[i],-1); // residue of eta at i-th point of infinity 
    for j:=1 to #infplacesKinf-1 do
      A[j,i]:=Coefficient(omegax[i][2*g+j],-1); // residue of omega_{2g+j} at i-th point at infinity
    end for;
  end for;
  //

  eta:=Solution(A,Vector(v)); // solve for eta
  eta:=ChangeRing(eta,K);

  gx:=[]; // functions g_x
  for i:=1 to #infplacesKinf do
    dgxi:=omegaZOmega[i]; 
    for j:=1 to (#infplacesKinf-1) do
      dgxi +:= eta[j]*omegax[i][2*g+j]; 
    end for;
    gx[i]:=Integral(dgxi);
  end for;

  OmegaZs2Omega:=[];
  for i:=1 to #infplacesKinf do
    OmegaZs2OmegaP:=0;
    for j:=1 to 2*g do
      for k:=g+1 to 2*g do
        OmegaZs2OmegaP +:= Omegax[i][j]*Z[j,k]*Omegax[i][k]; 
      end for;
    end for;
    OmegaZs2Omega[i]:=OmegaZs2OmegaP;
  end for;

  for i:=1 to #infplacesKinf do
    gx[i] +:= OmegaZs2Omega[i];
  end for;

  poleorder:=0;
  for i:=1 to #infplacesKinf do
    for j:=1 to 2*g do
      poleorder:=Minimum(poleorder,Valuation(Omegax[i][j]));
    end for;
  end for;
  poleorder_Omegax := poleorder;

  for i:=1 to #infplacesKinf do
    val:=Valuation(gx[i]);
    poleorder:=Minimum(poleorder,val);
  end for;

  done:=false;
  degx:=0;
  while not done do // try larger and larger degree in x
    
    for i:=1 to #infplacesKinf do
      for j:=1 to d do
        poleorder:=Minimum(poleorder,Valuation(b0funx[i][j])+degx*Valuation(xfunx[i]));
      end for;
    end for;

    v:=[]; // coefficients of principal parts of all gx 
    cnt:=0;
    for i:=1 to #infplacesKinf do
      for j:=poleorder to -1 do
        cnt +:= 1;
        v[cnt]:=Coefficient(gx[i],j);
      end for;
    end for;

    rows:=[];

    for i:=1 to g do
      row:=[];
      cnt:=0;
      for j:=1 to #infplacesKinf do
        for k:=poleorder to -1 do
          cnt +:= 1;
          row[cnt]:=Coefficient(Omegax[j][i+g],k); // coefficients of principal part of Omegax_{i+g} at jth point at infinity
        end for; 
      end for;
      Append(~rows,row);
    end for;

    for i:=1 to d do
      for j:=0 to degx do
        row:=[];
        cnt:=0;
        for k:=1 to #infplacesKinf do
          for l:=poleorder to -1 do
            cnt +:= 1;
            row[cnt]:=Coefficient(b0funx[k][i]*xfunx[k]^j,l); // coefficients of principal part of x^j*b^0_i at kth point at infinity
          end for;
        end for;
        Append(~rows,row);  
      end for;
    end for;   
      
    suc,sol:=IsConsistent(Matrix(rows),Vector(v));
    if suc then
      done:=true;
    else // if no success, increase the degree in x
      degx +:= 1;
    end if;
  
  end while;

  // read off beta from solution

  beta:=[];
  for i:=1 to g do
    beta[i] := K!sol[i];
  end for;

  // read off gamma from solution

  Kt := PolynomialRing(K);
  gamma:=[];
  cnt:=g;
  for i:=1 to d do
    poly:=Kt!0;
    for j:=0 to degx do
      cnt +:= 1;
      poly +:= (K!sol[cnt])*Kt.1^j;
    end for;
    Append(~gamma, poly);
  end for;

  b0fun:=[]; // 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;
    b0fun[i]:=b0i;
  end for;

  // substract constant such that gamma(bpt)=0

  gamma_FF:=FF!0;
  for i:=1 to d do
    gamma_FF +:= Evaluate(gamma[i],Kx.1)*b0fun[i];
  end for;
  gamma[1] -:= Evaluate(gamma_FF,bpt); 
  //"gamma_FF", gamma_FF;
  //"gamma", gamma;
  //Divisor(gamma_FF);
  //b0fun;


  return Vector(eta),Vector(beta),Vector(gamma),Integers()!poleorder_Omegax;

end intrinsic;