// import "/home/magma/package/Geometry/Algaem/localdata.m":GetLocalData,StoreLocalData;
import "localdata.m":GetLocalData,StoreLocalData;
// import "/home/magma/package/Geometry/Arith/loctools.m":
//        Minim, Maxim, IntNInf, MinValuation, MinPrec, ShiftVal,
//       StripPrecisionlessTerms, FlattenPrecision, pAdicEvaluate;
       
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 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);
  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;  
