#########################################
# biraF for MapleV, releases 4, 5 and 5.1
#########################################

readlib('sturm'):

g_:=a_*X^4+b_*X^3+c_*X^2+d_*X+e_:g_d:=diff(g_,X):

#when p=2, Zp_sol(p) will test for a solution
#X,Y in Z_2 of g_=Y^2 by the algor. in Cremona, Algorithms for
#Elliptic Curves, 2nd ed. sec. 3.6. For p>2, Zp_sol(p) calls
#ZPSOLV which implements Siksek's fast algor. --- see Merriman, Sisek and 
#Smart, Acta Arith. 77('96), p.401.

#gh_:=a_*U^4+b_*U^3*W+c_*U^2*W^2+d_*U*W^3+e_*W^4:
crem_D:=4*I_^3-J_^2:
p_:=3*b_^2-8*a_*c_:s_:=p_^3-48*I_*a_^2*p_-64*J_*a_^3:
DD_:=4*II_^3-JJ_^2:
II_:=12*a_*e_-3*b_*d_+c_^2:
JJ_:=72*a_*c_*e_+9*b_*c_*d_-27*a_*d_^2-27*e_*b_^2-2*c_^3:

Crem:=proc() global  iee;
  iee:=1;crem(args);iee:=0;NULL;
end:

crem:=proc() local dil1, i, l_, r1_0;
global  cq_, crem_dat, den_, dil, dil0, eq_, flag_q, heq_, ncur_, neq_,\
npo2, num_, p, q_, r4, rc, SFE, trou;
  if bouv then aM();RETURN();fi;
  if rc>3 and r4=0 then
    lprint(`Already found the rank r4=0 (using no conjectures)`,qtr);
    RETURN();
  fi;
  ncur_:=ncur;r1_0:=r1;dil:=30;flag_q:=0;
  for i to min(2,nargs) do
    l_:=args[i];
    if type(l_,realcons) then dil:=floor(evalf(l_));
    elif l_=A then flag_q:=1;
    else 
      print(`(Each) args in crem must be a real constant or the letter A.`);
      lprint(`See Menu(crem).`,qtr);
      RETURN();
    fi;
  od;
  if crem_dat=0 then 
    cremo(args);
    if flag_=5 then RETURN();fi;
  else 
    q_:=crem_dat[1];cq_:=crem_dat[2];eq_:=crem_dat[3];neq_:=crem_dat[4];
    dil0:=crem_dat[5];npo2:=crem_dat[6];flag_q:=1;
    if dil<=dil0 then
      fprint(`Quartics already searched for rational points n/d with`);
      eprint(`max(|n|,d) <`=dil0);
      fprint(`To extend search call Crem(x) with x >`=dil0+1);
      RETURN();
    elif iee>0 then
      lprint(`Search resumed for rational points n/d on quartics with`);
      print(dil0<`max(|n|,d) <=`,dil);
    fi;
  fi;
  if r1=0 then heq_:=1;  
  elif nops(RR)=r1 or nops(q_)/2^r1<1/2 then heq_:=1+flag_q;
  else heq_:=0;
  fi;
  dil1:=dil0;
  for den_ to dil while heq_<>1 do
    if den_<=dil0 then
      for num_ from dil0+1 to dil while heq_<>1 do lfs_();od;
    else
      for num_ from 0 to dil while heq_<>1 do lfs_();od;
      if heq_<>1 then dil1:=den_;fi;
    fi;
  od;
  dil0:=dil1;
  if heq_>0 then 
    r4:=r1;rc:=7;i:=(-1)^modp(r4,2);
    if SFE<>0 then
      if SFE<>i then 
        print(QT,` crem`);trou:=[op(trou),cur,r4,SFE];RETURN();
      fi;
    else SFE:=i;
    fi;
    if iee>0 then
      if r1>0 then
        if nops(q_)=0 then 
          lprint(`All quartics were found to be "elliptic", hence`);
        elif nops(q_)/2^r1<1/2 then
          lprint(`Sufficiently many quartics were found to be "elliptic"`);
          lprint(`(the undecided cases are stored in q_), hence`);
        else 
          lprint(`Not all quartics were found to be "elliptic"`);
          lprint(`(the undecided cases are stored in q_), but`);
          lprint(`we have found r1 independent points (stored in RR), hence`);
        fi; 
      fi;
      pr4();
    fi;
    upis();
  else 
    if iee>0 then
      lprint(`Crem (with search limit dil`=dil,`) obtained only the \
upper bound`);
      print(` r1`=r1);
      lprint(`for the rank. The unresolved quartics are stored in q_.`);
      if MM>0 then
        lprint(`The independent point(s) we have so far:`);
        print(` RR`=RR);
      fi;
    fi;
    if r1_0=-1 or r1<r1_0 then upis();fi;
  fi;
  crem_dat:=[q_,cq_,eq_,neq_,dil0,npo2];r1;
end:  

cremo:=proc() local hex,i,l_,ll_, Ng_, p__, r__, r1_, s__;
global cq_, dil0, eq_, flag_, heq_, I_0, J_0, list_, neq_, npo2, p, q_,\
r1, trou, X, x_x;
  if rc>5 and flag_q=0 then
    if r4>0 then lprint(`The rank r4 already determined to be`,r4);fi; 
    if not NNeven then
      lprint(`If you wish to look for the quartics and points on the quartics \
anyway,`);
      lprint(`use Crem(A) or Crem(A,dil).`);i:=`Alternatively y`;
    else i:=Y;
    fi;
    if nops(RR)<r4 and nops(nisog)>0 then
      lprint(i.`ou can use Tri to transfer points from an isog. curve.`);
    fi;
    flag_:=5;RETURN();
  fi;
  do_tor();dil0:=0;
  if NNeven then npo2:=nops(q_1);else npo2:=0;fi;
  if NNeven then
    if flag_q=0 then
      if iee>0 then 
        lprint(`Since E(Q) has a point of order 2, we use RkNC.`);
        lprint(`To overide this, and carry on with crem, use Crem(A) or`);
        lprint(`Crem(A,dil) -- see Menu(Crem).`);
        RkNC(dil);
      else rkNC(dil);
      fi;
      flag_:=5;RETURN();
    else
      fprint(`Overiding the usual choice of 2-descent by 2-isogeny --`);
      fprint(`   will find rep.'s of 2-coverings by quartics.`);
    fi;
  fi;
  I_0:=c4;J_0:=2*c6;flag_:=1:
  if modp(I_0,16)=0 and modp(J_0,64)=0 then I_0:=I_0/16;J_0:=J_0/64;fi;
  crem_(I_0,J_0);Ng_:=1+nops(list_1)+nops(list_2)+nops(list_3);
  if modp(I_0,4)>0 or modp(J_0,8)>0 or modp(2*I_0+J_0,16)>0 then 
    if modp(I_0,4)=0 and modp(J_0,2)=0 then flag_:=3;
    else flag_:=2;
    fi;
    for i to 3 do list_.(i+3):=list_.i;od;
    crem_(16*I_0,64*J_0);
    for i to 3 do
      list_:=list_.i;list_.i:=[];
      for l_ in list_ do
        hex:=0:
        for ll_ in list_.(i+3) do
          p__:=(2/3)*(512*l_[1]*ll_[1]*I_0+4*l_[6]*ll_[6]);
          r__:=8*l_[7]*ll_[7];
          s__:=1024*I_0*(16*l_[1]^2*ll_[6]^2+ll_[1]^2*l_[6]^2+4*l_[1]*ll_[1]*\
l_[6]*ll_[6])-16384*l_[1]*ll_[1]*J_0*(4*l_[1]*ll_[6]+ll_[1]*l_[6])-\
(4*l_[6]*ll_[6])^2:
          if nops(roots(X^4-p__*X^2-8*r__*X+s__/27))>0 then 
            hex:=1:break;
          fi;
        od;
        if hex=0 then list_.i:=[op(list_.i),l_];
        elif iee>0 then
          print(l_[1..5],` discarded --- equiv. to a quartic with the \
smaller I_,J_`);
        fi;
      od;
    od;
    Ng_:=Ng_+nops(list_1)+nops(list_2)+nops(list_3);
  else 
    for i to 3 do list_.(i+3):=[];od;
  fi;
  p:=2;r1_:=vam(Ng_);
  if zz<>1 then 
    print(QT,` crem`);trou:=[op(trou),cur,list_.(1..6)]:RETURN(qR):
  fi;
  if r1_=0 then
    eprint(`The order of the Selmer group S^{(2)} is 1`);
  else
    eprint(`The order of the Selmer group S^{(2)} is 2`^r1_);
  fi;
  r1_:=r1_-min(2,npo2);
  fprint(`This gives the (conjecture-free) upper bound`,r1_,`for the rank.`);
  if r1<>-1 and r1<r1_ then
    if iee>0 then
lprint(`But the upper bound`,r1,`was previously determined by another method.`);
      if rc<4 then 
        lprint(`(Since rc`=rc,`this depends on standard conjecture(s).)`);
      fi;
lprint(`Thus we expect the 2-dimension of the Shafarevich-Tate group`);
lprint(`to be at least`,r1_-r1,`and some of the loc. sol. quartics`);    
lprint(`will not be "elliptic" --- will not have a rational point.`);
    fi;
  else r1:=r1_;
  fi;
  q_:=[];cq_:=[];eq_:=[];neq_:=[];
  if r1=0 then RETURN();
  elif r1=MM then 
    if flag_q=1 then
fprint(`We now know (using no conj.'s) that rank`=r1,`, and RR contains`);
fprint(`that many indep. points. But we look for points on the quartics`);
      fprint(`anyway, in an attempt to improve RR.`);
    else RETURN();      
    fi;
  elif iee>0 then
    lprint(`We now attempt to determine which quartics are "elliptic", \
i.e., have`);
    lprint(`a rational point and so are elliptic curves iso. with E.`);
  fi;
  for i to 6 do 
    for x_x in list_.i do
      set_(op(x_x[1..5]));
      if essqr(e_) then quar1(0,esqrt);
      elif issqr(a_) then quar1(0);
      else
        cq_:=[op(cq_),x_x[1..5]];q_:=[op(q_),g_];
      fi;  
    od;
  od;
end:

crem_:=proc(a,b) local i;global  I_, J_, phi_;
  I_:=a:J_:=b;
  for i to 3 do list_.i:=[];od;
  phi_:=sort([fsolve(X^3-3*I_*X+J_)]);
  if crem_D<0 then crem_3();
  else crem_1();crem_2();
  fi;
end:

crem_1:=proc() local k_,rk_;global  type_, list_1;
  type_:=1;k_:=(4*I_-phi_[3]^2)/3:rk_:=sqrt(k_);  
  list_1:=step_2(1,trunc((k_+rk_*phi_[3])/(3*rk_+phi_[3]+2*phi_[2])+.0001));
end:

crem_2:=proc() local k_;global  type_, list_2;
  type_:=2;k_:=(I_-phi_[2]^2)/3;
  list_2:=step_2(ceil(k_/(phi_[2]-phi_[3])-.0001),floor(k_/(phi_[2]-phi_[1])+.0001));
end:

crem_3:=proc() local k_;global  type_, list_3;
  type_:=3;k_:=sqrt((4/27)*(phi_[1]^2-I_));
  list_3:=step_2(ceil(phi_[1]/3-k_-.0001),floor(phi_[1]/3+k_+.0001));
end:

Cubic_fit:=proc()
global  iee;
;
  iee:=1;cubic_fit(args);iee:=0:NULL;
end:

cubic_fit:=proc() 
global  gg_, s_10, x_, y_, f_;
local e,i,N,p,sols:
  if nargs=0 then print(`Cubic_fit requires args.`):RETURN():fi:
  N:=nargs:gg_:=args[nargs]:prep_():s_10:=0:
  if type(gg_,{equation,list(equation),set(equation)}) then
     N:=N-1:e:=traperror(assign(gg_)):
     if e=lasterror then 
       prep_();ERROR(`PLEASE REPEAT THE SAME cubic_fit COMMAND.`);
     fi;
  fi:
  if N>0 then
    gg_:={}:
    for i to N do
      p:=args[i]:
      if not type(p,list) or nops(p)<>2 then print(qts):RETURN():fi:
      x_:=p[1]:y_:=p[2]:
      gg_:=gg_ union {f_}:
    od:
    sols:={solve(gg_)}:
    if sols={} then print(`No solution.`):f_:=NULL;RETURN():fi:
    assign(sols):
  fi;
  x_:='x_':y_:='y_':
  if iee=1 then
    lprint(`The general cubic satisfying the requirements is`);
    print(f_);
lprint(`One can now do arithmetic on the cubic using Third_pt, etc., or,`);
lprint(`after specializing the remaining s_.i (if any) to rational values,`);
lprint(`one can call Gcub(s_.(1..10)), to change to the min. Weier. equation.\
`);
  fi;
  f_;
end:

eub:=proc(fub) local dub, i, j;global  Uo, Vo;
#(internal) used by gcub 3X
  dub:=collect(sex_(fub),XR):
  for i from 0 to 3 do
    for j from 0 to i do
      s.(10-j-i*(i+1)/2):=coeff(collect(coeff(dub,XR,j),YR),YR,i-j):
    od;
  od;
  Uo:='Uo':Vo:='Vo':
end:

Gcub:=proc() global  iee;
  iee:=1:gcub(args);iee:=0;NULL;
end:

gcub:=proc() local ape, gc_2, i, j, nergs, s;
global  a1, a2, a3, a4, a5, a6, bec, bic, boc, ca, czb, d, e1, e2, e3, e4, e5, 
GC_2, GC_p, iee, K_, new_bec, new_bic, p_at_inf_c, p_at_inf_w, r, s11, s12, 
sBb, stoRE, storE, t, T0, ta, tim, tom, U_, Uo, uV, V_, Vo, Weq, Wr, Ws, Wt, 
Wu, X, X_, xww, Y, Y_, yww, z, zQ;
;
#Input (general) cubic and convert to weierstrass form
  if type(GC_2,boolean) and GC_2 then gc_2:=false;GC_2:=false;
  else gc_2:=true;
  fi;
  if nargs<>10 and nargs<>12 and nargs<>13 then
    print(`gcub requires 10, 12 or 13 args --- see Menu(gcub)`);RETURN();
  fi;
  s:=[];nergs:=nargs;
  for i to nergs do s.i:=sex_(args[i]):s:=[op(s),s.i];od;
  if type(s,list(rational)) then K_:=1;
  elif nergs=10 then 
print(`With irrational args, gcub requires that a point be given; see Menu(\
gcub)`);
    RETURN();
  elif type(s,list(ratpoly(rational,T))) then K_:=T;
  elif type(m_,integer) and type(s,list(ratpoly(rational,w_))) then K_:=m_;
  else K_:=0;
  fi;
  if U_=X and V_=Y then
    Uo:=U:Vo:=V:czb:=cub:
    fprint(`The equation of the curve is`);eprint(sort(cub,[U,V])=0);
  fi;
  Uo:='Uo':Vo:='Vo':
  if cu3=0 then print(`curve is of genus 0`,qtr);RETURN();fi;
  if nergs=13 then
    Uo:=s11:Vo:=s12:
    if s13<>0 or (Uo=0 and Vo=0) then 
      print(`point at infinity must be specified by homogeneous coordinates`);
      print(`a,b,0 as args[11..13]`);RETURN();
    fi;
    if cu3<>0 then print(`[Uo,Vo,0] not on curve `,qtr);RETURN();fi;
  fi;
  if (nergs=10 and s1=0) or (nergs=13 and Vo=0) then
    unk(1/XR,YR/XR);xnk(1/UR,VR/UR);
    fprint(`The point at infinity with homogeneous coordinates [1,0,0]`);
    fprint(`will be taken as the group O.`);GC_2:=true;GC_p:=[1,0,0];
    gcub(s10,s9,s7,s4,s8,s6,s3,s5,s2,s1,0,0);RETURN();
  elif nergs<>12 then
    if nergs=10 then Uo:='Uo':Vo:=1:e1:=[solve(cu3,Uo)]:
    else e1:=[Uo/Vo]:
    fi;
    for i to nops(e1) do
      if type(e1[i],rational) then 
        Uo:=e1[i]:unk(XR/YR,1/YR);xnk(UR/VR,1/VR);GC_p:=[Uo,1,0];
        fprint(`The point at infinity with homogeneous coordinates`,GC_p);
        fprint(`will be taken as the group O.`);GC_2:=true;
        gcub(s1,s5,s8,s10,s2,s6,s9,s3,s7,s4,Uo,0);RETURN();
      fi;
    od;
    fprint(`Beginning search for rational point.`);
    if s10=0 then z:=[0,0]:else z:=NULL:fi;
    if tim>0 then eprint(`Clock started`);tom:=time():fi;
    for t from 0 to lgcub while z=NULL do
      for i from -t to t while z=NULL do
        Uo:=i:Vo:=t-abs(Uo):e1:=[solve(cuh,Z)]:
        for j to nops(e1) do
          d:=e1[j]:
          if type(d,rational) and d<>0 then z:=[Uo/d,Vo/d]:break;fi;
        od;
        if z=NULL and Vo<>0 then
          Vo:=-Vo:e1:=[solve(cuh,Z)]:
          for j to nops(e1) do
            d:=e1[j]:
            if type(d,rational) and d<>0 then z:=[Uo/d,Vo/d]:break;fi;
          od;
        fi;
      od;
      if tim>0 then clk();
        if tim<0 then RETURN();fi;
      fi;
    od;
    tim:=0:
    if z=NULL then
      print(`Found no rational point up to |num(U)|+|num(V)|`=lgcub,qtr);
      RETURN();
    else
      fprint(`Found point`,z);s11:=z[1]:s12:=z[2]:
    fi;
  fi;
  Uo:=s11:Vo:=s12:z:=[Uo,Vo]:
  if gc_2 then GC_p:=z;fi;
  if U_=X and V_=Y then
    if cub<>0 then
      lprint(`[U,V] `=z,noc,qtr);RETURN();
    fi;
    fprint(`The point P0 `=z,`will become the point at infinity on the`);
    fprint(`birationally equivalent Weierstrass form.`);
    fprint(`Let the tangent at P0 meet the cubic in the third point Q.`);
  fi;
  if cbv=0 then
    if cbu=0 then
      lprint(`P0 `=z,`is a singular point`,qtr);RETURN();
    else 
      unk(YR,XR);xnk(VR,UR);gcub(s4,s3,s2,s1,s7,s6,s5,s9,s8,s10,s12,s11);
      RETURN();
    fi;
  fi;
  Uo:=sex_(s11+XR):Vo:=sex_(s12+YR):unk(Uo,Vo);xnk(UR-s11,VR-s12);eub(cub);
  Uo:=s9:Vo:=-s8:e2:=cu2:e3:=cu3:T0:=sex_(-s8/s9):
  if e3=0 then
    if e2=0 then 
      lprint(`Curve is reducible - tangent at P0 is a component`,qtr);    
      U_:=X:V_:=Y:X_:=U:Y_:=V:RETURN();
    fi;
    fprint(`Q is at infinity.`);Uo:=sex_(1-s9/XR):Vo:=sex_((YR+s8)/XR):
    unk(Uo,Vo);xnk(s9/(1-UR),VR*s9/(1-UR)-s8);eub(XR^3*cub);
  else
    e4:=sex_(-s9*e2/e3):e5:=sex_(s8*e2/e3):
    if e2=0 then fprint(`P0=Q is a flex.`);
    else
      Uo:=denom(U_):Vo:=denom(V_):X:=e4:Y:=e5:
      if Uo=0 or Vo=0 then zQ:=`point at infinity`:
      else zQ:=[U_,V_]:
      fi;
      X:='X':Y:='Y':fprint(`Q = `,zQ);
      Uo:=XR+e4:Vo:=YR+e5:unk(Uo,Vo);xnk(UR-e4,VR-e5);eub(cub);
    fi;
  fi;
  Uo:=1:ta:='ta':Vo:=T0+1/ta:r:=collect(sex_(ta^4*(cu2^2-4*cu1*cu3)),ta):
  ca:=coeff(r,ta,3):
  if ca=0 then print(`curve is of genus 0`,qtr);RETURN();fi;
  a1:=0:a2:=coeff(r,ta,2):a3:=0:
  a4:=sex_(ca*coeff(r,ta,1)):a6:=sex_(ca^2*coeff(r,ta,0)):a5:=a6:
  if K_=1 then con9();
  else Wu:=1;Wr:=0;Ws:=0;Wt:=0;uV:=1;
  fi;
  xww:=(Wu^2*XR+Wr)/uV^2:ta:=xww/ca:
  yww:=(Wu^3*YR+Wu^2*Ws*XR+Wt)/uV^3;
  Uo:=((ca*yww/xww^2)-cu2)/(2*cu3):unk(Uo,Vo*Uo);
  ta:=VR/UR:i:=1/(ta-T0):Uo:=1:Vo:=VR/UR:
  xnk((ca*i*uV^2-Wr)/Wu^2,(ca*i^2*(2*cu3*UR+cu2)*uV^3-Ws*ca*i*uV^2+Ws*Wr-Wt\
)/Wu^3):
  if K_=1 then con0();
  else 
    ape:=iee:iee:=0:stoRE:=storE:storE:=true:uV:=1:
    ell(a1,a2,a3,a4,a6,K_);store(stoRE):iee:=ape:
    if DD=0 then print(`curve has genus 0`,qtr);
    else
      X_:=simplify(X_,{czb},{U,V});Y_:=simplify(Y_,{czb},{U,V});
      U_:=normal(U_):V_:=normal(V_):X_:=normal(X_):Y_:=normal(Y_):
      Weq:={Y^2=-a1*X*Y-a3*Y+X^3+a2*X^2+a4*X+a6};
      boc:=[simplify(numer(U_),Weq,[Y,X]),denom(U_),simplify(numer(V_),Weq\
,[Y,X]),denom(V_),numer(X_),denom(X_),numer(Y_),denom(Y_),GC_p]:
      boc:=sex_(boc);
      if uV<>1 then 
        boc:=subs(X=X/uV^2,Y=Y/uV^3,boc);
        boc:=[op(boc[1..4]),uV^2*boc[5],boc[6],uV^3*boc[7],boc[8],boc[9]];
      fi;
      bic:=[op(bic),boc]:new_bic:=bic;
      eprint(cus,cur=yo);
      buc(op(boc));bec:=[op(bec),czb]:new_bec:=bec;sBb:=0;
      p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]];
    fi;
    U_:=X:V_:=Y:X_:=U:Y_:=V:   
  fi;
  if iee=0 then print(cus,cur=yo);fi;
  NULL;
end:

hens:=proc(La,Lb,Ah,Bh) global bh,bbh,hex,Lc,Ld,vh,Lu,r4,R4;
  Lc:=Lb:Lu:=La union Lb;
  while (r3<r4 or fleg=1) and nops(Lc)>0 do
    bh:=Lc[1];bbh:=Bh/bh;
    if bh<0 and Ah<=0 and Bh>0 then hex:=0;
    else hex:=1;set_(bh,0,Ah,0,bbh);
    fi;
    for vh in SS while hex=1 do
      if (vh>2 and modp(bh,vh)>0 and L(bh,vh)=1) or \
         (vh=2 and (modp(bh,8)=1 or modp(bh+4*Ah,8)=1)) then
        next;
      fi;
      hex:=Zp_sol(vh);
    od;
    Ld:=map(proc(x) x*bh/igcd(x,bh)^2 end,La):Lc:=Lc minus Ld:
    if hex=0 then
      if fleg=0 then R4:=R4/nops(Lu):fi;
      Lu:=Lu minus Ld:
      if fleg=0 then R4:=R4*nops(Lu):r4:=min(r4,sam(R4)):fi;
    fi;
  od;
end:

lemma7:=proc(x,n) global  X, a, b, l, c, v;
  X:=x;a:=g_:b:=g_d;X:='X';
  if a=0 then RETURN(true);fi;
  l:=vam(a);
  if modp(l,2)=0 and modp(zz,8)=1 then RETURN(true);fi;
  if n=1 then c:=1;else c:=modp(zz,4);fi;
  if b=0 then v:=n:
  else
    vam(b);
    if l>2*v or (n>v and l=v+n-1 and modp(l,2)=0) then RETURN(true);
    elif n>v and l=v+n-2 and c=1 and modp(l,2)=0 then RETURN(true);
    fi;
  fi; 
  if v>=n and (l>=2*n or (l=2*n-2 and c=1)) then RETURN(FAIL);fi;
  false;
end:

lfs_:=proc() global  heq_;
  if igcd(num_,den_)=1 then
    lft_(num_/den_);
    if num_<>0 and heq_<>1 then lft_(-num_/den_);fi;
    if den_=dil then heq_:=min(heq_,1);fi;
  fi;
end:

lft_:=proc() local ergs,i,j,k;
global  X, ss_, iee_, iee, trou, neq_, eq_, x, y, q_, cq_, heq_;
;
  ergs:=args;X:=args;ss_:=[];
  for j to nops(eq_) do
    if essqr(eq_[j]) then
      k:=iee:iee:=0;trcw([X,esqrt],neq_[j]);iee:=k;
    fi;
  od;
  for k to nops(q_) do
    if essqr(q_[k]) then 
      iee_:=iee;iee:=-1;quar(op(cq_[k]),X,esqrt);iee:=iee_;
      if ncur<>ncur_ then
        X:='X';print(QT,` lfs_`);trou:=[op(trou),ncur_,q_[k]]:
        RETURN(qR):
      else 
        neq_:=[op(neq_),nops(bic)];X:='X';eq_:=[op(eq_),q_[k]];X:=ergs;
      fi;
    else ss_:=[op(ss_),k];
    fi;
  od;
  X:='X';
  if nops(ss_)<nops(q_) then
    x:=NULL;y:=NULL;
    for i in ss_ do x:=x,q_[i];y:=y,cq_[i];od;
    q_:=[x];cq_:=[y];
  fi;
  if nops(RR)=r1 or nops(q_)/2^r1<1/2 then heq_:=1+flag_q;fi;
end:

loc_sol:=proc() local x;global  P__, P_;
  if issqr(a_) or issqr(e_) then RETURN(true);
  elif a_<0 and sturm(g_,X,-infinity,infinity)=0 then 
    P__:=infinity;RETURN(false);
  fi;
  for P_ in {op(badp_)} union {2,3} do
    x:=Qp_sol(P_);
    if x=false then P__:=P_;RETURN(false);fi;
  od;
  true;
end:

prep_:=proc()
global  s_1, s_2, s_3, s_4, s_5, s_6, s_7, s_8, s_9, s_10, x_, y_, f_;
;
s_1:='s_1':
s_2:='s_2':
s_3:='s_3':
s_4:='s_4':
s_5:='s_5':
s_6:='s_6':
s_7:='s_7':
s_8:='s_8':
s_9:='s_9':
s_10:='s_10':
x_:='x_':y_:='y_':
f_:=s_1*x_^3+s_2*x_^2*y_+s_3*x_*y_^2+s_4*y_^3+s_5*x_^2+s_6*x_*y_+\
   s_7*y_^2+s_8*x_+s_9*y_+s_10:
end:

qoo:=proc() local ie;
global  cuc, iee, new_bec, new_bic, p_at_inf_c, p_at_inf_w, sBb, U_, V_, 
Weq, X_, Y_;
  U_:=normal(sex_(U_)):
  V_:=normal(sex_(V_)):
  X_:=normal(sex_(X_)):
  Y_:=normal(sex_(Y_)):
  Weq:={Y^2=-a1*X*Y-a3*Y+X^3+a2*X^2+a4*X+a6};
  cuc:=[simplify(numer(U_),Weq,[Y,X]),denom(U_),simplify(numer(V_),Weq\
,[Y,X]),denom(V_),numer(X_),denom(X_),numer(Y_),denom(Y_)];
  cuc:=sex_(cuc);#Yes, another simplify is sometimes needed when K_=m_
  U_:=X:V_:=Y:X_:=U:Y_:=V:
  new_bec:=[sex_(czb)];new_bic:=[cuc];
  ie:=iee:iee:=0:ell(args);iee:=ie:
  sBb:=0;p_at_inf_c:=[op(p_at_inf_c),[]];p_at_inf_w:=[op(p_at_inf_w),[]];
  if K_=T then iee:=0;tesT();iee:=ie;fi;
  if iee=1 then 
    lprint(`Weierstrass equation:`);
    print(sort(Y^2+a1*X*Y+a3*Y,[Y,X])=sort(X^3+a2*X^2+a4*X+a6,X));
    print(`Field label K_`=K_);
  fi;
  iee:=ie;buc(op(cuc));
end:

Qp_sol:=proc(P_) global qP_;
  if Zp_sol(P_)=1 then RETURN(true);fi;
  set_(e_,d_,c_,b_,a_);qP_:=Zp_sol(P_);set_(e_,d_,c_,b_,a_);
  if qP_=1 then qP_:=true;
  else qP_:=false;
  fi;
end:

Quar:=proc() global  iee;
  iee:=1:quar(args);iee:=0:NULL;
end:

quar:=proc() local ergs,hex,j,k,k2,qboc,qum,quT,quV,qu_2,sI,x,y;
global a1,a2,a3,a4,a5,a6,a_3,a_4,a_5,a_6,a_7,badp_,bic,boc,czb,d,h,i,iee,
K_,quxr,quyr,quzr,QU_2,QU_p,RR,s,S6,S7,S8,t,T2,T3,T4,T5,tim,tom,
Uo,Us,Ut,uV,Vt,Weq,Wr,Ws,Wt,Wu,Xt,Yt,z_q;
  if nargs<5 or nargs>7 then
    print(`quar and quar0 require 5, 6 or 7 args`,qtr);RETURN(-2);
  fi;
  if type(QU_2,boolean) and QU_2 then qu_2:=false;QU_2:=false;
  else qu_2:=true;
  fi;
  quxr:=true:ergs:=args;S6:='S6';S7:='S7';uV:=1;
  if nargs=6 then ergs:=ergs[1..5],0;fi;
  for i to nargs do
    S.i:=sex_(ergs[i]):quxr:=quxr and type(S.i,rational):
  od;
  sI:=[S.(1..nargs)];
  quyr:=quzr and quxr:qum:=false;
  if not quyr then K_:=0;fi;
  if type(m_,integer) and type(sI,list(polynom(rational,w_))) then
    if quarm_flag or has(sI,w_) then 
      qum:=true;K_:=m_;
    else 
      if K_=m_ and quxr then
        eprint(`Changing field label to K_=1; to keep label`,m_);
        eprint(`for quartic with rational coef.'s, use commands Qfin, Quarm.`);
        K_:=1;
      fi;
    fi;
  fi;
  if type(sI,list(ratpoly(rational,T))) and has(sI,T) then quT:=true;K_:=T;
  else quT:=false;
  fi;  
  j:=`The quartic has discriminant 0, hence the curve has genus 0`,qtr;
  if S1=0 then
    if S2=0 then print(j);
    elif quyr then 
      fprint(`Since coeff. of U^4 is 0, calling Ein0:`);ein0(S2,0,S3,0,S4,S5);
    else 
      fprint(`Since coeff. of U^4 is 0, calling Ell with field label K_`=K_);
      ell(0,S3/S2,0,S4/S2,S5/S2,K_);
    fi;
    RETURN([0]);
  fi;
  set_(S.(1..5));
  if sex_(DD_)=0 then
    print(j);RETURN(0);
  fi;
  if U_=X and V_=Y then 
    Uo:=U:eprint(`curve is   `,V^2=sort(qu,U));czb:=qu-V^2:
  fi;
  Uo:='Uo':
  if nargs=5 then
    hex:=true;
    if essqr(S5) then hex:=false:S6:=0:S7:=esqrt:fi;
    if hex then
      if qum then k:=w_;else k:=NULL;fi;
      for j in factors(numer(normal(qu)),k)[2] while hex do
        k:=j[1];
        if degree(k,Uo)=1 then
          hex:=false;S6:=sex_(-coeff(k,Uo,0)/coeff(k,Uo,1));S7:=0;
        fi;
      od;
    fi;
    if hex then
      if essqr(S1) then
        eprint(`The places at infinity are rational - we take one as O.`);
        QU_2:=true;QU_p:=[0,1,0];unk(1/XR,YR/XR^2);xnk(1/UR,VR/UR^2);
        if not quzr then quar0(S5,S4,S3,S2,S1,0,esqrt);
        elif quarm_flag then quarm(S5,S4,S3,S2,S1,0,esqrt);
        else quar(S5,S4,S3,S2,S1,0,esqrt);
        fi;
        RETURN([]);
      elif quxr then
        k:=ilcm(seq(denom(S.j),j=1..5));
        if k>1 then
          k:=k/ofad(k,2):k2:=k^2;
          for j to 5 do S.j:=k2*S.j;od;
          if iee>0 then quar2(U,V/k);fi;
          if quarm_flag then RETURN(quarm(S.(1..5)));
          else RETURN(quar(S.(1..5)));
          fi;
        fi;     
        fprint(`Beginning search for rational point.`);
        badp_:=factorset(abs(DD_));
        if tim>0 then eprint(`Clock started`);tom:=time():fi;
        if not loc_sol() then 
          eprint(`The quartic has no point defined over Q since`);
          eprint(`it is not locally solvable at`,P__,qtr);
          RETURN(P__);
        fi;
        for i from 3 to 7 do a_.i:=false:od;
        if modp(a_,4)>1 then a_4:=true;fi;
        for i in {3,5,7} do 
          if L(a_,i)=-1 then a_.i:=true:fi;
        od;
        for h to lquar while hex do
          for d to h while hex do
            i:=h+1-d;
            if igcd(i,d)>1 then next;fi;
            for j from 3 to 7 do
              if modp(d,j)=0 and a_.j then next;fi;
            od;          
            Uo:=i/d:s:=qu:Uo:=-Uo:t:=qu:Uo:='Uo':
            if essqr(s) then S6:=i/d:S7:=esqrt:hex:=false:
            elif essqr(t) then S6:=-i/d:S7:=esqrt:hex:false:
            fi;
          od;
          if tim>0 then clk();
            if tim<0 then RETURN(FAIL);fi;
          fi;
        od;
        tim:=0:
        if hex then
          if not quzr then
            hex:=false;S6:=0:S7:=simplify(sqrt(S5),sqrt):
          else 
            fprint(`No point with U=a/b found in the range |a|+b<=lquar`);
            fprint(`where the present value of lquar (which you can reset`);
            fprint(`by a command lquar:=...) is`,lquar,qtr);
            RETURN(FAIL);
          fi;
        fi;
      fi;
      if hex then
        if qum or quT then
          fprint(`Didn't notice a rational point on quartic`);
          if quarm_flag then fprint(qtr);RETURN(FAIL);
          else 
fprint(`so in the transition to Weier. form, the present field label K_`=K_);
fprint(`stands a good chance of being reset to 0.`);
          fi;
        fi;
        S6:=0:S7:=simplify(sqrt(S5),sqrt,symbolic):
      fi;
    fi;
    eprint(`Taking point on quartic to be `,[S6,S7]);QU_p:=[S6,S7];
  elif nargs=6 then
    if essqr(S1) then
      unk(1/XR,YR/XR^2);xnk(1/UR,VR/UR^2);QU_2:=true;QU_p:=[0,1,0];
      quar(S5,S4,S3,S2,S1,0,esqrt);
      RETURN([]);
    else
      fprint(`Places at infinity not rational`,qts);RETURN(FAIL);
    fi;
  fi;
  T5:=S7^2:Uo:=S6:s:=qu:Uo:='Uo':
  if sex_(s-T5)<>0 then 
    print([S6,S7],noc,qtr);RETURN(-1);
  fi;
  T2:=S2+4*S1*S6:T3:=S3+3*S2*S6+6*S1*S6^2:
  T4:=S4+2*S3*S6+3*S2*S6^2+4*S1*S6^3:
  qboc:=boc;boc:=[S6,S7];quV:=uV;
  if qu_2 then QU_p:=boc;fi;
  if S7<>0 then
    a1:=T4/S7:a2:=T3-T4^2/(4*T5):a3:=2*S7*T2:a4:=-4*T5*S1:a6:=a2*a4:S8:=2*S7:
    a5:=a6:Ut:=UR-S6:Vt:=VR+S7:
    Xt:=(S8*Vt+T4*Ut)/Ut^2:Yt:=(4*T5*Vt+S8*(T4+T3*Ut)*Ut-(T4*Ut)^2/S8)/Ut^3:
    x:=-a2;y:=a1*a2-a3;
    if quyr then 
      con9();x:=x*uV^2-Wr;z_q:=[x/Wu^2,(y*uV^3-Ws*x-Wt)/Wu^3];
    else 
      Wu:=1:Wr:=0:Ws:=0:Wt:=0:uV:=1:
      for i to 6 do a.i:=sex_(a.i);od;    
      z_q:=sex_([x,y]);
    fi;
    xnk((Xt*uV^2-Wr)/Wu^2,(Yt*uV^3-Ws*Xt*uV^2+Ws*Wr-Wt)/Wu^3);
    i:=(Wu^2*XR+Wr)/uV^2:Us:=(S8*(i+T3)-T4^2/S8)*uV^3/(Wu^3*YR+Wu^2*Ws*XR+Wt):
    unk(Us+S6,-S7+Us*(Us*i-T4)/S8);boc:=qboc;uV:=quV;
    if quyr then con0();else qoo(a.(1..5),K_):fi;
  else
    unk(S6+1/XR,YR/XR^2);xnk(1/(UR-S6),VR/(UR-S6)^2);
    boc:=qboc;uV:=quV;z_q:=[];
    if quyr then ein0(T4,0,T3,0,T2,S1);
    else unk(XR/T4,YR/T4):xnk(T4*UR,T4*VR):qoo(0,T3,0,T2*T4,S1*T4^2,K_):
    fi;      
  fi;
  z_q:=sex_(z_q);
  if quyr then ap0(1);
  elif z_q<>[] and (K_=T or K_=m_) then z_q:=[z_q[1]*uV^2,z_q[2]*uV^3];
  fi;
  if z_q<>[] and not member(z_q,uR) then
    eprint(`We consider the point `,z_q);ap1_(z_q);
  fi;
#  fprint(`Present field label K_`=K_);
  if NN=-1 and nops(RR)>0 then 
eprint(`(but this list may contain torsion point(s) since Tor() hasn't bee\
n called).`);
  fi; 
#  If K_=m_ and uV<>1 then the adjustments to RR and bic are now 
#    done in oin.
  j:=nops(bic);bic:=[op(bic[1..(j-1)]),[op(bic[j]),QU_p]];
  if iee=0 then print(cus,cur=yo);fi;
  [S6,S7];
end:

Quarm:=proc() global iee;
  iee:=1;quarm(args);iee:=0;NULL;
end:

quarm:=proc() local q;global quarm_flag, quzr;
  if not type(m_,integer) then 
    lprint(`Use Qfin to initialize apecs`,qts);RETURN();
  elif not type(sex_([args]),list(polynom(rational,w_))) then
    lprint(`args not defined over Q(sqrt(m_), where m_`=m_,qts);
    RETURN();
  fi;
  quarm_flag:=true;quzr:=false;q:=quar(args);quarm_flag:=false;quzr:=true;q;
end:

Quar0:=proc() global  iee;
  iee:=1:quar0(args):iee:=0:NULL;
end:

quar0:=proc() local qu;global  quzr;
  quzr:=false:qu:=quar(args):quzr:=true:qu;
end:

quar1:=proc() local nergs,iee_;global iee,eq_,neq_;
  nergs:=nargs;iee_:=iee:iee:=-1;quar(op(x_x[1..5]),args);iee:=iee_;
  if nergs=1 then set_(op(x_x[1..5]));fi;
  eq_:=[op(eq_),g_];neq_:=[op(neq_),nops(bic)];
end:

quar2:=proc() local x;
  x:=args;
  lprint(`Replacing U,V with `,x,` -- see new input quartic below.`);
  lprint(`N.B. the commands trwc and trcw will refer to points [U,V] o\
n the new`);
  lprint(`quartic. The coord.'s on the old quartic above -- call them U1,V1 -\
--`);
  lprint(`can be calc. from `,[`U1`,`V1`] = [x]);
  lprint(`New quartic (now referred to as "the original equation"):`);
end:

set_:=proc(a,b,c,d,e)
global  a_, b_, c_, d_, e_;
  a_:=a;b_:=b;c_:=c;d_:=d;e_:=e;
end:

step_2:=proc(a_1,a_2) local  a__, b_1, b_3, b__, c_1, c_2, c__, c_3, p__, l_;
global list_, a_, b_, c_;
  list_:=[];c_3:=3;
  if flag_=1 then b_3:=1;else b_3:=4;fi;
  if iee>0 then
    lprint(`Searching for non-triv., inequiv., loc. sol. quartics \
a_X^4+b_X^3+...+e_`);
    lprint(`of type `.type_,` with a_1<=a_<=a_2 where I_, J_, a_1, a_2 =`);
    print(I_, J_, a_1, a_2);
  fi;
  for a__ from a_1 to a_2 do
    if a__=0 then next;fi;
    a_:=a__:
    if flag_>1 then
      if modp(a_,4)=0 or (flag_=3 and modp(a_,2)=0) then next;
      elif modp(a_,4)=2 then c_3:=12;
      else c_3:=6;
      fi;
    fi;
    b_1:=-2*abs(a_);
    if flag_>1 and modp(b_1,4)=2 then b_1:=b_1+2;fi;
    for b__ from b_1 to -b_1 by b_3 do
      b_:=b__;
      if type_=1 then 
        c_1:=ceil((4*a_*phi_[2]+3*b_^2)/(8*a_)-.0001);
        c_2:=floor((4*a_*phi_[3]+3*b_^2)/(8*a_)+.0001);
      elif type_=2 then
        if a_>0 then
          c_1:=ceil((4*a_*phi_[2]-(4/3)*(I_-phi_[2]^2)+3*b_^2)/(8*a_)-.0001);
          c_2:=floor((4*a_*phi_[1]+3*b_^2)/(8*a_)+.0001);
        else
          c_1:=ceil((4*a_*phi_[3]+3*b_^2)/(8*a_)-.0001);     
          c_2:=floor((4*a_*phi_[2]-(4/3)*(I_-phi_[2]^2)+3*b_^2)/(8*a_)+.0001);
        fi;
      else
        if a_>0 then 
          c_1:=ceil((9*a_^2-2*a_*phi_[1]+(1/3)*(4*I_-phi_[1]^2)+3*b_^2)/\
(8*a_)-.0001);
          c_2:=floor((4*a_*phi_[1]+3*b_^2)/(8*a_)+.0001);
        else          
          c_1:=ceil((4*a_*phi_[1]+3*b_^2)/(8*a_)-.0001);
          c_2:=floor((9*a_^2-2*a_*phi_[1]+(1/3)*(4*I_-phi_[1]^2)+3*b_^2)/\
(8*a_)+.0001);
        fi;
      fi;
      c_1:=c_1-1+c_3-modp(c_1-4*b2-1,c_3);
      if c_1>c_2 then next;fi;
      for c_ from c_1 to c_2 by c_3 do
        if issqr(s_/27) then
          step_2A(isqrt(s_/27));#step_2A(-isqrt(s_/27));
        fi;
      od;
    od;
  od;
  if iee>0 then
    if nops(list_)=0 then lprint(`None found.`);
    else 
      lprint(`Found`);p__:=NULL;
      for l_ in list_ do p__:=p__,l_[1..5]:od;
      print(p__);
    fi;
  fi;
  list_;
end: 

step_2A:=proc() local hex, l_, p__, r__, s__;
global  badp_, d_, e_, list_, r_;
  r_:=args;d_:=(r_-b_^3+4*a_*b_*c_)/(8*a_^2);
  if not type(d_,integer) then RETURN();fi;
  e_:=(I_+3*b_*d_-c_^2)/(12*a_);
  if (not type(e_,integer)) or nops(roots(g_))>0 then RETURN();fi;
  hex:=0;
  for l_ in list_ do
    p__:=(2/3)*(32*l_[1]*a_*I_+l_[6]*p_);
    r__:=l_[7]*r_;
    s__:=64*I_*(l_[1]^2*p_^2+a_^2*l_[6]^2+l_[1]*a_*l_[6]*p_)\
-256*l_[1]*a_*J_*(l_[1]*p_+a_*l_[6])-(l_[6]*p_)^2:
    if nops(roots(X^4-p__*X^2-8*r__*X+s__/27))>0 then 
      hex:=1:break;
    fi;
  od;
  if hex=0 then
    badp_:=factorset(abs(DD_));
    if loc_sol() then list_:=[op(list_),[a_,b_,c_,d_,e_,p_,r_]];fi;  
  fi;
end:

Third_pt:=proc()
global  iee;
;
  iee:=1;third_pt(args);iee:=0;NULL;
end:

third_pt:=proc() 
global  x_1, y_1, x_2, y_2, i, x_, y_, gg_, dx_, dy_, x_3, y_3;
local m;
  if not type(f_,polynom) or not has(f_,x_) then
    print(`First use Cubic_fit(..)`):RETURN():
  fi:
  if nargs=4 then
    x_1:=args[1]:y_1:=args[2]:x_2:=args[3]:y_2:=args[4]:
  elif nargs=2 then
    x_1:=args[1]:
    if type(x_1,list) then
      x_1:=args[1][1]:y_1:=args[1][2]:x_2:=args[2][1]:y_2:=args[2][2]:
    else
      x_1:=args[1]:x_2:=x_1:y_1:=args[2]:y_2:=y_1:
    fi:
  else print(qts):RETURN():
  fi:
  for i to 2 do
    x_:=x_.i:y_:=y_.i:
    if f_<>0 then 
      print([x_,y_],` not on f_`):
      x_:='x_':y_:='y_':
      RETURN():
    fi:
  od:
  x_:='x_':y_:='y_':
  if x_1<>x_2 then 
    y_:=y_1+(y_2-y_1)*(x_-x_1)/(x_2-x_1):
    gg_:=collect(expand(f_),x_):
    if degree(gg_,x_)<3 then 
      eprint(`Third point at infty`):
      y_:='y_':
      RETURN([]):
    fi:
    x_:=-x_1-x_2-coeff(gg_,x_,2)/coeff(gg_,x_,3):
  elif y_1<>y_2 then
    x_:=x_1+(x_2-x_1)*(y_-y_1)/(y_2-y_1):
    gg_:=collect(expand(f_),y_):
    if degree(gg_,y_)<3 then 
      eprint(`Third point at infty`):
      x_:='x_':
      RETURN([]):
    fi:
    y_:=-y_1-y_2-coeff(gg_,y_,2)/coeff(gg_,y_,3):
  else
    dx_:=diff(f_,x_):dy_:=diff(f_,y_):x_:=x_1;y_:=y_1;
    if dy_<>0 then
      m:=-dx_/dy_;x_:='x_';y_:=m*(x_-x_1)+y_1:
      gg_:=collect(expand(f_),x_):
      if degree(gg_,x_)<3 then 
        eprint(`Third point at infty`):y_:='y_':RETURN([]):
      fi:
      x_:=-2*x_1-coeff(gg_,x_,2)/coeff(gg_,x_,3):
    elif dx_<>0 then
      y_:='y_';x_:=x_1:gg_:=collect(f_,y_):
      if degree(gg_,y_)<3 then 
        eprint(`Third point at infty`):x_:='x_':RETURN([]):
      fi:
      y_:=-2*y_1-coeff(gg_,y_,2)/coeff(gg_,y_,3):
    else 
      print([x_,y_],` is a singular point`):
      x_:='x_':y_:='y_':RETURN():
    fi:
  fi:
  x_3:=x_:y_3:=y_:x_:='x_':y_:='y_':
  eprint(`Third point of intersection is `,[x_3,y_3]);
  [x_3,y_3]:
end:

trcw1:=proc() local e, h, i, j, P, Q, rr; 
  global iee, p, p_at_inf_c, p_at_inf_w;
  P:=p;
  if nargs=1 then 
    h:=nops(p_at_inf_c);
    p_at_inf_c:=[op(p_at_inf_c[1..(nW-1)]),inf_c,op(p_at_inf_c[(nW+1)..h])];
    p_at_inf_w:=[op(p_at_inf_w[1..(nW-1)]),inf_w,op(p_at_inf_w[(nW+1)..h])];
    e:=iee;iee:=0;rr:=RR;map(ap1_,inf_w);iee:=e;nowRR(rr,cur);
  fi;
  if P=3 then
    h:=0;
    for i to nops(inf_c) do
      if q=inf_c[i] then 
        Q:=inf_w[i];
        if Q=[`?`] then j:=`  --- sorry`;else j:=NULL;fi;
        if Q=[] then Q:=O;fi;
        eprint(q,` |--> `,Q,j);h:=1;
      fi;
    od;
    if h=0 then print(q,noc,qtr);p:=NULL;fi;
  else
    if nargs=0 then fprint(`Trcw([]) already done:`);fi;
    for i to nops(inf_c) do
      Q:=inf_w[i];
      if Q=[`?`] then j:=`  --- sorry`;else j:=NULL;fi;
      if Q=[] then Q:=O;fi;
      eprint(inf_c[i],` |--> `,Q,j);
    od;
    Q:=inf_w;
  fi;
  Q;
end:

Trcw:=proc() global  iee;
  iee:=1:trcw(args);iee:=0:NULL;
end:

trcw:=proc() local i, j, l, Q;
global  a, inf_c, inf_w, k, pt, s, so, U, V, Wd, Wf, x, xd, xn, y;
#  if K_=T then tesT();fi; # At present, we can have K_=T and nops(bec)>0
#  only via gcub or quar, and both of these have done tesT
  trw(args);
  if p=0 then RETURN();fi;
  if nops(inf_c)>0 and p>1 then RETURN(trcw1());
  elif p=1 and nops(b)=9 and q=b[9] then eprint(q,` |--> `,O);RETURN([]);
  fi;
  a:=collect(sex_(bec[nW]),{U,V},distributed);
  if p>1 then 
    if degree(a,U)<4 then 
      s:=0;
      for i from 0 to 3 do
        s:=s+coeff(collect(coeff(a,U,i),V),V,3-i)*U^i*V^(3-i);
      od;
      if degree(s,U)<3 then 
        inf_c:=[[1,0,0]];s:=simplify(s/V^(3-degree(s,U)));
      fi;
      if ldegree(s,U)>0 then
        inf_c:=[op(inf_c),[0,1,0]];s:=simplify(s/U^(ldegree(s,U)));
      fi;
      U:=1:so:=simplify([solve(s,V)]):U:='U':
      for i to nops(so) do inf_c:=[op(inf_c),[1,so[i],0]]:od:
      for i to nops(inf_c) do
        Wf:='Wf';Wd:='Wd';Q:=inf_c[i];
        if Q[1]=1 then U:=Wd:V:=Wf:k:=2:
        else U:=Wf:V:=Wd:k:=1:
        fi;
        so:=[solve(a,Wf)]:
        for j to nops(so) do
          Wf:=so[j]:
          if sex_(limit(Wf/Wd,Wd=infinity))=Q[k] then break:fi;
        od;
        pt:=sex_(limit(b[5]/b[6],Wd=infinity)):
        if has(pt,infinity) then pt:=[]:
        elif has(pt,undefined) or has(pt,unevaluated) then pt:=[`?`]:
        else
          pt:=[pt,sex_(limit(b[7]/b[8],Wd=infinity))];
        fi;
        inf_w:=[op(inf_w),pt];U:='U';V:='V';
      od;
      RETURN(trcw1(1));
    else
      Wd:='Wd':Wf:='Wf':
      i:='i';U:=1/i:a:=sex_(i^4*(a+V^2));j:=coeff(a,i,0);a:=a/j;
      a:=sex_(simplify(sqrt(j),sqrt,symbolic)*convert(taylor(sqrt(a),i,5),\
polynom)/i^2):
      i:=1/Wf:a:=collect(expand(a),Wf);U:=Wf:V:=Wd:l:=[];
eprint(`The two places at infinity on the quartic, both centered at the`);
eprint(`point [0,1,0], map to the following points on the Weier. curve:`);
      inf_c:=[[0,1,0],[0,1,0]];
      for Wd in [a,-a] do
        if b[6]=0 then
          if b[5]<>0 then pt:=infinity:else pt:=undefined:fi:
        else
          pt:=sex_(limit(b[5]/b[6],Wf=infinity)):
        fi;
        if has(pt,infinity) then l:=[op(l),O]:
        else 
          if b[8]=0 then
            if b[7]<>0 then pt:=infinity:else pt:=undefined:fi:
          else
            pt:=[pt,sex_(limit(b[7]/b[8],Wf=infinity))]:
          fi;
          if has (pt,infinity) then l:=[op(l),O]:
          elif has(pt,undefined) then l:=[op(l),[`?`]]:
          else l:=[op(l),pt]:
          fi:
        fi;
      od;
      U:='U';V:='V';inf_w:=l;
    fi;   
    RETURN(trcw1(1));
  else
    U:=x:V:=y:i:=simplify(a):U:='U':V:='V':
    if i<>0 then 
      print([x,y],` not on curve`,qtr):RETURN():
    fi:
    pt:=trq(b[5],b[6]):
    if pt=infinity then pt:=[]:eprint([x,y],`|----> O`):
    else 
      pt:=sex_([pt,trq(b[7],b[8])]):eprint([x,y],`|---->`,pt);Q:=pt;
      do_tor(1);ap1_(Q);pt:=Q;
    fi;
    pt;
  fi;
end:

trq:=proc(pn,pd) global  xn, xd, U, V, vn, vd, fu, fv, i;
#simplify is needed in vn:= and vd:= (also in trqw) because, e.g.,
#when m_=6 otherwise 6-w_^2 won't be recognized as 0.
  xn:=pn:xd:=pd:U:=x:V:=y:vn:=simplify(xn):vd:=simplify(xd):U:='U':V:='V':
  if vd<>0 then RETURN(vn/vd):
  elif vd=0 and vn<>0 then RETURN(infinity):
  else
    fu:=diff(a,U):fv:=diff(a,V):
    while vn=0 and vd=0 do
      i:=sex_((diff(xn,U)*fv-diff(xn,V)*fu)/(diff(xd,U)*fv-diff(xd,V)*fu)):
      xn:=numer(i):xd:=denom(i):
      U:=x:V:=y:vn:=simplify(xn):vd:=simplify(xd):U:='U':V:='V':
    od;
  fi;
  if vd=0 then infinity;else vn/vd;fi;
end:      

trqw:=proc(pn,pd) global  xn, xd, X, Y, vn, vd, fx, fy, i;
;
  xn:=pn:xd:=pd:X:=x:Y:=y:vn:=simplify(xn):vd:=simplify(xd):X:='X':Y:='Y':
  if vd<>0 then RETURN(vn/vd):
  elif vd=0 and vn<>0 then RETURN(infinity):
  else
    fx:=diff(a,X):fy:=diff(a,Y):
    while vn=0 and vd=0 do
      i:=sex_((diff(xn,X)*fy-diff(xn,Y)*fx)/(diff(xd,X)*fy-diff(xd,Y)*fx)):
      xn:=numer(i):xd:=denom(i):
      X:=x:Y:=y:vn:=simplify(xn):vd:=simplify(xd):X:='X':Y:='Y':
    od;
  fi;
  if vd=0 then infinity;else vn/vd;fi;
end:      

trw:=proc() 
global  b, inf_c, inf_w, l, nW, p, p_at_inf_c, p_at_inf_w, q, r, x, y;
  p:=0:nW:=1:
  if nops(bec)=0 then
    print(`There are no curves listed in bec`.qtr);RETURN();
  fi;
  if nargs>0 then
    q:=args[1]:
    if type(q,list) then
      if nops(q)=0 then p:=2:
      else 
        x:=q[1]:y:=q[2]:
        if nops(q)=3 then
          if q[3]=0 then 
            p:=3:
            if q[1]=0 then q:=[0,1,0];
            else q:=[1,q[2]/q[1],0];
            fi;
          elif member(q,[PP]) then p:=1;
          fi;
        else p:=1:
        fi;
      fi;
      if nargs=2 then nW:=args[2]:
      elif nops(bec)>1 then trx();
      fi:
    elif nargs=1 then
      if type(q,posint) then nW:=q:
      else print(qts);RETURN();
      fi;
    else
      x:=q:y:=args[2]:
      if nargs=3 then
        if args[3]<>0 then nW:=args[3]:p:=1:else p:=3:fi;
      elif nargs=4 then 
        if args[3]=0 then nW:=args[4]:p:=3:else print(qts);RETURN();fi;
      else 
        p:=1:
        if nops(bec)>1 then trx();fi;
      fi;
    fi;
  fi;
  if nW>nops(bec) then
    print(`There are only`,nops(bec),`curves listed in bec`,qtr);
    p:=0:RETURN():
  fi;
  b:=bic[nW]:inf_c:=p_at_inf_c[nW];inf_w:=p_at_inf_w[nW];
  if p=0 then
    if iee=0 then RETURN(b);fi;
    print(`curve is`,sort(bec[nW],[U,V])=0);
    l:=sort(Y^2+a1*X*Y+a3*Y,[Y,X]):r:=sort(X^3+a2*X^2+a4*X+a6,X);
    print(`with Weierstrass form`,l=r);
    lprint(`and the transformations are`);
    print(U=b[1]/b[2]);print(V=b[3]/b[4]);
    print(X=b[5]/b[6]);print(Y=b[7]/b[8]);
  fi;
end:

Trwc:=proc() global  iee;
  iee:=1:trwc(args);iee:=0:NULL;
end:

trwc:=proc() local  e, j, xx, yy;
 global  a, c, i, iee, p, pt, q, trou, x, X, y, Y;
  if K_=T then tesT();fi;
  trw(args);
  if p=3 then 
    if x=0 and y<>0 then 
      p:=2:
    else print(`point not on E`,qtr);p:=0:
    fi;
  fi;
  if p=0 then RETURN();fi;
  c:=sex_(b[1..4]):
  if p=2 then 
    if nops(b)=9 then eprint(`O |---->`,b[9]):RETURN(b[9]):fi;
    Y:=[solve(Y^2+a1*X*Y+a3*Y-X^3-a2*X^2-a4*X-a6,Y)]:Y:=Y[1]:
    pt:=sex_([limit(c[1]/c[2],X=infinity),limit(c[3]/c[4],X=infinity)]):
    if has(pt,infinity) then 
      q:=b[1]*b[4]/(b[2]*b[3]):pt:=sex_(limit(q,X=infinity)):
      if has(pt,infinity) then pt:=sex_(limit(1/q,X=infinity)):pt:=[1,pt,0]:
      else pt:=[pt,1,0];
      fi;
      eprint(`Image is at infinity:`);
    fi;
    Y:='Y':eprint(`O |---->`,pt):RETURN(pt):
  fi;
  if not on(x,y) then print([x,y],noc,qtr);RETURN();fi;
  if inf_c=[] then
    e:=iee:iee:=0;xx:=x;yy:=y;trcw([],nW);iee:=e;x:=xx;y:=yy;
  fi;
  if member([x,y],inf_w,j) then
    pt:=inf_c[j];eprint([x,y],`|---->`,pt,`  (at infinity)`):RETURN(pt):
  fi;
  a:=sex_(Y^2+a1*X*Y+a3*Y-(X^3+a2*X^2+a4*X+a6));
  pt:=trqw(c[1],c[2]);j:=trqw(c[3],c[4]);
  if pt=infinity or j=infinity then 
    lprint([x,y],`|----> point at infinity`);
    lprint(`but why wasn't this detected earlier in inf_c ----> inf_w`);
    lprint(`where inf_c,inf_w`-inf_c,inf_w);
    RETURN();
  fi;
  pt:=normal(sex_([pt,trqw(c[3],c[4])]));
  if has(pt,infinity) then trou:=[op(trou),TRWC];RETURN(trouble);fi;
  eprint([x,y],`|---->`,pt);pt;
end:

trx:=proc();
  if iee=1 then
    lprint(nops(bec),` curves birationally equivalent to`,cur,` are listed in b\
ec,`);
    lprint(`yet no code n was passed as a parameter to Trcw in the form Trcw(\
z,n).`);
    print(`The first one (n=1) is taken by default.`);
  fi;
end:

Weier_fit:=proc() global  iee;
  iee:=1;weier_fit(args);iee:=0:NULL;
end:

weier_fit:=proc() 
global  gg_, x_, y_, f_, a_5;
local e,i,N,p,sols:
  if nargs=0 then print(`Weier_fit requires args.`):RETURN():fi:
  N:=nargs:gg_:=args[nargs]:wprep_():
  if type(gg_,{equation,list(equation),set(equation)}) then
     N:=N-1:e:=traperror(assign(gg_)):
     if e=lasterror then 
       wprep_();ERROR(`PLEASE REPEAT THE SAME weier_fit COMMAND.`);
     fi;
  fi:
  if N>0 then
    gg_:={}:
    for i to N do
      p:=args[i]:
      if not type(p,list) or nops(p)<>2 then print(qts):RETURN():fi:
      x_:=p[1]:y_:=p[2]:
      gg_:=gg_ union {f_}:
    od:
    sols:={solve(gg_)}:
    if sols={} then print(`No solution.`):f_:=NULL;RETURN():fi:
    assign(sols):
  fi;
  x_:='x_':y_:='y_':a_5:=a_6;
  if iee=1 then
    lprint(`The general Weierstrass satisfying the requirements is`);
    print(f_);
lprint(`After specializing the remaining a_.i (if any) to rational values,`);
lprint(`one can call Ein(a_.(1..5)) (note a_5:=a_6 is automatic).`);
  fi;
  f_;
end:

wprep_:=proc()
global  a_1, a_2, a_3, a_4, a_5, a_6, x_, y_, f_;
;
  a_1:='a_1';
  a_2:='a_2';
  a_3:='a_3';
  a_4:='a_4';
  a_5:='a_5';
  a_6:='a_6';  
  x_:='x_';y_:='y_';
  f_:=x_^3+a_2*x_^2+a_4*x_+a_6-y_^2-a_1*x_*y_-a_3*y_;
end:

Zp_sol:=proc() local lst_,k_,lem_,lsu_,x_,t_,t1,t2,x,code;global p;
#This implements Cremona's alg. for p=2 and Siksek's alg. for
#p>2. See the #-lines at the top of this file.
  p:=args:
  if p>2 then RETURN(ZPSOLV());fi;
  lst_:=[0];k_:=0;
  while nops(lst_)>0 do
    k_:=k_+1:lsu_:=NULL;
    for x_ in lst_ do
      for t_ from 0 to 1 do
        x:=x_+t_*2^(k_-1);
        code:=lemma7(x,k_);
        if code=true then RETURN(1);
        elif code=FAIL then lsu_:=lsu_,x;
        fi;
      od;
    od;
    lst_:=[lsu_];
  od;
  0;
end:  

ZPSOLV:=proc() local a,c,d,i,f_,f,f1,f2,ff,epi,g,h,r,s;
#This is Siksek's fast alg., called by Zp_solv, to test for solvability 
#in Q_p of y^2=g_=a_X^4+b_X^3+c_X^2+d_X+e_; p and g_ are preset.
  f_:={g_};
  while nops(f_)>0 do
    f:=f_[1];f_:=f_ minus {f};
    while f mod p^2=0 do f:=f/p^2;od;
    ff:=f mod p;
    if ff<>0 then
      d:=degree(ff);
      if d=0 then
        if L(ff,p)=1 then RETURN(1);else next;fi;
      fi;
      if d=1 or d=3 then RETURN(1):fi;
      a:=coeff(ff,X,d);
      if L(a,p)=1 then RETURN(1);fi;
      ff:=ff/a mod p;
      for i from 0 to d-1 do c[i]:=coeff(ff,X,i);od;
      if d=2 then 
        r:=c[1]/2 mod p;
        if modp(r^2-c[0],p)>0 then RETURN(1);fi;
        g:=X+r:epi:=[-r];
      elif c[3]=0 then
        if c[1]<>0 then RETURN(1);fi;
        s:=c[2]/2 mod p;
        if modp(s^2-c[0],p)>0 then RETURN(1);fi;
        if L(-s,p)=-1 then next;fi;
        d:=msqrt(-s,p);epi:=[d,-d];g:=X^2+s;
      else
        r:=c[3]/2 mod p;s:=(c[2]-r^2)/2 mod p;
        if modp(c[1]-2*r*s,p)>0 or modp(c[0]-s^2,p)>0 then RETURN(1);fi;
        g:=X^2+r*X+s;d:=r^2-4*s;
        if modp(d,p)=0 then epi:=[-r/2 mod p];
        elif L(d,p)=-1 then RETURN(0);
        else d:=msqrt(d,p);epi:=[(-r+d)/2,(-r-d)/2] mod p;
        fi;
      fi; 
      h:=collect(expand((f-a*g^2)/p),X);
      for i in epi do
        if modp(subs(X=i,h),p)=0 then
          f_:=f_ union {collect(expand(subs(X=p*X+i,f)/p^2),X)};
        fi;
      od;
    else #this is Siksek's step II
      f1:=f/p;f2:=Factor(f1) mod p;
      for i in factors(f2)[2] do
        if degree(i[1],X)=1 then
          if i[2]=1 then RETURN(1);
          else f_:=f_ union {collect(expand(subs(X=p*X-coeff(i[1],X,0),f1)/p\
),X)};
          fi;
        fi;
      od;
    fi;
  od;
  0;
end: