// intrinsic not present

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;

//problem with projective spaces of dimension 0

intrinsic TestEquations(Eqs::[RngPowLazElt],d::RngIntElt)->BoolElt
  {Test equations}
  R:=Universe(Eqs);
  r:=Rank(R);
  Qp:=BaseRing(R);
  p:=Prime(Qp);
  dpad:=d-FactorialValuation(d,p);
  v:=Minimum( [Valuation(c):c in &cat[Coefficients(f,d,d):f in Eqs]]);
  d2:=v+1;
  while d2-FactorialValuation(d2,p) le v do
    v2:=Minimum( [Valuation(c):c in &cat[Coefficients(f,d2,d2):f in Eqs]]);
    if v2 le v then
       vprint EllChab:"Homogeneous part not decisive for this equation";
       return false;
    end if;
    d2+:=1;
  end while;
  Fp,toFp:=ResidueClassField(IntegerRing(Qp));
  if d eq 1 then
    M:=Matrix([[toFp(c/p^v):c in Coefficients(l,1,1)]:l in Eqs]);
    vprint EllChab:"Linear approximation gives matrix:",M;
    return Rank(M) eq NumberOfColumns(M);
  else
    P:=Proj(PolynomialRing(Fp,r));
    V:=MonomialsOfDegree(CoordinateRing(P),d);
    L:=[&+[toFp(Coefficient(f,Exponents(m))/p^v)*m:m in V]:f in Eqs];
    X:=Scheme(P,L);
    W:=RationalPoints(X);
    return #W eq 0;
  end if;
end intrinsic;

//NormEquation fix

intrinsic MyHasPoint(X::CrvCon)->BoolElt,Pt
  {}
  K:=BaseRing(X);
  fd,M:=DiagonalForm(X);
  fd:=LCM([Denominator(c):c in Coefficients(fd)])*fd;
  A:=K!Coefficient(fd,1,2);
  B:=K!Coefficient(fd,2,2);
  C:=K!Coefficient(fd,3,2);
  dn:=Denominator(A)*Denominator(B)*Denominator(C);
  g:=Polynomial(K,[A*B,0,1]);
  bl,rt:=HasRoot(g);
  if bl then
    sol:=[rt,A,0];
  else
    L:=NumberField(Polynomial(K,g));
    bl,sol:=NormEquation(L,-C*A:Nice);
    if bl then
      sol:=Eltseq(sol[1]) cat [1];
      sol:=[sol[1],-A*sol[2],A*sol[3]];
    end if;
  end if;
  if bl then
    return true,X!Eltseq(Vector(sol)*M);
  else
    return false,_;
  end if;
end intrinsic;

//Seems to be missing

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;

///Included functions to support Chabauty

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
    cf:=[c div pi^(-v):c in cf];
  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;

intrinsic MonomialCoefficient(f::RngMPolElt,v::[RngIntElt])->RngElt
  {Returns the coefficient in f of the monomial defined by the exponent vector
  v}
  require #v eq Rank(Parent(f)):"Vector must be of length equal to the rank of the polynomial ring";
  return CoefficientByExponentVector(f,v);
end intrinsic;

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;

function toint(a)
  Zp:=Parent(a);
  p:=Prime(Zp);
  N:=p^AbsolutePrecision(a);
  an:=Integers()!a;
  if Abs(an-N) lt Abs(an) then
    return an-N;
  else
    return an;
  end if;
end function;

///end of included functions to support Chabauty

//Chabauty included because it could not find CosetIntersection

intrinsic Chabauty(MWmap::Map, Ecov::MapSch,
            p::RngIntElt:Cosets:=false,Aux:={},Precision:=10,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;
  if Cosets cmpne false then
    require #Aux eq 0:"Either specify Cosets OR Aux, but not both.";
    fibp:=CosetIntersection(fibp,Cosets:Weak);
    vprint EllChab:
      "Insersecting with supplied cosets cuts down cosets to",fibp[2];
  else
    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;
  end if;
  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;

  if {fibp[1](pt):pt in MWfib} ne fibp[2] then
    vprint EllChab:
    "Cannot find representatives for all cosets. We'll see what happens.";
  end if;
  rem:=fibp[2];
  kerlev:=Infinity();
  nogos:=[];
  NrPoints:=0;
  for g in MWfib do
    P0:=MWmap(g);
    vprint EllChab:"P0 =",P0;
    vprint EllChab:"Image under cover:",EvaluateByPowerSeries(Ecov,P0);
    Eqs,order,IDXq:=ChabautyEquations(P0,Ecov,MWmap,Prs:
                                  Centred,Precision:=Precision);
    vprint EllChab:"Order of vanishing:",order;
    IDX:=LCM(IDX,IDXq);
    bl:=TestEquations(Eqs,order);
    if bl then
      vprint EllChab:"P0 is the only point in its fibre.";
      NrPoints+:=1;
      rem:=rem diff {fibp[1](g)};
    else
      vprint EllChab:"P0 is not necessarily unique in its fibre.";
      if order gt 1 then      
        vprint EllChab:
          "P0 is a ramification point.",
          "We do not get a bound on the number of rational points";
      end if;      
    end if;
  end for;
  if #rem eq 0 then
    vprint EllChab:"Done!!";
  else
    vprint EllChab:"Will try to deal with remaining points";
  end if;
  candidates:=[];
  for g in rem do
    bl:=exists(gr){gr:gr in MWfib| fibp[1](gr) eq g};
    if not bl then
      gr:=ShortLift(fibp[1],g);
    end if;
    P0:=MWmap(gr);
    vprint EllChab:"P0 =",P0;
    Eqs,order,IDXq:=ChabautyEquations(P0,Ecov,MWmap,Prs:Precision:=Precision);
    IDX:=LCM(IDX,IDXq);
    V1,V2:=SolveEquations(Eqs);
    vprint EllChab:"Solutions to p-adic system of equations:",V1;
    vprint EllChab:"Unresolved neighbourhoods:",V2;
    NrPoints+:=#V1;
    if #V2 ne 0 then
      NrPoints:=Infinity();
    end if;
    kerlev:=Minimum(kerlev,MinPrec(&cat Setseq(V1)));
    kerlev:=Minimum(kerlev,MinPrec(&cat Setseq(V2)));
    for v in V1 do
      Append(~candidates,gr+Kernel(fibp[1])![toint(c):c in v]);
    end for;
    for v in V2 do
      Append(~nogos,gr+Kernel(fibp[1])![toint(c):c in v]);
    end for;
  end for;
  for a in candidates do
    if forall{c: c in Eltseq(a)| Abs(c) lt 20} then
      vprint EllChab:"Point",a,"looks promisingly small. Let's try it.";
      b:=EvaluateByPowerSeries(Ecov,MWmap(a));
      if b[1] in QQ and b[2] in QQ then
        vprint EllChab:"Found point",a,"with rational image",b;
        MWfib join:={a};
      else
        vprint EllChab:"Didn't work out";
        Append(~nogos,a);
      end if;
    else
      vprint EllChab:"Point",a,"is supposedly a sole candidate,",
        "but looks too big to really try.";
      Append(~nogos,a);
    end if;
  end for;
  if #nogos eq 0 then
    vprint EllChab:
      "Returned set gives all points with rational image if given MWmap is",IDX,
      "saturated.";
    newtup:=<fibp[1],{Universe(fibp[2])|}>;
  else  
    vprint EllChab:
      "There are unresolved neighbourhoods.";
    _,toquo:=quo<Domain(MWmap)|
     [p^kerlev*g:g in OrderedGenerators(Kernel(fibp[1]))]>;
    newtup:=<toquo,{toquo(a):a in nogos}>;
  end if;
  return NrPoints,MWfib,IDX,newtup;
end intrinsic;

//part of HasRatXPoint

function HasSquareIntegralRatValues(Fs : Neighbourhood:=false)

/*Tests for integral solutions to y^2=F(x). The second return value
is an approximation of the x of such a point, if it exists. If the optional
parameter Neighbourhood is given, then only x that match up with Neighbourhood
(given the precision of Neighbourhood) are tried. If the flag NoRoots is
supplied, then F is not tested for having a p-adic root. Since p-adic
factorisation is slow, this saves some time, but may lead to very poor
performance and failure (loss of precision) if there is such a root*/
  
  vprint LocSol,2:"Looking for x in Zp giving square values in all of",Fs;
      
  // now test the conditions 
  
  Op:=BaseRing(Universe(Fs));
  p:=Prime(Op);
  Zp:=pAdicRing(p);
  pi:=UniformizingElement(Op);
  error if not(ISA(Type(Op),RngPad)), "Polynomial needs to be over a local ring";
  if Neighbourhood cmpeq false then
    Neighbourhood:=O(Op!p^0);
  end if;
  error if not( Parent(Neighbourhood) cmpeq Op),
    "Neighbourhood (if given) must be a member of the local ring";
  resfld,toresfld:=ResidueClassField(Op);
  resZp,toresZp:=ResidueClassField(Zp);
  resfldreps:={ u+O(Op!p): u in [0..p-1] };
//  print resfldreps;
  if (Prime(Op) ne 2) then
    // this is the test-function for odd primes
    // this function takes an element (an origin) and
    // recursively looks in smaller neigbourhoods around it (while
    // perturbing it within its given precision)
    // to decide if F takes square values somewhere near origin.
    // if true is returned, then the second value gives a certificate
    // (i.e. an x for which F takes square values for all x+d within
    // the given precision of x)
    u:=elt<Op|UniformizingElement(Op),2>;
    test:=function(origin)
    
        // first, if we've really run out of precision, give an error
	
	if (Precision(Op) le Precision(origin)) then
          error "Insufficient precision for decision";
	end if;
	prc:=Precision(origin);

	//increase precision of origin in order to represent pertubations
	origin:=elt<Op|origin,prc+Valuation(Op!p)>;
	vprint LocSol,2:"Looking for point in",origin;
	// loop through pertubations	
	for perturb in resfldreps do
	  assert prc mod Valuation(Op!p) eq 0;
          x:=origin+p^(prc div Valuation(Op!p))*perturb+
	                           O(u^(prc+Valuation(Op!p)));
//	  print "perturb=",perturb,"trying",x;
          bl:=1;
	  for F in Fs do
	    y:=Evaluate(F,x);
	    vly:=Valuation(y);
	    if vly ge AbsolutePrecision(y) then
	      bl:=0;
	    elif ((vly mod 2) eq 0) and IsSquare(resfld!(y/u^vly)) then
            else
	      bl:=-1;
	      break;
	    end if;
	  end for;
          if bl eq 1 then
	    return true,x;
	  elif bl eq 0 then
	    pred,val:=$$(x);
	    if pred then
	      return true,val;
	    end if;
	  end if;
	end for;
	  
	//by now, we've tested everything and haven't found any square value
	//therefore, there are none.
        vprint LocSol,2:"Found none.";
	return false,_;
      end function;
  else
  
    // we construct the square units. That is, we take all units with precision
    // Valuation(2)+1 and square them. This gives us precision Valuation(4)+1
    // which is enough precision to tell square units from non-square ones.
    sqr:={
      (&+[pi^(i-1)*(u[i]@@toresfld)+O(pi^(#u)) : i in [1..#u]])^2
      +O(pi^(Valuation(Op!4)+1))  :
           u in Flat(car<{u : u in resfld | u ne 0},
                      CartesianPower(resfld,Valuation(Op!2))>)};
    v4p1:=Valuation(Op!4)+1;
    u:=pi+O(pi^(v4p1+1));
//    print "<pi,v4p1,pi+O(pi^(v4p1+1))>:",<pi,v4p1,pi+O(pi^(v4p1+1))>;
    fs:=[Eltseq(F): F in Fs];
    test:=function(origin)
    
        // first, if we've really run out of precision, give an error
	// It's not clear how this generalizes to the new locals ...
	//if (Precision(Op) lt AbsolutePrecision(origin)+v4p1) then
        //  error "Insufficient precision for decision";
	//end if;
	prc:=AbsolutePrecision(origin);

	//increase precision of origin in order to represent pertubations
	origin:=elt<Op|origin,prc+Valuation(Op!p)>;
	vprint LocSol,2:"Looking for point in",origin;
	// loop through pertubations	
	for perturb in resfldreps do
          x:=origin+p^(prc div Valuation(Op!p))*perturb+O(u^(prc+Valuation(Op!p)));
          bl:=1;
	  for f in fs do
	    y:=&+[f[i]*(x^(i-1)):i in [1..#f]];
	    vly:=Valuation(y);

	    // is there *ANY* information in y?
	    if RelativePrecision(y) gt 0 then

	      // if y has odd valuation, there's no points around this x.
	      if (vly mod 2) eq 0 then

		// otherwise, multiply y into a unit (by a square!)
		// and throw away any possibly excess precision.
//        	print "y,u,vly,u^vly",<y,u,vly,u^vly>;
		yrd:=y/u^vly+elt<Op|v4p1,0,v4p1>;

		// now see if yrd is a square up to its precision
		// (because we've thrown out excess precision, ALL digits
		// of yrd should correspond to a listed square in sqr)

		if exists{sq: sq in sqr|
                             RelativePrecision(yrd-sq) eq 0} then
	          // if yrd had full precision, then we're really done.
	          if AbsolutePrecision(yrd) ge v4p1 then

	          // Otherwise, look closer (i.e. recurse around this x)
	          else
		    bl:=0;
		  end if;
		else
		  bl:=-1;
		  break;
        	end if;
              else
	        bl:=-1;
		break;
	      end if;

	    // if there was no info in y, look closer (recurse around this x).
	    else
	      bl:=0;
            end if;
	  end for;
	  if bl eq 1 then
	    return true,x;
	  elif bl eq 0 then
	    pred,val:=$$(x);
	    if pred then
	      return true,val;
	    end if;
	  end if;
	end for;

 	//by now, we've tested everything and haven't found any square value
	//therefore, there are none.
        vprint LocSol,2:"Found none.";
	return false,_;
      end function;
   
  end if;
  vprint LocSol,2:"Looking for point in",Neighbourhood;
  IndentPush();
  a,b:=test(Neighbourhood);
  IndentPop();
  if a then
    vprint LocSol,2:"Found",b;
    return a,b;
  else
    vprint LocSol,2:"Found none.";
    return a,_;
  end if;
end function;

//missing intrinsic

intrinsic HasRatXpoint(Fsin::[RngUPolElt])->SeqEnum
  {takes a list of polynomials and decides if there is a locally rational
  value x where all polynomials take a square value simultaneously}
  KpX:=Universe(Fsin);
  Kp:=BaseRing(KpX);
  Op:=IntegerRing(Kp);
  pi:=UniformizingElement(Op);
  OpX:=PolynomialRing(Op);
  p:=Prime(Kp);
  Qp:=pAdicField(p);
  Zp:=IntegerRing(Qp);
  Fint:=[];
  for f in Fsin do
    v:=2*Floor(Minimum([Valuation(c):c in Eltseq(f)])/2);
    Append(~Fint,OpX![Op!(c/(Kp!pi)^v):c in Eltseq(f)]);
  end for;
  bl,val:=HasSquareIntegralRatValues(Fint);
  if bl then
    return true, [val,1];
  end if;
  for i in [1..#Fint] do
    Fint[i]:=OpX!Reverse(Eltseq(Fint[i]) cat
                   (IsOdd(Degree(Fint[i])) select [0] else []));
  end for;
  bl,val:=HasSquareIntegralRatValues(Fint:Neighbourhood:=0+O(pi^Valuation(Op!p)));
  if bl then
    return true,[1,val];
  else
    return false;
  end if;
end intrinsic;
