## Maple V.5

EC:=proc(ec::polynom,indetlist::list(name),npt::list)
  local x,y,X,Y,T,elliptic_curve,
        normalizationproc,reductionproc,isequal,nullpt,
        thirdpoint1,thirdpoint2,
	elliptic_add,double,subt,opposite,ext1,ext2,
	i,grad_pt;
  global Nullpt,Isequal,Add,Double,Opposite,Subt,Ext1,Ext2;
	
  if nargs<3 then ERROR("EC expects to receive at least 3 arguments, but received",args) fi;
  if nops(indetlist)=2 then
    x:=indetlist[1];
    y:=indetlist[2];
    elliptic_curve:=collect(numer(subs(x=X/T,y=Y/T,ec)),{X,Y,T})
  elif nops(indetlist)=3 then
    X:=indetlist[1];
    Y:=indetlist[2];
    T:=indetlist[3];
    elliptic_curve:=collect(ec,{X,Y,T})
  else
    ERROR("EC expects its second argument to be a list of 2 or 3 names, but received",args)     
  fi;
  if degree(elliptic_curve,{X,Y,T})<>3 then
    ERROR("EC expects its first argument to define a cubic curve, but received",args)
  fi;
  nullpt:=npt;
  if nops(nullpt)=2 then nullpt:=[nullpt[1],nullpt[2],1]
  elif nops(nullpt)>3 then
    ERROR("EC expects its third argument to be a list with 2 or 3 components, but received",args)
  fi;
  if nargs>3 then normalizationproc:=eval(args[4])
  else  normalizationproc:=proc(expr) expr end
  fi;
  if nargs>4 then reductionproc:=eval(args[5])
  else  reductionproc:=proc(expr) expr end
  fi;
  
  isequal:=
  subs(_normalizationproc=normalizationproc,
  proc(pt1,pt2)
    local pt;
    _normalizationproc(pt1[2]*pt2[3]-pt1[3]*pt2[2])=0 and
    _normalizationproc(pt1[3]*pt2[1]-pt1[1]*pt2[3])=0 and
    _normalizationproc(pt1[1]*pt2[2]-pt1[2]*pt2[1])=0
  end);
  
  thirdpoint1:=
  subs(_elliptic_curve=elliptic_curve,
  proc(pt1,pt2)
    local eq,a,b;
    eq:=subs(X=lambda*pt1[1]+mu*pt2[1],
             Y=lambda*pt1[2]+mu*pt2[2],
  	     T=lambda*pt1[3]+mu*pt2[3],
  	   _elliptic_curve);
    eq:=collect(eq/lambda/mu,{lambda,mu});
    a:=coeff(eq,lambda,1);
    b:=coeff(eq,mu,1);
    subs(lambda=-b,mu=a,lambda=1,mu=1,[
            lambda*pt1[1]+mu*pt2[1],
            lambda*pt1[2]+mu*pt2[2],
  	    lambda*pt1[3]+mu*pt2[3]]);
  end);

  grad_pt:=[seq(diff(elliptic_curve,i),i=[X,Y,T])];

  thirdpoint2:=
  subs([_X=X,_Y=Y,_T=T,_elliptic_curve=elliptic_curve,_grad_pt=grad_pt],
  proc(pt1)
    local pt,eq,pt2,a,b,i,lambda,mu;
    pt:=subs([_X=pt1[1],_Y=pt1[2],_T=pt1[3]],_grad_pt);
    pt2[1]:=pt1[2]*pt[3]-pt1[3]*pt[2];
    pt2[2]:=pt1[3]*pt[1]-pt1[1]*pt[3];
    pt2[3]:=pt1[1]*pt[2]-pt1[2]*pt[1];
    eq:=subs(_X=lambda*pt1[1]+mu*pt2[1],
             _Y=lambda*pt1[2]+mu*pt2[2],
  	     _T=lambda*pt1[3]+mu*pt2[3],
  	   _elliptic_curve);
    eq:=collect(eq/mu^2,{lambda,mu}); 
    a:=coeff(eq,lambda,1);
    b:=coeff(eq,mu,1);
    subs(lambda=-b,mu=a,lambda=1,mu=1,[
            lambda*pt1[1]+mu*pt2[1],
            lambda*pt1[2]+mu*pt2[2],
  	    lambda*pt1[3]+mu*pt2[3]]);
  end);
  
  elliptic_add:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _reductionproc=reductionproc,
  proc(pt1,pt2)
    local pt;
    if _isequal(pt1,pt2) then
      pt:=_thirdpoint2(pt1)
    else
      pt:=_thirdpoint1(pt1,pt2)
    fi;
    if _isequal(pt,_nullpt) then
      _nullpt
    else
      _reductionproc(_thirdpoint1(pt,_nullpt))
    fi
  end);
    
  double:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _reductionproc=reductionproc,
  proc(pt1)
    local pt;
    pt:=_thirdpoint2(pt1);
    if _isequal(pt,_nullpt) then
      _nullpt
    else
      _reductionproc(_thirdpoint1(pt,_nullpt))
    fi
  end);
  
  subt:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _reductionproc=reductionproc,
  proc(pt1,pt2)
    local pt;
    if _isequal(pt1,_nullpt) then
      pt:=_nullpt
    else
      pt:=_thirdpoint1(pt1,_nullpt)
    fi;
    if _isequal(pt,pt2) then
      _reductionproc(_thirdpoint2(pt)) 
    else
      _reductionproc(_thirdpoint1(pt,pt2))
    fi
  end);

  opposite:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _reductionproc=reductionproc,
  proc(pt1)
    if _isequal(pt1,_nullpt) then
      _nullpt
    else
      _reductionproc(_thirdpoint1(pt1,_nullpt))
    fi
    end);
  
  ext1:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _elliptic_add=elliptic_add,_double=double,_opposite=opposite,
       _reductionproc=reductionproc,
  proc(n::integer,pt)
    local L,pt1,pt2;
    if n=0 then _nullpt
    elif n<0 then ext1(-n,_opposite(pt))
    else
      L:=convert(n,base,2);
      pt1:=_nullpt;
      pt2:=pt;
      for i to nops(L) do
        if L[i]=1 then pt1:=_elliptic_add(pt1,pt2) fi;
	pt2:=_double(pt2)
      od;
      _reductionproc(pt1)
    fi
  end);

  ext2:=
  subs(_isequal=isequal,_thirdpoint1=thirdpoint1,
       _thirdpoint2=thirdpoint2,_nullpt=nullpt,
       _elliptic_add=elliptic_add,_double=double,
       _subt=subt,_opposite=opposite,
       _reductionproc=reductionproc,
  proc(n::integer,pt)
    local L,pt1,pt2,state; # pt1=P, pt2=Q, state=b
    if n=0 then _nullpt
    elif n<0 then ext2(-n,_opposite(pt))
    else
      L:=convert(n,base,2);
      pt1:=_nullpt; 
      pt2:=pt;
      state:=0;
      for i to nops(L) do
        if L[i]=0 then 
	  if state=0 then
	    pt2:=_double(pt2)
	  elif state=1 then
	    pt1:=_elliptic_add(pt1,pt2);
	    pt2:=_double(_double(pt2)); 
	    state:=0 
          else
	    state:=1
	  fi;
        else
	  if state=0 then 
	    state:=1
          else 
	    if state=1 then
	      pt1:=_subt(pt1,pt2);
	      pt2:=_double(pt2);
	      state:=11
	    fi;
	    pt2:=_double(pt2)
	  fi
        fi;
      od;
      _reductionproc(_elliptic_add(pt1,pt2))
    fi
  end);
      
  table([
  Equ=ec,
  Nullpt=nullpt,
  Isequal=eval(isequal),
  Add=eval(elliptic_add),
  Double=eval(double),
  Subt=eval(subt),
  Opposite=eval(opposite),
  Ext1=eval(ext1),
  Ext2=eval(ext2)
  ])
end: