/****************************************************************************
routines.m

by: Nils Bruin

date: August 15, 2003

Routines to accompany the article:

Nils Bruin, The primitive solutions to the equation x^3+y^9=z^2.
****************************************************************************/

/****************************************************************************
Some routines copied in from the standard source tree:
****************************************************************************/

//from Geometry/Algaem/localdata:

declare attributes Str:StrLocalData;

StrLocalDataRecFormat:=recformat<
    
    //for isogenies (2 and mult-by-2) of elliptic curves
    xlist,twotor,twotorimg,
    
    //for mwmaps (maps from abelian group to an elliptic curve)
    reduction,
    
    //for absolute algebras
    TwoSelmerMap,TwoSelmerMapComponents
    >;

function GetLocalData(S,p)
//  {Retrieve information stored on a structure, indexed by a prime ideal}
  if not(assigned(S`StrLocalData)) then
    S`StrLocalData:=[];
  end if;
  if exists(t){u[2]: u in S`StrLocalData | u[1] eq p} then
    return t;
  else
    return rec<StrLocalDataRecFormat|xlist:=[]>;
  end if;
end function;

procedure StoreLocalData(S, p, data) // (S::Str,p::RngOrdIdl,data::Rec)
//  {Store information on a structure, indexed by a prime ideal}
  if not(assigned(S`StrLocalData)) then
    S`StrLocalData:=[];
  end if;
  if exists(j){i:i in [1..#S`StrLocalData] | S`StrLocalData[i][1] eq p} then
    S`StrLocalData[j]:=<p,data>;
  else
    Append(~S`StrLocalData,<p,data>);
  end if;
end procedure;

//from Geometry/Arith/loctools.m:
function Minim(L)
//  {Computes minimum}
  M:=Infinity();
  for a in L do
    if ISA(Type(Universe(L)),Cop) then
      if Retrieve(a) lt M then
        M:=Retrieve(a);
      end if;
    else
      if a lt M then
        M:=a;
      end if;
    end if;
  end for;
  return M;
end function;

function Maxim(L)
//  {Computes maximum}
  M:=-Infinity();
  for a in L do
    if ISA(Type(Universe(L)),Cop) then
      if Retrieve(a) gt M then
        M:=Retrieve(a);
      end if;
    else
      if a gt M then
        M:=a;
      end if;
    end if;
  end for;
  return M;
end function;

function IntNInf()
//  {Structure that contains both the integers and infinity}
  return cop<PowerStructure(Infty),Integers()>;
end function;

function MinValuation(f)
//  {The minimum valuation of the coefficients of f}
  if ISA(Type(f),SeqEnum) then
    return Minim([IntNInf()|Valuation(c):c in f]); 
  elif ISA(Type(f),RngUPolElt) or ISA(Type(f),RngMPolElt) then  
    return Minim([IntNInf()|Valuation(c):c in Coefficients(f)]);
  else
    error "Argument not supported";
  end if;
end function;

function MinPrec(f)
//  {The minimum absolute precision in the coefficients of f}
  if ISA(Type(f),SeqEnum) then
    return Minim([IntNInf()|AbsolutePrecision(c):c in f]); 
  elif ISA(Type(f),RngUPolElt) or ISA(Type(f),RngMPolElt) then  
    return Minim([IntNInf()|AbsolutePrecision(c):c in Coefficients(f)]);
  else
    error "Argument not supported";
  end if;
end function;

function ShiftVal(F,v)
//  {Shift content}
  if F cmpeq 0 then
    return F;
  end if;
  R:=BaseRing(Parent(F));
  pi:=UniformizingElement(R);
  cf:=Coefficients(F);
  mn:=Monomials(F);
  if v lt 0 then
    if ISA(Type(Universe(cf)),FldLoc) then
      cf:=[c / pi^(-v):c in cf];
    else  
      cf:=[c div pi^(-v):c in cf];
    end if;
  else
    cf:=[c*pi^v:c in cf];
  end if;
  return &+[cf[i]*mn[i]:i in [1..#mn]];
end function;

function StripPrecisionlessTerms(P)
//  {Equate precisionless terms to 0}
  //nice to get polynomials that lost all precision to actually compare 0.
  cf:=Coefficients(P);
  if forall{c:c in cf| RelativePrecision(c) eq 0} then
    //#cf eq 0 then
    return Zero(Parent(P));
  else
    return &+[cf[i]*mn[i]:i in [1..#cf]|RelativePrecision(cf[i]) ne 0] where mn:=Monomials(P);
  end if;
end function;

function CoefficientByExponentVector(P,v)
  error if #v ne Rank(Parent(P)),
    "Vector must be of length equal to the rank of the polynomial ring";
  cfs:=Coefficients(P);
  mons:=[Exponents(m):m in Monomials(P)];
  if exists(i){i: i in [1..#mons] | mons[i] eq v} then
    return cfs[i];
  else
    return BaseRing(Parent(P))!0;
  end if;
end function;

function FlattenPrecision(P)
//  {flatten the absolute precision on the coefficients of P}
  //It's not always a good idea to have different precisions floating around on
  //the coefficients of a polynomial. This one flattens the precision to
  //the lowest value.
  cf:=Coefficients(P);
  mn:=Monomials(P);
  pi:=UniformizingElement(BaseRing(Parent(P)));
  minprec:=MinPrec(P);
  if minprec eq Infinity() then
    return P;
  else
    cf:=[c+O(pi^minprec):c in cf];
    return &+[cf[i]*mn[i]:i in [1..#cf]];
  end if;
end function;

function pAdicEvaluate(P,v)
  return &+[cf[i]*v^(i-1):i in [1.. #cf]] where cf:=Coefficients(P);
end function;

/****************************************************************************
Routines for Chabauty-arguments:
****************************************************************************/
       
declare verbose EllChab,1;

function DiscreteLog(EpGrpMap,g)
  grp:=Domain(EpGrpMap);
  if Ngens(grp) eq 1 then
    return Log(EpGrpMap(grp.1),g);
  elif Ngens(grp) eq 2 then
    error if not( Order(grp.2) mod Order(grp.1) eq 0), "Unexpected group order";
    g1:=EpGrpMap(grp.1);
    g2:=EpGrpMap(grp.2);
    n:=Order(grp.1); m:=Order(grp.2) div n;
    e:=Log(n*g2,n*g);
    mg2:=m*g2;
    bl:=exists(a){grp![i,m*j]:i,j in [0..n-1]|
                       i*g1+j*mg2 eq g0} where g0:=g-e*g2;
    error if not(bl), "Element does not lie in the specified group";
    assert EpGrpMap(e*grp.2+a) eq g;
    return e*grp.2+a;
  else
    error "Impossible group structure for elliptic curve";
  end if;
end function;

function ProjectiveReduction(L,p,Rp)
  pi:=SafeUniformiser(p);
  v:=Minimum([Minimum([Valuation(c,p):c in Coefficients(l)]):l in L]);
  L:=[pi^(-v)*l:l in L];
  
  R:=Universe(L);
  K:=BaseRing(R);
  k,tok:=ResidueClassField(p);
  assert BaseRing(Rp) eq k;
  assert IntegerRing(K) eq Order(p);
  assert Rank(R) eq Rank(Rp);

  toRp:=hom<Universe(L)->Rp|map<K->Rp|c:->Rp!(tok(c))>,OrderedGenerators(Rp)>;

  return [Rp|toRp(l):l in L];
end function;

intrinsic Reduction(MWmap::Map,P::RngOrdIdl)->Map,Map
  {Computes the map between Abelian groups induced by the reduction map.
   The codomain of the returned map is isomorphic to the group of rational
   points on the curve in reduction. Presently, the elliptic curve must have
   good reduction at P.}
   
  require ISA(Type(Domain(MWmap)),GrpAb) and
          (ISA(Type(Codomain(MWmap)),CrvEll) or
             (ISA(Type(Codomain(MWmap)),SetPtEll) and
              Ring(Codomain(MWmap)) cmpeq BaseRing(Scheme(Codomain(MWmap))))):
    "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
    E:=Codomain(MWmap);
    if ISA(Type(Codomain(MWmap)),SetPtEll) then
      E:=Scheme(E);
    end if;
    Ep,toEp:=Reduction(E,P);
    k,tok:=ResidueClassField(P);
    EpGrp,EpGrpMap:=AbelianGroup(Ep);
    InvEpGrpMap:=map<Codomain(EpGrpMap)->Domain(EpGrpMap)|
                                 g:->DiscreteLog(EpGrpMap,g)>;
    rc:=GetLocalData(MWmap,P);
    rc`reduction:=<
        hom<Domain(MWmap)->EpGrp|
         [g@MWmap@toEp@InvEpGrpMap:g in OrderedGenerators(Domain(MWmap))]>,
        EpGrpMap>;
    StoreLocalData(MWmap,P,rc);
  end if;
  rc:=GetLocalData(MWmap,P);
  return rc`reduction[1],rc`reduction[2];
end intrinsic;

intrinsic IsPSaturated(mwmap::Map,p::RngIntElt)->BoolElt
  {Tests if given MWmap is p-saturated}
  G:=Domain(mwmap);
  pG:=sub<G|[p*g:g in OrderedGenerators(G)]>;
  GModpG,ModMap:=quo<G|pG>;
  return forall{g:g in GModpG| g eq Zero(GModpG) or
    #DivisionPoints(mwmap(g@@ModMap),p) eq 0};
end intrinsic;

intrinsic FormalLog(E::CrvEll:Precision:=10)->RngSerLaurElt,PtEll
  {Returns a formal power series log(z) of E, together with a point P(z),
  such that the map P(z):->log(z) describes the formal logarithm.}
  R:=BaseRing(E);
  a1,a2,a3,a4,a6:=Explode(aInvariants(E));
  Rz<z>:=LaurentSeriesRing(R:Precision:=Precision);
  Rzw<w>:=PolynomialRing(Rz);
  rt:=Roots(w-(z^3+a1*z*w+a2*z^2*w+a3*w^2+a4*z*w^2+a6*w^3));
  rt:=[r:r in rt | Valuation(r[1]) eq 3];
  assert #rt eq 1 and rt[1][2] eq 1;
  Wz:=rt[1][1];
  assert Valuation(Wz) eq 3;
  omega:=(Wz-z*Derivative(Wz))/Wz/(a1*z+a3*Wz-2);
  log:=Integral(omega);
  Pz:=E(Rz)![z/Wz,-1/Wz];
  return log,Pz;
end intrinsic;

intrinsic RelevantCosets(MWmap::Map,Ecov::MapSch,Prs::{RngOrdIdl})->Tup
  {Computes useful information for Chabauty methods using elliptic covers}
  MWgrp:=Domain(MWmap);
  E:=Domain(Ecov);
  require Codomain(MWmap) eq E:"MWmap should correspond to cover";
  require ISA(Type(Codomain(Ecov)),Prj) and Dimension(Codomain(Ecov)) eq 1:
    "cover should be to a projective line";
  require BaseRing(Codomain(Ecov)) eq Rationals():
    "The cover should map to a projective line defined over the rationals.";
  require forall{P:P in Prs|IsPrime(P)}: "Ideals should all be prime";
  L:={Minimum(P):P in Prs};
  require #L eq 1: "Primes should all be of same residue characteristic";
  p:=Rep(L);
  Prs:=[LookupPrime(P):P in Prs];
  
  Ep,toEp:=Explode(<[v[1]:v in V],[v[2]:v in V]>) where
               V:=[<a,b>where a,b:=Reduction(E,P):P in Prs];
  k,tok:=Explode(<[v[1]:v in V],[v[2]:v in V]>) where
               V:=[<a,b>where a,b:=ResidueClassField(p):p in Prs];
  P1p:=ProjectiveSpace(GF(p),1);
  Ecovp:=[map<Ep[i]->P1p|ProjectiveReduction(
     DefiningEquations(Ecov),Prs[i],CoordinateRing(Ambient(Ep[i])))>:i in
     [1..#Prs]];
  MWred,EpGrpMap:=Explode(<[v[1]:v in V],[v[2]:v in V]>) where
                        V:=[<a,b>where a,b:=Reduction(MWmap,P):P in Prs];
  IDX:=1;
  fibs:=[* *];
  for i in [1..#Prs] do
    P1k:=P1p(k[i]);
    V:={P1k!Eltseq(pt):pt in RationalPoints(P1p)};
    IDX:=LCM(IDX,#(Codomain(MWred[i])/Image(MWred[i])));
    V:={<p,P1p!Eltseq(v)>:p in Image(MWred[i])| v in V where
         v:=EvaluateByPowerSeries(Ecovp[i],EpGrpMap[i](p))};
    Vim:={v[2]:v in V};
    fibs[i]:={<v,{p[1]:p in V|p[2] eq v}>:v in Vim};
  end for;

  fib1:=fibs[1];
  m1:=MWred[1];
  G:=FreeAbelianGroup(Ngens(MWgrp));
  m1:=hom<G->Image(m1)|[m1(g):g in OrderedGenerators(Domain(m1))]>;

  for i in [2..#Prs] do
    fib2:=fibs[i];
    m2:=MWred[i];
    m2:=hom<G->Image(m2)|[m2(g):g in OrderedGenerators(Domain(m2))]>;

    gcd,togcd:=quo<G|Generators(Kernel(m1))join Generators(Kernel(m2))>;
    lcm,tolcm:=quo<G|Kernel(m1) meet Kernel(m2)>;
    lcmto1:=hom<lcm->Codomain(m1)|[m1(a@@tolcm):a in OrderedGenerators(lcm)]>;
    lcmto2:=hom<lcm->Codomain(m2)|[m2(a@@tolcm):a in OrderedGenerators(lcm)]>;
    m1togcd:=hom<Codomain(m1)->gcd|
        [togcd(a@@m1):a in OrderedGenerators(Codomain(m1))]>;
    m2togcd:=hom<Codomain(m2)->gcd|
        [togcd(a@@m2):a in OrderedGenerators(Codomain(m2))]>;
    function CRT(e1,e2)
      if m1togcd(e1) eq m2togcd(e2) then
        v:=Eltseq(e1@@m1-e2@@m2);
        A1:=[Eltseq(G!a):a in OrderedGenerators(Kernel(m1))];
        A2:=[Eltseq(G!a):a in OrderedGenerators(Kernel(m2))];
        w:=Solution(Matrix(A1 cat A2),Vector(v));
        crt:=tolcm(e1@@m1-Kernel(m1)![w[i]:i in [1..Ngens(Kernel(m1))]]);
        return {crt+a:a in Kernel(lcmto1)meet Kernel(lcmto2)};
      else
        return {lcm|};
      end if;
    end function;

    fib1:={ <f1[1],W>: f1 in fib1,f2 in fib2| f1[1] eq f2[1] and W ne {}
                      where W:=&join[CRT(e1,e2):e1 in f1[2],e2 in f2[2]]};
    m1:=tolcm;
  end for;
  return <hom<MWgrp->Codomain(m1)|[m1(g):g in OrderedGenerators(G)]>,
                                                   &join[f[2]:f in fib1]>,IDX;
end intrinsic;

intrinsic CosetIntersection(V::Tup,W::Tup:Weak:=false)->Tup
  {Computes the intersection of two collections of cosets of an abelian group.
   If the parameter Weak is supplied, returns cosets of the kernel associated
   with V, otherwise returns cosets of the intersection of the kernels
   associated with V and W}
  require #V eq 2 and #W eq 2 and ISA(Type(V[1]),Map) and ISA(Type(W[1]),Map)
  and ISA(Type(V[2]),SetEnum) and ISA(Type(W[2]),SetEnum) and
  ISA(Type(Domain(V[1])),GrpAb) and Domain(V[1]) eq Domain(W[1]) and
  Universe(V[2]) eq Codomain(V[1]) and Universe(W[2]) eq Codomain(W[1]):
    "V and W must describe coset collections of the same abelian group";
  m1:=V[1];m2:=W[1];
  if Weak then
    G:=Domain(m1);
    gcd,togcd:=quo<G|Generators(Kernel(m1))join Generators(Kernel(m2))>;
    m1togcd:=hom<Codomain(m1)->gcd|
        [togcd(a@@m1):a in OrderedGenerators(Codomain(m1))]>;
    m2togcd:=hom<Codomain(m2)->gcd|
        [togcd(a@@m2):a in OrderedGenerators(Codomain(m2))]>;
    gcdW:={m2togcd(w):w in W[2]};
    return <V[1],{Universe(V[2])|v:v in V[2]| m1togcd(v) in gcdW}>;
  else
    G:=FreeAbelianGroup(Ngens(Domain(V[1])));
    m1:=hom<G->Image(m1)|[m1(g):g in OrderedGenerators(Domain(m1))]>;
    m2:=hom<G->Image(m2)|[m2(g):g in OrderedGenerators(Domain(m2))]>;
    gcd,togcd:=quo<G|Generators(Kernel(m1))join Generators(Kernel(m2))>;
    lcm,tolcm:=quo<G|Kernel(m1) meet Kernel(m2)>;
    lcmto1:=hom<lcm->Codomain(m1)|[m1(a@@tolcm):a in OrderedGenerators(lcm)]>;
    lcmto2:=hom<lcm->Codomain(m2)|[m2(a@@tolcm):a in OrderedGenerators(lcm)]>;
    m1togcd:=hom<Codomain(m1)->gcd|
        [togcd(a@@m1):a in OrderedGenerators(Codomain(m1))]>;
    m2togcd:=hom<Codomain(m2)->gcd|
        [togcd(a@@m2):a in OrderedGenerators(Codomain(m2))]>;
    function CRT(e1,e2)
      if m1togcd(e1) eq m2togcd(e2) then
        v:=Eltseq(e1@@m1-e2@@m2);
        A1:=[Eltseq(G!a):a in OrderedGenerators(Kernel(m1))];
        A2:=[Eltseq(G!a):a in OrderedGenerators(Kernel(m2))];
        w:=Solution(Matrix(A1 cat A2),Vector(v));
        crt:=tolcm(e1@@m1-Kernel(m1)![w[i]:i in [1..Ngens(Kernel(m1))]]);
        return {crt+a:a in Kernel(lcmto1)meet Kernel(lcmto2)};
      else
        return {lcm|};
      end if;
    end function;
    //This creates #V[2] * #W[2] complexity. Especially if
    //#gcd is comparatively large, one can do quite a bit better by first
    //sorting W[2] according to class in gcd and only CRT-ing with
    //elements of V[2] of the same class in gcd.
    return <hom<Domain(V[1])->lcm|[tolcm(g):g in OrderedGenerators(G)]>,
                  &join[CRT(v,w):v in V[2],w in W[2]]>;
  end if;
end intrinsic;

function CoeffList(f)
  Lp:=Parent(f);
  if ISA(Type(Lp),FldPad) then
    return Eltseq(f);
  end if;
  OLp:=IntegerRing(Lp);
  Qp:=pAdicField(Lp);
  if Valuation(f) ge 0 then
    return [Qp!c:c in &cat(Locseq(OLp!f))];
  else
    p:=UniformizingElement(Qp);
    v:=Valuation(f) div Valuation(p);
    return [c*p^v:c in &cat(Locseq(OLp!(f/p^v)))];
  end if;
end function;  

intrinsic ChabautyEquations(P0::PtEll,Ecov::MapSch,
                    MWmap::Map,Prs::{RngOrdIdl}:Precision:=10)->SetEnum
  {Returns a set of equations used in Chabauty methods}

  E:=Codomain(MWmap);
  require Codomain(MWmap) eq E:"MWmap should correspond to cover";
  require ISA(Type(Codomain(Ecov)),Prj) and Dimension(Codomain(Ecov)) eq 1:
    "cover should be to a projective line";
  require BaseRing(Codomain(Ecov)) eq Rationals():
    "The cover should map to a projective line defined over the rationals.";
  require forall{P:P in Prs|IsPrime(P)}: "Ideals should all be prime";
  L:={Minimum(P):P in Prs};
  require #L eq 1: "Primes should all be of same residue characteristic";
  require P0 in E: "Point should be on elliptic curve";

  Elog,Pz:=FormalLog(Domain(Ecov):Precision:=Precision-1);
  assert AbsolutePrecision(Elog) ge Precision;

  Prs:=[LookupPrime(P):P in Prs];
  
  covP0z:=a[1]/a[2] where a:=Ecov(P0+Pz);
  Kz:=Parent(covP0z);
  if Valuation(covP0z) lt 0 then covP0z:=1/covP0z; end if;

  covP0:=Evaluate(covP0z,0);
  if not(covP0 in Rationals()) then
    error "P0 has no rational image";
  end if;

  vals:=[IntNInf()|Valuation(covP0,p):p in Prs];
  
  if Rep(vals) lt 0 then
    covP0z:=1/covP0z;
  end if;
  
  vprint EllChab:"Power series expansion:",covP0z;
  covP0z:=covP0z-Evaluate(covP0z,0);
  covP0z:=covP0z/UniformizingElement(Parent(covP0z))^(Valuation(covP0z)-1);
  vprint EllChab:"After translation and renormalization:",covP0z;
  
  //assert AbsolutePrecision(covP0z) ge Precision;
  
  KerPrs:=Domain(MWmap);
  for P in Prs do
    KerPrs meet:=Kernel(Reduction(MWmap,P));
  end for;

  KerPrsVec:=[Eltseq(Domain(MWmap)!g):g in OrderedGenerators(KerPrs)];

  Qp:=pAdicField(Minimum(Prs[1]));
  Qpn:=LazyPowerSeriesRing(Qp,Ngens(KerPrs));

  Eqs:=[];
  zBmat:=[];

  for p in Prs do
    Kp,toKp:=MyCompletion(p);
    error if Valuation(Kp!Minimum(p))/(Minimum(p)-1) ge 1,
      "Prime has too high ramification compared to residue characteristic.";
    
    pi:=UniformizingElement(Kp);
    EKp:=PointSet(E,toKp);
    prc:=Precision*RamificationDegree(Kp);
    repeat
      error if prc gt Kp`DefaultPrecision, "Insufficient precision available";
      Gp:=[EKp![toKp(c)+O(pi^prc):c in Eltseq(MWmap(g))]:
                   g in OrderedGenerators(Domain(MWmap))];
      zBp:=[-p[1]/p[2] where p:=&+[v[i]*Gp[i]:i in [1..#Gp]]:v in KerPrsVec];
      prc:=prc*2;
    until MinPrec(zBp) ge Precision*RamificationDegree(Kp);

    zBmat:=zBmat cat [[Integers()!(CoeffList(a)[i]/Minimum(p)): a in zBp]:
                    i in [1..Degree(Kp)]];

    Kpz:=LaurentSeriesRing(Kp:Precision:=Kz`Precision);
    toKpz:=map<Kz->Kpz|f:->elt<Kpz|Valuation(f),[toKp(a):a in Eltseq(f)]>>;
    Elogp:=toKpz(Elog);
    Eexpp:=Reverse(Elogp);
    lBp:=[Evaluate(Elogp,b)+O(pi^(Valuation(b)*AbsolutePrecision(Elogp))):b in zBp];
    Kpn:=LazyPowerSeriesRing(Kp,#lBp);
    lBpn:=&+[lBp[i]*Kpn.i:i in [1..#lBp]];

    covP0pz:=toKpz(covP0z);
    zBpn:=Evaluate(Eexpp,lBpn);
    covBpn:=Evaluate(covP0pz,zBpn);
    RcovBpn:=[elt<Qpn|map<CartesianPower(Integers(),Rank(Qpn))->Qp|
                         v:->(CoeffList(Coefficient(covBpn,[c:c in v])))[j]>>:
            j in [1..Degree(Kp)]];
    error if MinPrec(&cat[Coefficients(S,Precision):S in RcovBpn]) lt
        Precision-2, "Insufficient precision computed";
    Append(~Eqs,RcovBpn);
  end for;
  vprint EllChab:"Reduction of matrix of Z-coordinates:",
                     Matrix(GF(Minimum(Prs[1])),zBmat);
  Eq0:=Eqs[1][1];
  Eqs:=&cat[[i eq 1 select l[i]-Eq0 else l[i]:i in [1..#l]]:l in Eqs];
  if Rank(Matrix(GF(Minimum(Prs[1])),zBmat)) lt Ngens(KerPrs) then
    vprint EllChab:
      "Matrix does not have maximal rank. Provide separate proof that MW-group",
      "is saturated at",Minimum(Prs[1]),"if required";
    return Eqs[2..#Eqs],Minimum(Prs[1]);
  else
    vprint EllChab:
      "Matrix has full rank. Implicit proof that MW-group is saturated at",
      Minimum(Prs[1]);
    return Eqs[2..#Eqs],1;
  end if;
end intrinsic;

intrinsic IsIntegral(E::CrvEll,p::RngOrdIdl) ->BoolElt
  {returns true if the defining coefficients of E have non-negative valuation at
  the prime ideal p}
  cfs:=Eltseq(E);
  return forall{c:c in cfs|Valuation(c,p) ge 0};
end intrinsic;

intrinsic Chabauty(MWmap::Map, Ecov::MapSch,
            p::RngIntElt:Aux:={},Precision:=2,Bound:=5)->SetEnum,RngIntElt
  {Straightforward attempt for determining the set of points in the image
  of MWmap that have a rational image under Ecov}
  K:=BaseRing(Domain(Ecov));
  OK:=IntegerRing(K);
  Prs:=Support(p*OK);
  vprint EllChab:"Computing relevant cosets using the primes above",p;
  fibp,IDX:=RelevantCosets(MWmap,Ecov,Prs);
  vprint EllChab:"Found cosets",fibp[2];
  vprint EllChab:"Modulo kernel",Kernel(fibp[1]);
  vprint EllChab:"LCM of indices:",IDX;
  for q in Aux do
    vprint EllChab:"Computing relevant cosets using the primes above",q;
    fibq,IDXq:=RelevantCosets(MWmap,Ecov,Support(q*OK));
    vprint EllChab:"Found cosets:",fibq[2];
    vprint EllChab:"Modulo kernel:",Kernel(fibq[1]);
    vprint EllChab:"LCM of indices:",IDXq;
    fibp:=CosetIntersection(fibp,fibq:Weak);
    vprint EllChab:"For",p,"this cuts down cosets to",fibp[2];
    IDX:=LCM(IDX,IDXq);
  end for;
  vprint EllChab:
    "Testing linear combinations in Mordell-Weil group up to bound",Bound;
  V:={v: c in CartesianPower([-Bound..Bound],Ngens(Domain(MWmap))) |
         fibp[1](v) in fibp[2] where v:=Domain(MWmap)![d:d in c]};
  QQ:=Rationals();
  MWfib:={a: a in V| b[1] in QQ and b[2] in QQ where b:=
              EvaluateByPowerSeries(Ecov,MWmap(a))};
  vprint EllChab:"Found",MWfib;

  error if {fibp[1](pt):pt in MWfib} ne fibp[2],
    "Cannot find representatives for all cosets.",
    "Increase Bound or give more auxiliary primes";
  Qp:=pAdicField(p);
  GFp,toGFp:=ResidueClassField(IntegerRing(Qp));
  vprint EllChab:
    "Computing Chabauty-bound for each point using primes above",p;
  for g in MWfib do
    P0:=MWmap(g);
    vprint EllChab:"P0 =",P0;
    vprint EllChab:"Image under cover:",EvaluateByPowerSeries(Ecov,P0);
    Eqs,IDXq:=ChabautyEquations(P0,Ecov,MWmap,Prs);
    IDX:=LCM(IDX,IDXq);
    if MinValuation([Coefficients(l,0)[1]:l in Eqs]) lt 1 then
      error "Finding Chabauty bound 0 in a coset where we know a point";
    end if;
    M:=Matrix([[toGFp(c/Prime(Qp)):c in Coefficients(l,1,1)]:l in Eqs]);
    vprint EllChab:"Linear approximation gives matrix:",M;
    if Rank(M) lt NumberOfColumns(M) then
      error "Finding a Chabauty bound of more than 1";
    end if;
    vprint EllChab:"Chabauty bound for this point is 1";
  end for;
  vprint EllChab:
    "Returned set gives all points with rational image if given MWmap is",IDX,
    "saturated.";
  return MWfib,IDX;
end intrinsic;


/****************************************************************************
Some other routines:
****************************************************************************/

intrinsic SwapExtension(R::Rng)->Rng,Map
  {takes a polynomial ring over a number field and swaps the order of the
  two topmost extensions}
  if ISA(Type(R),RngUPol) then
    K:=BaseRing(R);
    assert ISA(Type(K),FldNum);
    k:=BaseField(K);
    Pk:=PolynomialRing(k);
    PkX:=PolynomialRing(Pk);
    PkK:=quo<PkX|Evaluate(DefiningPolynomial(K),PkX.1)>;

    KtoPkK:=hom<K->PkK|Coercion(k,PkK),PkK.1>;
    RtoPkK:=hom<R->PkK|KtoPkK, PkK!Pk.1>;

    PktoR:=hom<Pk->R|Coercion(k,R), R.1>;
    PkKtoR:=hom<PkK->R|PktoR, R!(K.1)>;

    return PkK,map<R->PkK|a:->RtoPkK(a),b:->PkKtoR(b)>;
  elif ISA(Type(R),RngMPol) then
    K:=BaseRing(R);
    if ISA(Type(K),FldNum) then
      k:=BaseField(K);
      modulus:=DefiningPolynomial(K);
    elif ISA(Type(K),RngUPolRes) then
      k:=BaseRing(K);
      modulus:=Modulus(K);
    else
      error "Can't do that type yet";
    end if;
    Pk:=PolynomialRing(k,Rank(R));
    PkX:=PolynomialRing(Pk);
    PkK:=quo<PkX|Evaluate(modulus,PkX.1)>;

    KtoPkK:=hom<K->PkK|Coercion(k,PkK),PkK.1>;
    RtoPkK:=hom<R->PkK|KtoPkK, [PkK!Pk.i: i in [1..Rank(R)]]>;

    PktoR:=hom<Pk->R|Coercion(k,R), [R.i: i in [1..Rank(R)]]>;
    PkKtoR:=hom<PkK->R|PktoR, R!(K.1)>;

    return PkK,map<R->PkK|a:->RtoPkK(a),b:->PkKtoR(b)>;
  elif ISA(Type(R),RngUPolRes) then
    K:=BaseRing(R);
    assert ISA(Type(K),FldNum);
    k:=BaseField(K);
    fR:=Modulus(R);
    Pk:=PolynomialRing(k);
    fk:=Polynomial(k,fR);//the extension must be descendable ...
    Ak:=quo<Pk|fk>;
    AkX:=PolynomialRing(Ak);
    AK:=quo<AkX|Evaluate(DefiningPolynomial(K),AkX.1)>;
    
    KtoAK:=hom<K->AK|Coercion(k,AK),AK.1>;
    RtoAK:=hom<R->AK|KtoAK,AK!(Ak.1)>;
    
    AktoR:=hom<Ak->R|Coercion(k,R), R.1>;
    AKtoR:=hom<AK->R|AktoR, R!(K.1)>;
    
    return AK,map<R->AK|a:->RtoAK(a),b:->AKtoR(b)>;
  else
    error "not implemented";
  end if;  
end intrinsic;

intrinsic PowerFreePart(n::FldRatElt,e::RngIntElt)->RngIntElt,RngIntElt
{returns the e-th power free part of n and the e-th root of the remaining part
(local intrinsic by Nils)}
  local R,L,u;
  if n eq 0 then
    return 1,0;
  end if;
  L,u:=Factorisation(n);
  if Abs(n) eq 1 then
    R:=[1,1];
  else
    R :=
    [&*{@ u[1]^(u[2] mod e) : u in L @},&*{@ u[1]^(u[2] div e) : u in L @}];
  end if;
  if e mod 2 eq 0 then
    return u*R[1],R[2];
  else
    return R[1],u*R[2];
  end if;
end intrinsic;

intrinsic Factorisation(n::FldRatElt) -> SeqEnum
{Factorisation of numerator and denominator. By abuse of notation, this is
"factorisation" of a rational number (local intrinsic by Nils)}
  L,u:=Factorisation(Numerator(n));
  L:= L cat [ <u[1],-u[2]> : u in Factorisation(Denominator(n))];
  return L,u;
end intrinsic;
