Fact:=proc(a);
  lprint(`The factorization of`,a,` in Z[w_]`);
  lprint(`up to a unit factor is`);
  print(fact(a));NULL
end:

fact:=proc(a) local a_,n,pp_;global facs_,iee,p,p_;
#Factorize a (an element of the ring of integers O of Q(sqrt(m_)) in O
#and let facs_ be the list of prime power divisors p_^{v_} in the
#format [p_,v_], where p_ is either a prime element, or, possibly when 
#the class number h_<>1, a prime ideal in the format P.p or Q.p. 
#(When p is ramified only P.p appears; if p is split then both
#P.p and Q.p appear with v_ = (resp.) N and zz-N where N denotes an
#unknown integer while zz has a particular value depending on a.)
  if m_>0 then n:=iee:iee:=0:fundunit();iee:=n;fi;
  facs_:=[];a_:=qf_conv(a);
  if not type(a_,polynom(integer,w_)) then
    print(`args for fact must be an alg. integer in Q(sqrt(m_))`,qtr);
    RETURN();
  fi;
  n:=nm_(a_):
  for pp_ in factorset(n) minus{1,-1} do
    for p_ in pri_(pp_) do
      if type(p_,algnum) then
        val_(a_);
        if v_>0 then facs_:=[op(facs_),[p_,v_]];fi;
      else
        p:=pp_:vam(n);
        if member(p,ram_) then facs_:=[op(facs_),[P_.p,v]];
        elif p_=P_.p then facs_:=[op(facs_),[p_,N]]
        else facs_:=[op(facs_),[Q_.p,v-N]];
        fi;
      fi;
    od;
  od;
  facs_
end:    

Fundunit:=proc() global iee;
  iee:=1:fundunit():iee:=0:
  if m_>0 then
    lprint(`(Apecs notation: the fundamental unit, i.e., the smallest unit >1\
,`);
    lprint(` = unit_ + vnit_*w_, with norm = nnit_ = +-1, and`);
    lprint(` absolute value (to Digits accuracy) enit_.)`);
  fi;
  NULL
end:

fundunit:=proc() local A,b,d,p,q,t,u1,u2,v1,v2;
  global enit_,nnit_,unit_,vnit_;
#Taken from Cohen's book, p.265; fund. unit = unit_+vnit_*w_ (except
#our w_ = sqrt(m_) or (-1+sqrt(m_))/2) with
#norm = nnit_ =+-1 and abs. value enit_.
  if unit_=0 then 
    if m_<0 then 
      print(`This proc. is for REAL quad. fields only.`);RETURN();
    fi;
    d:=floor(sqrt(D_)):
    if modp(d-D_,2)=0 then b:=d;else b:=d-1:fi:
    u1:=-b:u2:=2:v1:=1:v2:=0:p:=b:q:=2:
    do
      A:=floor((p+d)/q):t:=p:p:=A*q-p:
      if t<>p or v2=0 then
        t:=A*u2+u1:u1:=u2:u2:=t:t:=A*v2+v1:v1:=v2:v2:=t:t:=q:q:=(D_-p^2)/q:
        if q=t and v2<>0 then 
          unit_:=abs((u1*u2+D_*v1*v2)/(2*q)):vnit_:=abs((u1*v2+u2*v1)/q):
          break;
        fi;
      else
        unit_:=abs((u2^2+D_*v2^2)/(2*q)):vnit_:=abs((2*u2*v2)/q):
        break;
      fi;
    od;
    if ch_=4 then enit_:=evalf(unit_+vnit_*sqrt(m_)/2):unit_:=unit_+vnit_/2:
    else enit_:=evalf(unit_+vnit_*sqrt(m_)):  
    fi;
    nnit_:=nm_(unit_+vnit_*w_):
  fi;
  eprint(`Fundamental unit`=unit_+vnit_*'w_');
  fprint(`with norm`=nnit_,` and abs. value`=enit_);enit_
end:

L2_:=proc() local i,j;global a_,c4__,c6__,e_,Psi2_,Psi3_,xL;
  if (a_=0 and type(x2_[mods(-c6__,4)],list)) or \
   (a_>=4*v2_ and type(x2_[mods(c6__/8,4)],list)) then 
#OK when c4=0 and a_ has the fictitious value c_-4*min_ since min_>=1
    RETURN();
  fi;
  if v2_=2 then
    Psi2_:=xL^3-3*c4__*xL-2*c6__;
    Psi3_:=xL^4-6*c4__*xL^2-8*c6__*xL-3*c4__^2;
    for i from 0 to 1 do
      for j from 0 to 1 do
        xL:=(i+j*w_)^2;
        if val_(Psi3_)>15 and val_(Psi2_)>7 and type(x2_[mods(simplify\
(expand(-Psi2_/16)),4)],list) then
          RETURN();
        fi;
      od;
    od;
  fi;
  a_:=a_+4:e_:=e_-1;
  c4__:=simplify(c4__*p_^4);c6__:=simplify(c6__*p_^6);
  L2_()
end:
  
L3_:=proc() local i,j;global b_,c4__,c6__,e_,Psi2_,xL;
  b_:=b_-6*min_;
  if not member(3,ram_) then
    if b_=2 then e_:=e_-1;fi;
  elif a_>0 and b_<6 then
    Psi2_:=xL^3-3*c4__*xL-2*c6__;   
    for i from -1 to 1 do
      for j from -1 to 1 do
        xL:=i+j*w_;
        if val_(Psi2_)>5 then RETURN();fi;
      od;
    od;
    e_:=e_-1;
  fi
end:

Minf:=proc() global iee;
  iee:=1;minf();iee:=0;NULL
end: 

minf:=proc() 
  local aa_,bb_,cc_,i,j,old_PP,old_RR,p2_,p3_,p4_,pp_,x,y,y2,z;
  global a_,ap1,ap2,ap3,ap4,ap5,ap6,bec,bic,b_,bp2,bp4,bp6,c_,c4_,c4__,c6_,
  c6__,cp4,cp6,DD_,e_,iee,min_,new_,NN,old_,old_NN,p,PP,p_,p_at_inf_c,
  p_at_inf_w,RR,r_,s_,t_,tau_,U,u_,V,X,x2_,x3_,x4_,Y;
#After ell. curve E defined by Weier. equation over ring of integers O of
#Q(sqrt(m_)) is entered using ell --- the Weier. coeff.'s must be expressed
#in terms of w_ (see #remarks for qfin) ---  find global minimal Weier. 
#equation and call ell with new Weier. coeff.'s (which are stored in new_; 
#the original coeff.'s are in old_). tau_=[r_,s_,t_,u_] is the 
#transformation from the old to the new E.
  if m_<>K_ then aa_:=iee:iee:=0:ell(cur,m_);iee:=aa_:fi;
  if h_>1 then 
    eprint(`Not yet programmed for class number >1 --- sorry.`);RETURN();
  elif h_=0 then
    fprint(`apecs is not certain that the class number is 1;`);
    fprint(`if you are not certain you proceed at your own risk.`);
  fi;
  if x4_='x4_' then
    if m_>0 and unit_=0 then fundunit();fi;
    for i in ram_ do pri_(i);od;
    p4_:=[];x4_:='x4_';p2_:=[];x2_:='x2_';
    for i from 0 to 1 do
      for j from 0 to 1 do
        x:=i+j*w_:y2:=simplify(x^2);
        y:=simplify(y2^2);y2:=muds(y2,4);y:=muds(y,8);
        if not member(y,p4_) then
          p4_:=[op(p4_),y];x4_[y]:=[x];
        else x4_[y]:=[op(x4_[y]),x];
        fi;
        if not member(y2,p2_) then
          p2_:=[op(p2_),y2];x2_[y2]:=[x];
        fi;
      od;
    od;
    p3_:=[];x3_:='x3_';
    for i from -1 to 1 do
      for j from -1 to 1 do
        x:=i+j*w_:y:=simplify(x^3);y:=muds(y,3);
        if not member(y,p3_) then
          p3_:=[op(p3_),y];x3_[y]:=[x];
        else x3_[y]:=[op(x3_[y]),x];
        fi;
      od;
    od;
  fi;
  u_:=1;DD_:=simplify(DD);c4_:=simplify(c4);c6_:=simplify(c6);
  aa_:=abs(nm_(DD_));
  for pp_ in factorset(aa_) do
    p:=pp_;
    if val(aa_)<12 then next;fi;
    pri_(pp_);
    for p_ in pfac_[pp_] do
      c_:=val_(DD_);
      if c_<12 then next;fi;
      if c4_=0 then a_:=c_;else a_:=val_(c4_);fi;
      if c6_=0 then b_:=c_;else b_:=val_(c6_);fi;
      min_:=min(floor(a_/4),floor(b_/6),floor(c_/12));
      if min_=0 then next;
      elif pp_>=5 then e_:=min_;
      else 
        c4__:=simplify(c4_*p_^(-4*min_)):c6__:=simplify(c6_*p_^(-6*min_));
        e_:=min_;a_:=a_-4*min_;
        if pp_=3 then L3_();
        else L2_();
        fi;
      fi;
      u_:=u_*p_^e_;
    od;
  od;
  u_:=simplify(u_);
  cp4:=simplify(c4_*u_^(-4));cp6:=simplify(c6_*u_^(-6));
  aa_:=x4_[muds(cp4,8)];
  if nops(aa_)=0 then print(`Trouble in main_ at A, cp4 =`,cp4);RETURN();fi;
  for ap1 in aa_ do
    s_:=simplify((u_*ap1-a1)/2);
    if not type(s_,polynom(integer,w_)) then next;fi;
    bb_:=x3_[muds(simplify(-cp6-ap1^6),3)];
    if nops(bb_)=0 then next;fi;
    for ap2 in bb_ do
      r_:=simplify((u_^2*ap2-a2+s_*a1+s_^2)/3);
      if not type(r_,polynom(integer,w_)) then next;fi;
      bp2:=simplify(ap1^2+4*ap2);
      bp4:=simplify((bp2^2-cp4)/24);
      bp6:=simplify((-bp2^3+36*bp2*bp4-cp6)/216);
      if not type(bp4,polynom(integer,w_)) or not \
type(bp6,polynom(integer,w_)) then
        next;
      fi;
      cc_:=x2_[muds(bp6,4)];
      if nops(cc_)=0 then next; else ap3:=cc_[1];fi;
      t_:=simplify((u_^3*ap3-a3-r_*a1)/2);
      ap4:=simplify((bp4-ap1*ap3)/2);
      if not type(t_,polynom(integer,w_)) or not \
type(ap4,polynom(integer,w_)) then
        next;
      fi;
      ap6:=simplify((bp6-ap3^2)/4);
      if type(ap6,polynom(integer,w_)) then 
        ap5:=ap6:tau_:=[r_,s_,t_,u_];old_:=[a.(1..5)];new_:=[ap.(1..5)];
        if old_=new_ then
          fprint(` `);eprint(`No change -- E is already minimal.`);
        else
          if NN>0 then old_PP:=PP;old_NN:=NN;fi;
#We don't want ell below to print "Tor not done..." when NN>0
          old_RR:=RR;fprint(` `);
          fprint(`Data for the new min. Weierstrass equation :--`);
          ell(new_,m_);U:='U';V:='V';X:='X';Y:='Y';
bec:=[op(bec),collect(U^3+old_[2]*U^2+old_[4]*U+old_[5]-V^2-old_[1]*U*V-\
old_[3]*V,{U,V},distributed)];
bic:=[op(bic),[op(collect(simplify([u_^2*X+r_,1,u_^3*Y+s_*u_^2*X+t_,1]),\
{X,Y},distributed)),op(collect(simplify([(U-r_)/u_^2,1,(V-s_*U+s_*r_-t_)/\
u_^3,1]),{U,V},distributed)),[]]];
          p_at_inf_c:=[op(p_at_inf_c),[[]]];
          p_at_inf_w:=[op(p_at_inf_w),[[]]];
          RR:=map(qftr,old_RR);
          fprint(`The coeff.'s of the "old" (original) weier. equation are`);
          fprint();eprint(` old_`=old_);
          fprint(`Those of the min. weier. "new" equation are`);
          fprint();eprint(` new_`=new_);
          i:=nops(bec);
          if i=1 then j:=z;else j:=z,i;fi;
          fprint(cat(\
`The commands Trcw(`,j,`) and Trwc(`,j,`) are available during this`));
fprint(`apecs session to transfer points z between the "old" and "new" cur\
ves.`);
          if NN=-1 and type(old_NN,posint) then          
            PP:=op(map(qftr,[old_PP]));NN:=old_NN;print_PP();old_NN:='old_NN';
          fi;
        fi;
        DD_:=simplify(DD);RETURN();
      fi;
    od;
  od;
  print(`An error somewhere`)
end:      

muds:=proc(z,m) mods(coeff(z,w_,0),m)+mods(coeff(z,w_,1),m)*w_ end:

pri_:=proc(p) local n; global pfac_, SSo;
#ideal factor the prime p in Q(sqrt(m_) 
  if type(pfac_[p],list) then RETURN(pfac_[p]);fi;
  if member(p,ram_) or L(m_,p)=1 then
#For m=1 mod 4, L(m,2)=1 (resp. -1) for m=1 mod 8 (resp. m=5 mod 8)
#in Maple5 Release 5, correcting Release 3.
    SSo:=(ram_ union {p}) minus {2};
    if nrs(m_,p)=-1 and nrs(m_,-p)=-1 then n:=[0,0]:
    elif m_>0 then
      n:=uv_(ceil(-enit_/sD_),(p-1)/sD_,p);
#search bounds passed to uv_ from Hasse, Number Theory, p.578
      if n[1]=0 and nnit_=1 then n:=uv_(1,(p+enit_)/sD_,-p);fi;
    else n:=uv_(1,2*p/sD_,p):
    fi;
    if n[1]=0 then 
      if member(p,ram_) then pfac_[p]:=[P_.p];
      else pfac_[p]:=[P_.p,Q_.p];
      fi;
    else pfac_[p]:=n;
    fi;
  else pfac_[p]:=[p]:
  fi
end:

Qfin:=proc(m::integer) global iee;
  iee:=1;qfin(m);iee:=0;NULL
end:

qfin:=proc(m::integer) local p_;
  global ch_,D_,enit_,h_,m_,m1_,nm_,nnit_,pfac_,pol_,ram_,
sD_,Tabqf_,unit_,vnit_,v2_,w_,x2_,x3_,x4_;
#set m_:=m and initialize apecs to treat elliptic curves defined over 
#Q(sqrt(m_)). (Many of the names used in this package end in an underline _ 
#to make them distinctive). An integral basis is 1,w_ where w_= sqrt(m_)
#or (-1+sqrt(m_))/2. w_ is a root of pol_.  ram_ is the set of ramified 
#primes. Note: except when m_=-1, alpha is defined as an alias for w_, 
#and appears in print outs as the Greek letter, and in lprint-outs as 
#the word alpha. When m_=-1, w_ is printed and lprinted as I.
  if m_=m then 
    fprint(`Apecs is already pointing to Q(sqrt(m_)) where m_`=m_);
    RETURN();
  fi;
  if not issqrfree(m) or m=1 then eprint(qts);RETURN();fi;
  if type(m_,integer) then
    pfac.m1_:=copy(pfac_);
    if m_>0 then
      Tabqf_[m_,14]:=unit_; 
      Tabqf_[m_,16]:=enit_;
      if unit_>0 then
        Tabqf_[m_,15]:=vnit_; 
        Tabqf_[m_,17]:=nnit_;
      fi;
    fi;
    if not type(Tabqf_[m_,4],posint) then 
      Tabqf_[m_,1]:=w_;
      Tabqf_[m_,2]:=pol_;
      Tabqf_[m_,3]:=ram_;
      Tabqf_[m_,4]:=ch_;
      Tabqf_[m_,5]:=h_;
      Tabqf_[m_,6]:=m1_;
      x4.m1_:=copy(x4_);
      x2.m1_:=copy(x2_);
      x3.m1_:=copy(x3_);
      Tabqf_[m_,13]:=v2_;
    fi;
  fi;
  m_:=m;ch_:=Tabqf_[m_,4]:
  if type(ch_,posint) then
    D_:=m_*4/ch_:sD_:=evalf(sqrt(abs(D_))):
    w_:=Tabqf_[m_,1]:
    pol_:=Tabqf_[m_,2]:
    ram_:=Tabqf_[m_,3]:
    h_:=Tabqf_[m_,5]:
    m1_:=Tabqf_[m_,6];
    pfac_:=copy(pfac.m1_):
    if ch_=4 then
      nm_:=proc(a) local q,s,t;
        q:=simplify(a);
        s:=coeff(q,w_,0);t:=coeff(q,w_,1);s^2-s*t-(m_-1)*t^2/4:
      end:
    else 
      nm_:=proc(a) local q;
        q:=simplify(a);coeff(q,w_,0)^2-m_*coeff(q,w_,1)^2:
      end:
    fi;
    x4_:=copy(x4.m1_):
    x2_:=copy(x2.m1_):
    x3_:=copy(x3.m1_):
    v2_:=Tabqf_[m_,13]:
    if m_>0 then
      unit_:=Tabqf_[m_,14]:
      enit_:=Tabqf_[m_,16]:
      if unit_>0 then
        vnit_:=Tabqf_[m_,15]:
        nnit_:=Tabqf_[m_,17]:
      fi;
    fi;
    fprint(`Apecs returns to the quadratic field Q(sqrt(m_)) where m_`=m_);
    if m_<>-1 then alias(alpha=w_);fi;
    RETURN(m_);
  fi;   
  m1_:=2*abs(m_)+(1-sign(m_))/2:
  if m<0 then 
    if member(m,{-1,-2,-3,-7,-11,-19,-43,-67,-163}) then h_:=1:
    else h_:=2:
    fi;
  else
    if m>500 then h_:=0:
    else h_:=1:
    fi;
    if (m>2 and modp(m,8)=2) or (modp(m,2)=0 and nops(factorset(m/2))>1)\
    or member(m,{79,142,223,229,254,257,321,326,359,401,439,443,469,473}) then 
      h_:=2:
    elif modp(m,2)=1 then
      w_:=factorset(m);
      if nops(w_)>2 then h_:=2:
      elif nops(w_)=2 then
        for p_ in w_ do
          if modp(p_,4)=1 then h_:=2;fi;
        od;    
      fi;
    fi;
  fi;
  if modp(m,4)=1 then 
    pol_:=_Z^2+_Z-(m-1)/4;w_:=RootOf(pol_);ram_:={};ch_:=4;
    nm_:=proc(a) local q,s,t;
      q:=simplify(a);
      s:=coeff(q,w_,0);t:=coeff(q,w_,1);s^2-s*t-(m_-1)*t^2/4:
    end:
  else 
    pol_:=_Z^2-m;ram_:={2};ch_:=1;
    if m=-1 then w_:=I;else w_:=RootOf(pol_);fi;
    nm_:=proc(a) local q;
      q:=simplify(a);coeff(q,w_,0)^2-m_*coeff(q,w_,1)^2:
    end:
  fi;
  D_:=m_*4/ch_:sD_:=evalf(sqrt(abs(D_))):x4_:='x4_';
  unit_:=0:pfac_:='pfac_';ram_:=ram_ union factorset(m) minus {1,-1};
  if iee>0 then 
    lprint(`Apecs is pointing to the quadratic field Q(sqrt(m_)) where m_`=m_);
    lprint(`with integral basis 1, w_`=w_);
    lprint(`and ramified prime(s) ram_`=ram_);
    if h_=1 then lprint(`Class number = 1`);
    elif h_=2 then lprint(`Class number > 1`);
    else lprint(`Class number unknown`);
    fi;
    if m_=-1 then
      print(`Note that w_ appears in print-outs as I.`);
    else
      print(`Note that w_ appears in print-outs as the Greek letter`,alpha);
      lprint(`and in lprint-outs as alpha.`);
    fi;
  fi;
  if m_<>-1 then alias(alpha=w_);fi;
  p_:=abs(m_);
  if modp(m_,4)=1 then pfac_[p_]:=[1+2*w_];
  else pfac_[p_]:=[w_];
  fi;
  if member(2,ram_) then v2_:=2;else v2_:=1;fi;
  m_
end:

qftr:=proc(p) local b, t; global U, V;
  b:=bic[nops(bic)];U:=p[1];V:=p[2];
  t:=simplify(eval([b[5]/b[6],b[7]/b[8]]));
  U:='U';V:='V';
  if nops(p)=3 then t:=[op(t),p[3]];fi;
  t
end:

qf_conv:=proc(a) local b;
  b:=simplify(expand(a));
  if has(b,sqrt(m_)) then 
    b:=simplify(subs(sqrt(m_)=op(pfac_[abs(m_)]),b));
  fi;
  b
end:

subw_:=proc(a,ex);
  coeff(ex,w_,0)+coeff(ex,w_,1)*a
end:  

Torq:=proc() global iee;
  iee:=1:
  if NN<0 then
    if torq(args)=0 then RETURN();fi;
  fi;
  print(cat(`The order of the torsion subgroup over Q(sqrt(`,m_,`)) is `,NN));
  if NN>1 then
    lprint(`The nonzero torsion points, in the notation [x,y,n]`);
    lprint(`where n is the order of the point, form the sequence PP,`);
    lprint(`and, in the notation [x,y], form the list ouP.`);
  fi;
  iee:=0;NULL
end:

torq:=proc() local hex,i,ie,j,P,q,r,s,t,x_,y;
  global iee,lbNN,NN,Nsubg,p,PP,ouP,ubNN,x;
  if K_<>m_ then 
    ell(cur,m_);
    if K_=0 then RETURN(0);fi;
  fi;
  ub_tor(args);
  if NN>0 then NN_pair();RETURN(NN);fi;
  t:={seq(s[1],s=ouP)};
  if modp(ubNN,2)=0 then
    x:='x';
    for r in roots(prac,w_) do
      x_:=r[1];P:=simplify([x_,-(a1*x_+a3)/2]);
      if not member(P,ouP) then
        PP:=PP,[op(P),2];ouP:=[op(ouP),P];
        if nops(ouP)+1=ubNN then NN:=ubNN;NN_pair();RETURN(NN);fi;
      fi;
    od;
  fi;
  Nsubg:=5000;
  if nops(ouP)+1>lbNN and lbNN>1 then
      subg(ouP);PP:=op(GG);ouP:=[seq([i[1],i[2]],i=GG)];
  fi;
  lbNN:=nops(ouP)+1:
  if modp(lbNN,2)=1 then p:=2:vam(ubNN):ubNN:=zz;fi;
  if lbNN=ubNN then NN:=ubNN;NN_pair();RETURN(NN);fi;
  ie:=iee;iee:=0;ofac(ubNN);
  for i in ntab do
    p:=i[1];
    if vam(ubNN/lbNN)=0 then next;fi;
    for j to i[2] do
      q:=p^j;
      if q=2 then next;fi;
      div(q);s:=roots(al.q,w_);hex:=0;
      for r in s do
        x_:=r[1];
        if member(x_,t) then next;fi;
        t:=t union {x_};
        y:=roots(Y^2+(a1*x_+a3)*Y-x_^3-a2*x_^2-a4*x_-a6,w_):
        if nops(y)=0 then next;fi;
        P:=[x_,y[1][1],q]:PP:=PP,P,neg(P);P:=[P[1],P[2]]:
        ouP:=[op(ouP),P,neg(P)]:hex:=1;
      od;
      if hex=1 then 
        subg(ouP);PP:=op(GG);ouP:=[seq([i[1],i[2]],i=GG)]:
        lbNN:=1+nops(ouP);
        if lbNN=ubNN then NN:=ubNN;iee:=ie;NN_pair();RETURN(NN);fi;        
        t:={seq(i[1],i=ouP)};
      else
        ubNN:=ubNN/p^(i[2]-j+1);
        if lbNN=ubNN then NN:=ubNN;iee:=ie;NN_pair();RETURN(NN);fi;        
        break;
      fi;
    od;
  od;
  iee:=ie;NN:=nops(ouP)+1;NN_pair();NN
end:      

Ub_tor:=proc() global iee;
  iee:=1;ub_tor(args);iee:=0;
  if NN=0 then print(`NN divides ubNN =`,ubNN);
  elif ubNN=0 then print(` NN`=NN);
  fi;
  NULL
end:  

ub_tor:=proc() local hex,i,ue,PP_,ps,q,s,x;
  global iee,lbNN,nmD_,NN,ouP,p,PP,ubNN;
  ubNN:=0;ps:=op(yo),m_;
  if NN>0 then eprint(`Torq already done:`);RETURN(NN);fi;
  nmD_:=factorset(abs(nm_(DD)));
  if type(yo,list(rational)) and member(yo,et2,i) and eT[i][3]>0 then
    PP_:=eT[i][4];lbNN:=eT[i][3];
  else PP_:=[]:lbNN:=1:
  fi;
  eprint(`The "bad" primes are`,nmD_);
  if nargs=0 then hex:=4;else hex:=args;fi;
  for p in p_up_to_100 while hex>0 do
    if member(p,nmD_) then next;fi;
    hex:=hex-1;q:=1;
    if member(p,ram_) then 
      if p=2 then s:=modp(m_,2);q:=4;
      elif modp(m_,4)=1 then s:=(p-1)/2;
      else s:=0;
      fi;
      if p=3 then q:=3;fi;
    elif p=2 or L(m_,p)=-1 then hex:=hex+1;next;
    elif modp(m_,4)=1 then
      s:=modp((-1+msqrt(m_,p))/2,p);
    else s:=msqrt(m_,p);
    fi;
    for i to 5 do
      a_.i:=subw_(s,a.i);
    od;
    ue:=iee;iee:=0;bc(a_.(1..5));emoe(p);iee:=ue;
    if ubNN=0 then ubNN:=card*q;
    else ubNN:=igcd(ubNN,card*q);
    fi;
    if iee>0 then
      if q>1 then q:=cat(` and |Ker of reduc. mod p| divides `,q);
      else q:=NULL;
      fi;
      lprint(`For p =`,p,`,  card(E(F_p)) =`,card,q);
      lprint(`so upper b. for NN is ubNN =`,ubNN);
    fi;
    ue:=iee;iee:=0;bc(ps);iee:=ue;
    if ubNN=lbNN then hex:=0;NN:=lbNN;NN_pair();fi;
  od;
  ouP:=[seq([i[1],i[2]],i=PP_)];PP:=op(PP_):ubNN
end:

uv_:=proc(L,U,n) local n4,u,u1,u2,v,x;
#called by pri_ in 3 places to find u,v with L<=v<=U such that
#D_*v^2+4*n=u^2 and output as the norm nm_(u1+v*w_)=n, or
#if no sol'n. RETURN([0])
  n4:=n*4;
  for v from L to U do
    x:=D_*v^2+n4;
    if issqr(x) then
      u:=isqrt(x)/2;
      if ch_=1 then u1:=u:u2:=-u:
      else u1:=u+v/2;u2:=v/2-u:
      fi;    
      if member(abs(n),ram_) then RETURN([u1+v*w_]);
      else RETURN([u1+v*w_,u2+v*w_]);
      fi;
    fi;
  od;
  [0]
end:

val_:=proc(a) global v_,zz_;
#assumes that p_ is a prime in the quadratic field Q(w_) where
#w_=RootOf(pol_) where pol_=_Z^2-m_ or _Z^2+_Z-(m_-1)/4
  zz_:=simplify(a);
  if zz_=0 then RETURN(infinity);fi;
  v_:=-1;
  while type(coeff(zz_,w_,0),integer) and type(coeff(zz_,w_,1),integer) do
    v_:=v_+1;zz_:=simplify(zz_/p_);
  od;
  v_
end: