/*****************************************
 * Routines to work with divisor classes & jacobians of hyperelliptic curves
 *****************************************/

// import "/home/magma/package/Geometry/Algaem/localdata.m":GetLocalData,StoreLocalData;
declare attributes JacHyp:Curve;
import "localdata.m":GetLocalData,StoreLocalData;

intrinsic Reduction(C::CrvHyp,p::RngIntElt)->CrvHyp,Map
  {returns a reduction for a hyperelliptic curve}
  require BaseRing(C) eq Rationals():
    "C must be over the rationals";
  if not(assigned GetLocalData(C,p)`reduction) then
    Cr:=BaseChange(C,Coercion(Rationals(),GF(p)));
    w:=Gradings(C)[1][2];
    redmap:=function(a)
      if Valuation(a[3],p) gt Valuation(a[1],p) then
        return Cr![GF(p)!c: c in [1,a[2]/a[1]^w,a[3]/a[1]]];
      else
        return Cr![GF(p)!c: c in [a[1]/a[3],a[2]/a[3]^w,1]];
      end if;
    end function;
    rc:=GetLocalData(C,p);
    rc`reduction:=map<C(Rationals())->Cr(GF(p))| p:->redmap(p)>;
    StoreLocalData(C,p,rc);
  end if;
  rc:=GetLocalData(C,p);
  return Codomain(rc`reduction),rc`reduction;
end intrinsic;

intrinsic Reduction(J::JacHyp,p::RngIntElt)->JacHyp,Map
  {}
  require BaseRing(J) eq Rationals():
     "J must be over the rationals";
  if not(assigned GetLocalData(J,p)`reduction) then
    Cr:=Reduction(Curve(J),p);
    CrIsOdd:=IsOdd(Degree(HyperellipticPolynomials(Cr)));
    Jr:=Jacobian(Cr);
    FF:=BaseRing(Jr);
    PQ:=PolynomialRing(BaseRing(J));
    Pr:=PolynomialRing(FF);
    PQtoPr:=hom<PQ->Pr|
          Coercion(BaseRing(PQ),BaseRing(Pr))*Coercion(BaseRing(Pr),Pr),Pr.1>;
    red:=function(P)
      if P[3] eq 0 then
        return Jr!0;
      end if;
      xpol:=P[1];
      yrel:=P[2];

      x:=Parent(xpol).1;
      cfxpol:=Coefficients(xpol);
      v:=Infinity();
      ixv:=0;
      for i in [1..#cfxpol] do
        if Valuation(cfxpol[i],p) lt v then
          v:=Valuation(cfxpol[i],p);
          ixv:=i;
        end if;
      end for;
      xpol:=xpol/p^v;
      cfxpol:=Coefficients(xpol);
      print "ixv",ixv;
      //now we know that xpol has integral coefficients at p and that
      //x^(ixv-1) is the lowest degree monomial with coefficient that's a unit at p
      while true do
        cfyrel:=Coefficients(yrel);
        iyv:=0;
        v:=Infinity();
        for i in [1..#cfyrel] do
          if Valuation(cfyrel[i],p) lt v then
            v:=Valuation(cfyrel[i],p);
            iyv:=i;
          end if;
        end for;

        if v ge 0 then
          //we have an integral relation for y, so now xpol,yrel can be reduced
          xpolr:=PQtoPr(xpol);
          yrelr:=PQtoPr(yrel);
          if yrelr eq 0 then
            yrelr:=xpolr;
          end if;
          print "Trying to create",<xpolr,yrelr,P[3]>;
          if Degree(xpolr) eq 0 then
            print "Divisor in reduction entirely supported at infinity";
            f,h:=HyperellipticPolynomials(Cr);
            d:=Degree(f);
            if exists{i: i in [d-P[3]+1..d] |
                  Coefficient(yrelr^2+h*yrelr-f,i) ne 0} then
              print "but yrel does not have enough zeros there. Returning 0 divisor";
              return Jr!0;
            else
              print "yrel seems to define a legitimate point here";
            end if;
          end if;
          if CrIsOdd then
            return elt<Jr|xpolr,yrelr,Degree(xpolr)>;
          else
            return elt<Jr|xpolr,yrelr,P[3]>;
          end if;
        elif Valuation(cfxpol[P[3]+1],p) eq 0 then
          //something strange happens: xpol mod p has full degree,
          //so the reduction of the divisor is completely supported at
          //finite points. Yet, the relation for y is awkward.
          //We assume that the divisor in reduction is of the form
          //[(x,y)+(x,-y) - infty^+ - infty^-].
          //The least we can do is check that Discriminant(xpolr)=0
          xpolr:=PQtoPr(xpol);
          assert Discriminant(xpolr) eq 0;
          print "exceptional branch triggered.";
          return Jr!0;
        end if;

        //We have to do something about the coefficient of x^iyv in yrel
        //if iyv<ixy, we cannot reach the offending term.
        error if iyv lt ixv, "Cannot reach term";

        //the offending coefficient:
        c:=cfyrel[iyv];
        //take the unit part ...
        cunit:=c/p^Valuation(c,p);
        M:=(Integers()!(GF(p)!(cunit/cfxpol[ixv])))*p^Valuation(c,p);
        //*p^Valuation(c,p);
        print "M:",M;
        //this should kill one valuation level in the coefficient
        yrel:=yrel-M*xpol*x^(iyv-ixv);
        print "xpol,yrel:=",<xpol,yrel>;
      end while;
    end function;
    rc:=GetLocalData(J,p);
    rc`reduction:=map<J->Jr| P:-> red(P)>;
    StoreLocalData(J,p,rc);
  end if;
  rc:=GetLocalData(J,p);
  return Codomain(rc`reduction),rc`reduction;
end intrinsic;

intrinsic Reduction(MWmap::Map,P::RngIntElt)->Map,Map
  {}
  require ISA(Type(Domain(MWmap)),GrpAb) and
          ISA(Type(Codomain(MWmap)),JacHyp):
    "The map should be from an abelian group into an elliptic curve or its",
    "point set over the base ring";
  if not(assigned GetLocalData(MWmap,P)`reduction) then
    J:=Codomain(MWmap);
    Jr,toJr:=Reduction(J,P);
    JrGrp,JrGrpMap:=AbelianGroup(Jr);
    rc:=GetLocalData(MWmap,P);
    rc`reduction:=<
      hom<Domain(MWmap)->JrGrp|
        [g@MWmap@toJr@@JrGrpMap:g in OrderedGenerators(Domain(MWmap))]>,
      JrGrpMap>;
    StoreLocalData(MWmap,P,rc);
  end if;
  rc:=GetLocalData(MWmap,P);
  return rc`reduction[1],rc`reduction[2];
end intrinsic;
  
intrinsic Divisor(P::JacHypPt)->DivCrvElt
  {Returns a divisor representing the divisor class corresponding to P}
  C:=Curve(Parent(P));
  K:=FunctionField(C);
  f:=Evaluate(P[1],K.1);
  g:=K.2-Evaluate(P[2],K.1);
  Kr:=RationalExtensionRepresentation(K);
  Ofin:=MaximalOrderFinite(Kr);
  Oinf:=MaximalOrderInfinite(Kr);
  D:=Divisor(ideal<Ofin|f,g>);
  if Degree(P[1]) gt 0 and Degree(D) lt Degree(P[1]) then
    assert ((Degree(D)-Degree(P[1])) mod Degree(D)) eq 0;
    D:=((Degree(D)-Degree(P[1])) div Degree(D))*D;
  end if;
  if Degree(P[1]) lt P[3] then
    Z:=Divisor(ideal<Oinf|g> meet Oinf);
    Z:=Support(Z);
    assert #Z eq 1;
    Z:=Z[1];
    assert ((P[3]-Degree(D)) mod Degree(Z)) eq 0;
    D:=D+((P[3]-Degree(D)) div Degree(Z))*Z;
  end if;
  infty:=&+InfinitePlaces(Kr);
  assert (Degree(D) mod Degree(infty)) eq 0;
  return CurveDivisor(DivisorGroup(K)!(D-(Degree(D) div Degree(infty))*infty));
end intrinsic;

intrinsic JacobianPoint(D::DivCrvElt)->JacHypPt
  {Returns the point on the jacobian corresponding to the divisor class of D}
  C:=Curve(D);
  require ISA(Type(C),CrvHyp):
    "Only implemented for hyperelliptic curves";
  require Degree(D) eq 0:
    "Only degree 0 divisors can be mapped to Jacobian";
  g:=Genus(C);
  K:=FunctionField(C);
  Kr:=RationalExtensionRepresentation(K);
  Ofin:=MaximalOrderFinite(Kr);
  Oinf:=MaximalOrderInfinite(Kr);
  infty:=DivisorGroup(K)!(&+InfinitePlaces(Kr));
  require g mod 2 eq 0 or Degree(infty) eq 1:
    "For odd genus curves there must be a unique place of degree 1 at infinity";
  Dred:=Reduction(FunctionFieldDivisor(D)+(g div Degree(infty))*infty,infty);
  gens:=Generators(Ideals(DivisorGroup(Kr)!Dred));
  P:=PolynomialRing(BaseRing(C));

  if #gens eq 1 then
    a:=Evaluate(Numerator(RationalFunction(K!gens[1])),[P.1,0]);
    b:=Parent(a)!0;
  else
    a:=Evaluate(Numerator(RationalFunction(K!gens[1])),[P.1,0]);
    b:=Numerator(RationalFunction(K!gens[2]));
    assert Degree(b,2) eq 1;
    b:=Evaluate(-b+MonomialCoefficient(b,Parent(b).2)*Parent(b).2,[P.1,0]);
  end if;
  
  if Degree(a) eq Degree(Dred) then
    jpnt:=Jacobian(C)![a,b];
    assert FunctionFieldDivisor(Divisor(jpnt)) eq
                              Dred-(Degree(Dred) div Degree(infty))*infty;
    return jpnt;
  end if;
  
  assert #InfinitePlaces(K) eq 2;
  rrs,rrm:=RiemannRochSpace(-Dred+((Degree(Dred) div 2)+1)*infty);
  assert Dimension(rrs) eq 1;
  a:=rrm(rrs.1);
  a:=Evaluate(Numerator(RationalFunction(K!gens[1])),[P.1,0]);
  rrs,rrm:=RiemannRochSpace(-Dred+((Degree(Dred) div 2)+2)*infty);
  assert Dimension(rrs) eq 3;
  bl:=exists(b){b: f in Basis(rrs)| Degree(b,2) eq 1 where
          b:=Numerator(RationalFunction(rrm(f)))};
  require bl: "This should not happen.";
  b:=Evaluate(-b+MonomialCoefficient(b,Parent(b).2)*Parent(b).2,[P.1,0]);
  jpnt:=elt<Jacobian(C)|a,b,Degree(Dred)>;
  assert FunctionFieldDivisor(Divisor(jpnt)) eq
                         Dred-(Degree(Dred) div Degree(infty))*infty;
  return jpnt;
end intrinsic;
