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

lprint(`Reading E mod p procedures ..`); 

Allp:=proc(p1) 
global  iee;
local q,Q;
  iee:=1:allp(p1);iee:=0:q:=` on `.cur.` mod `.p:
  if allq<5 then RETURN();fi;
  if xs>=0 then
    lprint([xs,ys],` is a `,qs,q);Q:=` non-singular`:
  else
    Q:=NULL:
  fi;
  print(cat(`group of`,Q,` points`,q,` = O`),op(lisp));
  if Np=Np1 then q:=cyclic:
  else q:=`non-cyclic `.Np1.X.(Np/Np1):
  fi;
  lprint(`group order = `,Np,`, type: `.q);
end:

allp:=proc(p1) local b,i,j;
global allq, lisp, lisx, ml, mula, myl, Np, Np1, ord, p, pk, po, p_2, 
qs, q2, r, r2, x, xs, x1, x2, x4, y, ys, y1, y2;
#Find all points and the group structure of E mod p.
  if not type(p1,integer) or p1<=0 or not isprim(p1) then   
    eprint(`args must be a prime`);allq:=0:Np:=0:RETURN();
  elif cur='cur' then aM();allq:=0:Np:=0:RETURN();
  fi;
  if p1>100 and tFr3[1]=101 then frobr():fi:
  alP(p1);
  if allq<5 then Np:=0:RETURN();fi;
  lisp:={}:lisx:={}:x4:=0:p:=p1:Np1:=1:po:=p-1:r:=1:p_2:=po/2:
  if p<tFr3[1] then b:=p-emoe(p):else b:=p+trunc(evalf(2*sqrt(p)))-1:
  fi:
  for pk while modp(po,2)=0 do po:=po/2:od:
  if pk=2 then r:=2:
  elif pk>2 then
    for r from 3 while modp(power(r,p_2),p)=1 do od:
  fi;
  r2:=modp(power(r,po),p):q2:=(po+1)/2:
  if modp(DD,p)=0 then
    if modp(c4,p)=0 then
      qs:=cusp:
      if p=2 then xs:=modp(a4,2):ys:=modp(a2*a4+a6,2):
      elif p=3 then xs:=modp(-a3^2-a6,3):ys:=modp(a1*xs+a3,3):
      else xs:=modp(-b2/12,p):ys:=modp(-(a1*xs+a3)/2,p):
      fi;
    else
      qs:=node:xs:=modp((18*b6-b2*b4)/c4,p):
      if p>2 then ys:=modp(-(a1*xs+a3)/2,p):
      else ys:=modp(a3+a4,2):
      fi;
    fi;
  else
    xs:=-1:
  fi;
  while x4<p do
    nexp(x4);
    if x<p and x<>xs then
      x1:=modp(x,p):y1:=modp(y,p):x2:=x1:y2:=y1:ml:=[[x1,y1]]:myl:={x1}:
      for i do
        eadp(x1,y1,x2,y2);
        if x3=NULL then ord:=i+1:break;fi;
        x2:=x3:y2:=y3:ml:=[op(ml),[x2,y2]]:myl:=myl union {x2}:
      od;
      mula:={}:
      for j to i do
        x:=ml[j][1]:y:=ml[j][2]:mula:=mula union {[x,y,ord/igcd(j,ord)]}:
      od;
      lisp:=lisp union mula:Np1:=max(Np1,ord):lisx:=lisx union myl:
      if nops(lisp)=b then break;fi;
    fi;
    x4:=x+1:
  od;
  Np:=nops(lisp)+1:
end:

alP:=proc(palP)
global  allq, allP;
;
  allq:=5:
  if cur='cur' and not member(palP,allP) then 
    if traperror(we() mod palP)=lasterror then 
      allq:=0:eprint(`Curve not defined mod`,palP);
    else 
      allP:=allP union {palP}:
    fi;
  fi;
end:

Eadp:=proc() 
global  x3;
local q,z;
  eadp(args);q:=` on `.cur.` mod `.p:
  if allq<5 then RETURN();fi;
  if x3<>NULL then
    if x3=-1 or x3=-2 then 
      x3:=-x3:z:=[x.x3,y.x3]:print(z,` not`.q);RETURN();
    else
      z:=[x3,y3]:
    fi;
  else
    z:=O:
  fi;
  lprint(q,[x1,y1],`+`,[x2,y2],` = `,z);
end:

eadp:=proc()
global  x3, x2, y2, x1, y, pcn, z3, y3;
;
#add two points on elliptic curve mod p
  orgs(args);
  if not onp([x2,y2]) then x3:=-2:RETURN();
  else x2:=xonp:y2:=yonp:
  fi:
  if not onp([x1,y1]) then x3:=-1:RETURN();
  else x1:=xonp:y:=yonp:
  fi;
  if x1-x2 mod p=0 and sex_(y1+y2+a1*x1+a3) mod p=0 then 
    x3:=NULL:y3:=NULL;RETURN([]);
  fi;
  pcn:=1:z3:=eadd(args);pcn:=0:
  if x3=NULL or traperror(z3 mod p)=lasterror then 
    x3:=NULL:y3:=NULL:RETURN([]);
  fi;
  z3:=Normal(z3) mod p:x3:=z3[1]:y3:=z3[2]:z3;
end:

Mulp:=proc();
  print(mulp(args));
end:

mulp:=proc() 
global  zm, xm, ym, n, Nz;
local u,v;
#calc. m*z in E reduced mod p where args = m,z and p is preset
  if not type(p,integer) or p<=0 or not isprim(p) then   
    print(`p must be preset to a prime value`);RETURN();
  fi;
  if nargs=3 then zm:=[args[2],args[3]]:
  elif nargs=2 and type(args[2],list) then zm:=args[2]:
  else print(`Improper args in mulp`,qtr);RETURN();
  fi;
  xm:=NULL:ym:=NULL:n:=args[1]:
  if n=0 or zm=[] then RETURN([]):fi;
  if not onp(zm) then print(zm,noc.` mod `.p);RETURN();fi;
  if xonp=NULL then RETURN([]):fi;
  u:=xonp:v:=yonp:
  if nops(zm)=3 then n:=mods(n,zm[3]):fi;
  if n<0 then n:=-n:v:=sex_(-v-a1*u-a3):fi;  
  while n>=1 and u<>NULL do
    if modp(n,2)=1 then
      if xm=NULL then 
        xm:=u:ym:=v:
      else
        eadp(xm,ym,u,v);xm:=x3:ym:=y3:
      fi;
    fi;
    n:=trunc(n/2):eadp(u,v,u,v);u:=x3:v:=y3:
  od;
  Nz:=[xm,ym]:
end:

Negp:=proc()
global  Edd;
;
  Edd:=0:negp(args);
  if Edd=0 then lprint(`-`,pz,` = `,nz,` mod `.p);fi;
  NULL;
end:

negp:=proc()
global  pz, nz, Edd;
;
  if nargs>1 then pz:=[args[1..nargs]]:
  else pz:=args:
  fi;
  if pz=[] then nz:=[]:RETURN([]):fi;
  if not onp(pz) then Edd:=1:print(qts);RETURN();fi;
  if xonp=NULL then pz:=[]:nz:=[]:RETURN([]):fi;
  nz:=[xonp,sex_(-yonp-a1*xonp-a3 mod p),op(pz[3..nops(pz)])]:
end:

nexp:=proc(x5) 
global  b3, b5, b7, x, y;
local f,h,q;;
  if p>3 then 
    h:=modp(-p_2,p):q:=modp(h^2,p):b3:=modp(b2*q,p):b5:=modp(b4*h,p):
    b7:=modp(b6*q,p):
  fi;
  for x from x5 to p-1 do
    if member(x,lisx) then next;fi;
    if p<5 then
      for y from 0 to p-1 do
        if modp(y*(y+a1*x+a3)-x*(x*(x+a2)+a4)-a6,p)=0 then RETURN();fi;
      od;
    else
      f:=modp(x*(x*(x+b3)+b5)+b7,p):
      if f=0 or L(f,p)=1 then
        y:=modp(power(f,q2),p):
        while modp(y^2-f,p)>0 do
          y:=y*r2:
        od;
        y:=y-(a1*x+a3)*h:RETURN();
      fi;
    fi;
  od;
end:

Onp:=proc();
  print(onp(args));
end:

onp:=proc()
global  zonp, xonp, yonp;
;
  alP(p);
  if allq<5 then RETURN(false);fi;
  if not type(p,integer) or p<=0 or not isprim(p) then   
    print(`p must be preset to a prime value`);RETURN(false);
  fi;
  if nargs=1 then zonp:=args:
  else zonp:=[args]:
  fi;
  if args=[] or traperror(zonp mod p)=lasterror then 
    xonp:=NULL:yonp:=NULL:RETURN(true):
  fi;
  xonp:=sex_(zonp[1]) mod p:yonp:=sex_(zonp[2]) mod p:
  evalb(sex_(yonp*(yonp+a1*xonp+a3)-xonp*(xonp*(xonp+a2)+a4)-a6\
) mod p=0):
end: