\\ ell2coverP(E,P): 2-covers (from ell2cover) NOT already accounted for by the points in P.
\\ Method (Cremona-style filter): for good primes p with E(F_p) cyclic of even order,
\\ the image of C_i(F_p) in E(F_p)/2E(F_p) is a single class c_i(p), and the class of a
\\ rational point Q is its reduction mod p. Solve  Q = sum eps_i*C_i  over F_2.

\\ Class of a point of E(F_p) in E(F_p)/2E(F_p), as a vector of bits, one per even-order
\\ cyclic factor of E(F_p) (G = ellgroup(Ep,,1); works for cyclic and non-cyclic groups).
e2prep(Ep,G) =
{
  my(N=G[1], m=N>>valuation(N,2), gens=vector(#G[2],i,ellmul(Ep,G[3][i],m)),
     ord=vector(#G[2],i,2^valuation(G[2][i],2)), mp=Map(), keep=[], T);
  for(i=1,#G[2], if(ord[i]>1, keep=concat(keep,i)));
  [m,gens,ord,keep];
}

\\ brute force discrete log in the 2-part (orders are small for p <= a few thousand)
e2vec(Ep,pre,Q) =
{
  my(m=pre[1],gens=pre[2],ord=pre[3],keep=pre[4],Q2=ellmul(Ep,Q,m),k=#keep,mp=Map(),R,res,idx);
  if(k==1,
    my(g=gens[keep[1]],T=ellmul(Ep,g,0));
    for(a=0,ord[keep[1]]-1, if(T==Q2, return([a%2])); T=elladd(Ep,T,g));
    error("log failed"));
  \\ k==2
  my(g1=gens[keep[1]],g2=gens[keep[2]],T=[0]);
  for(a=0,ord[keep[1]]-1, mapput(mp,if(#T==1,[0],lift(T)),a); T=elladd(Ep,T,g1));
  T=[0];
  for(b=0,ord[keep[2]]-1,
    R=ellsub(Ep,Q2,T); idx=if(#R==1,[0],lift(R));
    if(mapisdefined(mp,idx), return([mapget(mp,idx)%2,b%2]));
    T=elladd(Ep,T,g2));
  error("log failed");
}

\\ class vector of cover Cv at p; returns [] if no usable point; 0 if inconsistent
covpt(Pm,xv,yv) = [subst(subst(Pm[1],'x,xv),'y,yv), subst(subst(Pm[2],'x,xv),'y,yv)];
covbit(Ep,pre,Cv,p) =
{
  my(R=Cv[1],Pm=Cv[2],b=[],xv,yv,Q,c,r,cnt=0);
  for(xi=0,p-1,
    r=Mod(subst(R,'x,xi),p);
    if(issquare(r),
      yv=sqrt(r); xv=Mod(xi,p);
      Q=iferr(covpt(Pm,xv,yv),err,0);
      if(type(Q)=="t_VEC" && #Q==2 && ellisoncurve(Ep,Q),
        c=e2vec(Ep,pre,Q);
        if(b==[],b=c, if(b!=c,return(0)));
        cnt++; if(cnt>=6,return(b)))));
  b;
}

ell2coverP(E,P,maxp=1000,extra=20) =
{
  my(cov,n,pts=P,T,D,rows=List(),rhs,k,G,Ep,pre,ok,v,Q,den,M,sol,eps,cur,res=List(),bq,nb);
  if(type(E)=="t_VEC" && #E<12, E=ellinit(E));   \\ accept a plain [a1,a2,a3,a4,a6] vector too
  cov=ell2cover(E); n=#cov; T=elltors(E); D=ellglobalred(E)[1]*2;
  for(i=1,#T[3], pts=concat(pts,[T[3][i]]));
  k=#pts; rhs=vector(k,j,List());
  forprime(p=3,maxp,
    if(D%p==0, next);
    Ep=ellinit(E,p); G=ellgroup(Ep,,1); pre=e2prep(Ep,G);
    nb=#pre[4]; if(nb==0, next);
    ok=1; v=vector(n);
    for(i=1,n,
      if(Cv_bad(cov[i],p), ok=0;break);
      v[i]=covbit(Ep,pre,cov[i],p);
      if(type(v[i])!="t_VEC" || #v[i]!=nb, ok=0;break));
    if(!ok, next);
    bq=vector(k);
    for(j=1,k,
      Q=pts[j];
      if(#Q==0, bq[j]=vector(nb); next);      \\ point at infinity
      den=denominator(Q[1])*denominator(Q[2]);
      if(den%p==0, ok=0;break);
      bq[j]=e2vec(Ep,pre,[Mod(Q[1],p),Mod(Q[2],p)]));
    if(!ok, next);
    for(t=1,nb,
      listput(rows,vector(n,i,v[i][t]));
      for(j=1,k,listput(rhs[j],bq[j][t])));
    \\ early stop: full column rank plus 'extra' spare equations (consistency checks)
    if(#rows>=n+extra && matrank(Mod(matrix(#rows,n,r,c,rows[r][c]),2))==n, break));
  if(#rows<n+5, error("too few usable primes (",#rows," equations for ",n," covers); increase maxp"));
  M=Mod(matrix(#rows,n,r,c,rows[r][c]),2);
  if(matrank(M)<n, error("local data cannot separate the covers; increase maxp"));
  eps=vector(k,j,
     sol=matinverseimage(M,Mod(Col(rhs[j]),2));
     if(#sol==0, error("inconsistent: a point is not in the span of the covers"));
     sol);
  cur=matrix(n,0);
  for(j=1,k, cur=matconcat([cur,Mod(eps[j],2)]));
  for(i=1,n,
    my(e=Mod(vector(n,t,t==i)~,2), r0=matrank(cur));
    if(matrank(matconcat([cur,e]))>r0, cur=matconcat([cur,e]); listput(res,cov[i])));
  Vec(res);
}

\\ crude test that p does not divide discriminant data of the quartic / denominators
Cv_bad(Cv,p) =
{
  my(R=Cv[1], d=poldisc(R));
  if(d==0, return(1));
  d%p==0 || denominator(content(R))%p==0 || poldegree(R)<4 || pollead(R)%p==0;
}