# Elliptic curve routines
# An elliptic curve is stored as the vector of coefficients of a
#   homogeneous cubic polynomial in three variables.   There are
#   nine coefficients  a[j],  -3 <= j <= 6 where  j = 6 - w,
#   w = weight in x,y  when  x  has weight 2 and  y  has weight 3
#   with the exception of  a_5, the coefficient of y^2.  Thus,
#   F(x,y) =  a[5] y^2 + a[1] x*y + a[3] y + a[-3] y^3 + a[-2] x*y^2
#        +  a[-1] x^2*y - a[0] x^3 - a[2] x^2 - a[4] x - a[6]
#
`help/text/elliptic_curves` := TEXT(`HELP FOR: Elliptic Curves`,` `,
`SYNOPSIS: Package of routines for the arithmetic of elliptic curves`,` `,
`FUNCTIONS: addell, affell, decell, encell, fpoell, homogell, isoncurve`,
`iswform, makell, negell, numell, powell, tgtell, winitell`,` `,
`NOTE: This package also loads uses the following global names:`,
`      X,  Y,  Z,  ellmonom`):
`help/text/elliptic` := `help/text/elliptic_curves`:
`help/text/ellipticcurves` := `help/text/elliptic_curves`:
print(`Elliptic curve routines; global names used:`);
print(`Variable names: X,Y,Z,ellmonom`);
print(`Funtion names: winitell,homogell,affell,isoncurve,tgtell,iswform,`);
print(`  negell,addell,numell,makell,powell,fpoell,encell,decell`);
ellmonom:=array(-3..6,[Y^3,X*Y^2,X^2*Y,-X^3,X*Y*Z,-X^2*Z,Y*Z^2,-X*Z^2,
   Y^2*Z,-Z^3]):
# winitell(v, a, l) returns the full vector  a  when  v  is the short vector
#    <a[0], a[1], a[2], a[3], a[4], a[6]> for a Weierstrass normal form
#    optional arg.  l  is a prime indicating arithmetic mod  l
`help/text/winitell` :=
TEXT(`HELP FOR: winitell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: winitell(v, a, l)`,
`       -- v  length 6 list of coefficients for Weierstrass normal form`,
`       -- a  OUTPUT array a[j], j=-3..6, of all coefficients for cubic`,
`       -- l  optional prime indicating field of coefficients`, ` `,
`SYNOPSIS:  Initialize the full array of coefficients for the general`,
`       cubic polynomial`,
`     f(x,y) =  a[5] y^2 + a[1] x*y + a[3] y + a[-3] y^3 + a[-2] x*y^2`,
`        +  a[-1] x^2*y - a[0] x^3 - a[2] x^2 - a[4] x - a[6]`,
`       from the list [a[0],a[1],a[2],a[3],a[4],a[6]] when the`,
`       cubic is to be in Weierstrass normal form.`,` `,
`NOTES:  "winitell" does not verify non-degeneracy.`,
`        "winitell" returns the empty expression sequence.`,` `,
`EXAMPLE:`,
`> winitell([1,0,0,0,-1,0],ea);`,
`> affell(ea);`,
``,
`                                    3        2`,
`                                 - X  + X + Y`,
` `,
`SEE ALSO:  iswform`):
winitell := proc(v, a, l) local j;
if linalg[vectdim](v) <> 6 then
  print(`winitell: First arg. must have six coords`);
  RETURN()
fi;
a := array(-3..6,[0,0,0,v[1],v[2],v[3],v[4],v[5],1,v[6]]);
if nargs > 2 then if l > 0 then
for j from 0 to 6 do a[j] := a[j] mod l; od;
fi; fi;
RETURN()
end:
# homogell(a) = homogeneous cubic in X,Y,Z with coefficient vector a[j]
`help/text/homogell` :=
TEXT(`HELP FOR: homogell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: homogell(a)  -- a  an array of type "-3..6"`,` `,
`SYNOPSIS: returns the homogeneous cubic polynomial given by forming`,
`    the "inner product" of the coefficient array  a  with the monomial`,
`    array  ellmonom ,  which is initialized as`,
`array(-3..6,[Y^3,X*Y^2,X^2*Y,-X^3,X*Y*Z,-X^2*Z,Y*Z^2,-X*Z^2,Y^2*Z,-Z^3])`,
` `,`SEE ALSO: affell`):
homogell := proc(a) local F,j;
F:=0;
for j from -3 to 6 do
F := F + a[j]*ellmonom[j];
od;
F
end:
# affell(a) = affine polynomial in X,Y with coefficient vector a[j]
`help/text/affell` :=
TEXT(`HELP FOR: affell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: affell(a)  -- a  an array of type "-3..6"`,` `,
`SYNOPSIS: returns the affine cubic polynomial given by forming`,
`    the "inner product" of the coefficient array  a  with the monomial`,
`    array  ellmonom ,  which is initialized as`,
`array(-3..6,[Y^3,X*Y^2,X^2*Y,-X^3,X*Y*Z,-X^2*Z,Y*Z^2,-X*Z^2,Y^2*Z,-Z^3]),`,
`    and evaluated with  Z = 1.`,
` `,`SEE ALSO: homogell`):
affell := proc(a);
subs(Z=1, homogell(a))
end:
# isoncurve(a,p,l) = true if  p  is a point of the elliptic curve  a
`help/text/isoncurve` :=
TEXT(`HELP FOR: isoncurve    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: isoncurve(a, p, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- p  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns "true" or "false" according to whether the given`,
`           point is on the curve in the specified field, which is the`,
`           rational field if  l  is either omitted or set to zero.`,
` `):
isoncurve := proc(a,p,l) local d,s;
d := linalg[vectdim](p);
if d < 2 or d > 3 then
  print(`isoncurve: Point -- second arg -- must have 2 or 3 coords`);
  RETURN()
fi;
if d = 3 then
  s := subs({X=p[1],Y=p[2],Z=p[3]}, homogell(a));
else
  s := subs({X=p[1],Y=p[2],Z=1}, homogell(a));
fi;
if nargs > 2 then if l > 0 then s := s mod l; fi; fi;
if s = 0 then true else false fi;
end:
# tgtell(e,p) = linear equation for the tangent to e at p
#   homogeneous if  p  has 3 coords, affine if  p  has 2 coords
`help/text/tgtell` :=
TEXT(`HELP FOR: tgtell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: tgtell(a, p, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- p  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns the linear polynomial for the tangent line at`,
`           the given point for the given elliptic curve.  Uses the`,
`           affine or homogeneous polynomial according to whether`,
`           the point is given with two (affine) or three (homogeneous)`,
`           coordinates.`,
` `):
tgtell := proc(e,p,l) local F,FX,FY,FZ,d,r;
d := linalg[vectdim](p);
F := homogell(e);
if d = 3 then
  FX := subs({X=p[1],Y=p[2],Z=p[3]}, diff(F,X));
  FY := subs({X=p[1],Y=p[2],Z=p[3]}, diff(F,Y));
  FZ := subs({X=p[1],Y=p[2],Z=p[3]}, diff(F,Z));
  r := FX*X+FY*Y+FZ*Z;
else
  FX := subs({X=p[1],Y=p[2],Z=1}, diff(F,X));
  FY := subs({X=p[1],Y=p[2],Z=1}, diff(F,Y));
  FZ := subs({X=p[1],Y=p[2],Z=1}, diff(F,Z));
  r := FX*X+FY*Y+FZ;
fi;
if nargs > 2 then if l > 0 then r := r mod l; fi; fi;
RETURN(r)
end:
# iswform(e) = true if e is in (generalized) Weierstrass form
`help/text/iswform` :=
TEXT(`HELP FOR: iswform    [elliptic curve routine]`,` `,
`CALLING SEQUENCE: iswform(a, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns "true" or "false" according to whether the given`,
`           elliptic curve is in (generalized) Weierstrass normal form`,
`           relative to the given field, which is the rational field if`,
`           l  is omitted or is zero.`,` `,
`SEE ALSO:  winitell`,
` `):
iswform := proc(e,l)  local a,b,c,d;
a := e[-3]; b := e[-2]; c := e[-1]; d := e[5];
if nargs > 1 then if l > 0 then
  a := a mod l;
  b := b mod l;
  c := c mod l;
  d := d mod l;
fi; fi;
if a = 0 and b = 0 and c = 0 and d <> 0 then true
else false fi;
end:
#negell(e,p) = negative of p in the group law on e if e is Weierstrass
`help/text/negell` :=
TEXT(`HELP FOR:  negell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  negell(a, p, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- p  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns, as a list, the negative, with respect to the`,
`           group law on the elliptic curve specified by  a ,  of the`,
`           point on that curve specified by  p .  The optional prime`,
`           l  specifies the field for this computation.  The field is`,
`           the rational field if  l  is omitted or zero.`,` `):
negell := proc(e,r,l) local dr,ll,p,q,s,z;
if nargs > 2 then ll := l; else ll := 0; fi;
if not iswform(e,ll) then
  print(`negell: Curve not in Weierstrass form`);
  RETURN()
fi;
dr := linalg[vectdim](r);
if dr < 2 or dr > 3 then
  print(`negell: Wrong no. of coords. in 2nd arg.`);
  RETURN()
fi;
if not isoncurve(e,r,ll)
  then print(`negell: Point is not on curve`);
  RETURN()
fi;
if dr = 2 then p := [r[1],r[2],1]; else p := [r[1],r[2],r[3]]; fi;
if nargs = 2 then
  z := p[3];
elif l > 0 then
  z := p[3] mod l;
else
  z := p[3];
fi;
if z = 0 then RETURN([0,1,0]) fi;
# s = sum of roots of quadratic eqn for y
s := -(e[1]*p[1]*p[3]+e[3]*p[3]^2)/(e[5]*p[3]);
q := [p[1],s-p[2],p[3]];
if nargs > 2 then if l > 0 then q := q mod l; fi; fi;
if dr = 3 then
  RETURN(q)
else
  RETURN([q[1],q[2]])
fi;
end:
# addell(e,p,q) = sum of p and q in the group law on e if e is Weierstrass
`help/text/addell` :=
TEXT(`HELP FOR:  addell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  addell(a, p, q, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- p  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- q  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns, as a list, the sum, with respect to the`,
`           group law on the elliptic curve specified by  a ,  of the`,
`           points on that curve specified by  p  and  q .  The optional`,
`           prime  l  specifies the field for this computation.  The field`,
`           is the rational field if  l  is omitted or zero.`,` `):
addell := proc(e,p,q,l) local ll,r,dp,dq,dr,eq,pp,qq,tf,j,t,u,mid,s,m,ans;
if nargs = 3 then ll := 0 else ll := l; fi;
if not iswform(e,ll) then
  print(`addell: Curve not in Weierstrass form`);
  RETURN()
  fi;
dp := linalg[vectdim](p);
if dp < 2 or dp > 3 then
  print(`addell: Wrong no. of coords. in 2nd arg.`);
  RETURN()
fi;
dq := linalg[vectdim](q);
if dq < 2 or dq > 3 then
  print(`addell: Wrong no. of coords. in 3rd arg.`);
  RETURN()
fi;
if not isoncurve(e,p,ll) then
  print(`addell: 2nd arg. is not on curve`);
  RETURN()
fi;
if not isoncurve(e,q,ll) then
  print(`addell: 3rd arg. is not on curve`);
  RETURN()
fi;
if dp = 2 then pp := array(1..3,[p[1],p[2],1]);
  else pp := array(1..3,[p[1],p[2],p[3]]);
fi;
if ll > 0 then pp := map(modp,pp,ll); fi;
if dq = 2 then qq := array(1..3,[q[1],q[2],1]);
  else qq := array(1..3,[q[1],q[2],q[3]]);
fi;
if ll > 0 then qq := map(modp,qq,ll); fi;
#  Finish quickly if either summand is the origin
dr := max(dp,dq);
if pp[3] = 0 then
  if dq = 3 then RETURN([qq[1],qq[2],qq[3]]) else RETURN([qq[1],qq[2]]) fi;
fi;
if qq[3] = 0 then
  if dp = 3 then RETURN([pp[1],pp[2],pp[3]]) else RETURN([pp[1],pp[2]]) fi;
fi;
#  Finish quickly if the two points have the same x-coordinate
if pp[1] = qq[1] and pp[2] <> qq[2] then RETURN([0,1,0]) fi;
#  Normalize third coordinate
if ll > 0 then
for j to 2 do pp[j] := pp[j]/pp[3] mod ll; od;
for j to 2 do qq[j] := qq[j]/qq[3] mod ll; od;
else
for j to 2 do pp[j] := pp[j]/pp[3]; od;
for j to 2 do qq[j] := qq[j]/qq[3]; od;
fi;
tf := 0;
for j to 2 do
  if pp[j] <> qq[j] then tf := 1; break; fi;
od;
if tf = 0 then
 # Use tangent at p
 eq := tgtell(e,[pp[1],pp[2]],ll);
 m := array(1..2,[diff(eq,X),diff(eq,Y)]);
 if ll > 0 then m := map(modp,m,ll); fi;
 if m[2] = 0 then RETURN([0,1,0]) fi;
 # Case where tangent is not vertical
 u:=array(1..2);
 u[1] := pp[1] - m[2]*t;
 u[2] := pp[2] + m[1]*t;
 eq := expand(subs({X=u[1],Y=u[2]}, affell(e)));
 if ll > 0 then eq := eq mod ll; fi;
 mid := quo(eq,t^2,t);
 if ll > 0 then mid := mid mod ll; fi;
else
 # Points  p  and  q  distinct
 u :=array(1..2);
 for j to 2 do u[j] := t*pp[j] + (1-t)*qq[j] od;
 eq := expand(subs({X=u[1],Y=u[2]}, affell(e)));
 if ll > 0 then eq := eq mod ll; fi;
 mid := quo(eq,t^2-t,t);
 if ll > 0 then mid := mid mod ll; fi;
fi;
s := solve(mid,t);
if ll > 0 then s := s mod ll; fi;
r := array(1..3);
r[3] := 1;
for j to 2 do r[j] := subs(t=s,u[j]); od;
if ll > 0 then for j to 2 do r[j] := r[j] mod ll od; fi;
## print(`addell: 3rd point on chord`, [r[1],r[2],r[3]]);
## if not isoncurve(e,[r[1],r[2],r[3]],ll) then
##   print(`addell: Warning -- 3rd point not on curve`);
## fi;
ans := negell(e,[r[1],r[2],r[3]],ll);
## print(`addell: negell returned`, ans);
if dr = 2 then RETURN([ans[1],ans[2]]) else RETURN(ans) fi;
end:
#numell(e,p) = number of point on  e  [in Weierstrass form for now]  over the
#              field of  p  elements,  p  prime for now
`help/text/numell` :=
TEXT(`HELP FOR:  numell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  numell(a, l)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:  Returns the number of points in the field of  l  elements`,
`           on the elliptic curve specified by  a .  This total includes`,
`           all points in the projective plane on the curve.`,` `):
numell := proc(e,p) local A,B,C,d,j,t;
if not iswform(e) then
  print(`Curve is not in Weierstrass form`);
  RETURN(-1)
fi;
if modp(p, 2) = 0 then
  print(`Looking for an odd prime`);
  RETURN(-2)
fi;
A := e[5];
if A mod p = 0 then
  print(`Curve is singular in given field`);
  RETURN(-3)
fi;
t:=1;
for j from 0 to p-1 do
B := e[1]*j + e[3]; 
C := -e[0]*j^3-e[2]*j^2-e[4]*j-e[6];
d := B^2 - 4*A*C;
if d mod p = 0 then
  t := t + 1;
elif numtheory[L](d, p) = 1 then
  t := t + 2;
fi;
od;
RETURN(t)
end:
# makell(a,x,y,e,p) = an elliptic curve  Y^2 = X^3 + aX + b ,  where  a  is
#                     the first arg, containing the point  (x, y)
#                     With optional arg p, a prime, checks non-singularity
#                     at the prime p.  e  is the computed curve.
`help/text/makell` :=
TEXT(`HELP FOR:  makell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  makell(a, x, y, e, l)`,
`       --  a, x, y   numbers`,
`       --  e  OUTPUT array of type -3..6 for the coefficient vector`,
`              of a cubic polynomial`,
`       --  l  optional prime indicating the field`,` `,
`SYNOPSIS:  Computes as  e  the length 10 array of cubic polynomial`,
`           coefficients for the elliptic curve`,
`                           Y^2 = X^3 + a*X + b ,`,
`           where  a  is the first argument containing the point`,
`           (X,Y) = (x,y).  This gives a method of constructing an`,
`           elliptic curve in Weierstrass form, as above, having`,
`           a given point (x,y) on it.`,` `):
makell := proc(a,x,y,e,p) local b,d,dp;
b:=y^2-x^3-a*x;
d:= 4*a^3-27*b^2;
if nargs = 5 then
 dp := d mod p;
 if dp = 0 then
   print(`Curve is singular at given prime`);
   RETURN()
 fi;
fi;
if d = 0 then
  print(`Curve is singular`);
  RETURN()
fi;
e := array(-3..6,[0,0,0,1,0,0,0,a,1,b]);
print(affell(e));
RETURN()
end:
# powell(e,n,p) = n-th multiple of point  p  on elliptic curve  e
`help/text/powell` :=
TEXT(`HELP FOR:  powell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  powell(a, m, p, l)`,` `,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- m  an integer`,
`       -- p  list or vector of coordinates of point in plane`,
`             (affine or homogeneous coordinates)`,
`       -- l  optional prime indicating field`,` `,
`SYNOPSIS:   Returns, as a list, the  m-th  multiple, in the group`,
`            law on the elliptic curve specified by  a ,  of the point`,
`            p .  The optional prime  l  indicates the field.  The field`,
`            is taken to be the rational field if  l  is omitted or`,
`            is zero.`,` `):
#  This code is inefficient.  See "powell" below.
powell1 := proc(e,m,p,l) local a,dp,j,ll,n,u,v;
if nargs = 3 then ll := 0 else ll := l fi;
u :=array(1..3);
dp := linalg[vectdim](p);
if dp < 2 or dp > 3 then
  print(`Wrong no. of coords. in 3rd arg.`);
  RETURN()
fi;
n := m;
if n = 0 then RETURN([0,1,0]) fi;
if n < 0 then
  a := negell(e,p,ll);
  n := -n;
else
  a := p;
fi;
u := a;
for j to n-1 do
 v := addell(e,a,u,ll);
 u := v;
od;
RETURN(u)
end:
#  The following code is relatively efficient.  It requires time that
#  is proportional to the number of digits in the 2-adic expansion of
#  the "power" m.
powell := proc(e,m,p,l) local dp,j,k,ll,n,q,s,u;
if nargs = 3 then ll := 0 else ll := l fi;
dp := linalg[vectdim](p);
if dp < 2 or dp > 3 then
  print(`Wrong no. of coords. in 3rd arg.`);
  RETURN()
fi;
if m = 0 then
  RETURN([0,1,0])
elif m > 0 then
  n := m;
  q := convert(p, list);
else
  n := -m;
  q := negell(e,p,ll);
fi;
s := convert(n,base,2);
k := nops(s);
u := convert(q,list);
for j from 2 to k do
if s[k+1-j] = 1 then
  u := addell(e,q,addell(e,u,u,ll),ll);
else
  u := addell(e,u,u,ll);
fi;
od;
RETURN(u)
end:
# fpoell(e,a,b) prints a list of all primes  p ,  a <= p <= b, for which
#       the elliptic curve  e  has a prime number of points mod p
`help/text/fpoell` :=
TEXT(`HELP FOR:  fpoell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  fpoell(a, m, n)    [elliptic curve routine]`,
`       -- a    array of type -3..6 of cubic polynomial coefficients`,
`       -- m,n  postive integers,  m < n`,` `,
`SYNOPSIS:  Prints a list of all primes  l  between the integers  m`,
`           and  n  for which the elliptic curve specified by the`,
`           vector  a  has a prime number of points mod l.`,` `):
fpoell := proc(e,a,b) local p,t,u;
if not type(a,integer) then print(`2nd arg must be integer`); RETURN() fi;
if not type(b,integer) then print(`3rd arg must be integer`); RETURN() fi;
if not a > 0 then print(`1st arg must be positive`); RETURN() fi;
if not b > a then print(`2nd arg must be greater than 1st arg`); RETURN() fi;
p := nextprime(a);
u := 0;
while p < b do
t := numell(e,p);
if type(t, integer) and t > 0 then
  if isprime(t) then u := u+1; print(p, t); fi;
else
  print(`Failed to find number of points for prime`,p);
fi;
p := nextprime(p);
od;
RETURN(u)
end:
# encell(v,e,p,k) = point on elliptic curve  e  over the integers
#     mod  p  representing the number  v  with failure rate
#     1/2^k
`help/text/encell` :=
TEXT(`HELP FOR:  encell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  encell(t, a, l, k)`,
`       -- t  an integer (corresponding to a text unit)`,
`       -- a  array of type -3..6 of cubic polynomial coefficients`,
`       -- l  prime denoting a finite field`,
`       -- k  integer: large enough so that 1/2^k is an acceptable`,
`             "probability" of failure`,` `,
`SYNOPSIS:  Returns, as a list, the coordinates of a point on the`,
`           elliptic curve specified by  a  that matches the number`,
`           t  corresponding to a unit of text.  The curve must be`,
`           in Weierstrass form.  This routine computes a point (x, y)`,
`           on the curve where  x = kt + j, for the smallest positive`,
`           integer  j  such that some point  (x, y)  on the curve`,
`           exists in the field of  l  elements.  Each value of  j`,
`           from  1  to  k  is tried until a suitable  y  is found.`,` `,
`           Note that  t  may be recovered from the point  (x, y)`,
`           by the formula  t = [(x-1)/k] .`,
`SEE ALSO:  decell`,` `):
encell := proc(v,e,p,k) local a,b,c,d,f,g,j,r,u,x,y;
if not iswform(e,p) then
  print(`encell: curve is not in Weierstrass form`);
  RETURN()
fi;
if (k+1)*v > p-1 then
  print(`encell: field too small for number`,v);
  RETURN(0)
fi;
f := affell(e);
u := 0;
for j from 1 to k do
x := k*v+j;
g := subs(X=x,f) mod p;
a := coeff(g,Y,2);
b := coeff(g,Y,1);
c := coeff(g,Y,0);
d := b^2 - 4*a*c mod p;
if d = 0 then
  y := -b/(2*a) mod p;
  u := 1;
  break;
elif numtheory[L](d,p) = 1 then
  r := numtheory[msqrt](d,p);
  y := (r-b)/(2*a) mod p;
  u := 1;
  break;
fi;
od;
if u = 1 then RETURN([x,y]) else RETURN(0) fi;
end:
# decell(p,l,k) = unit of text recovered from a point in the plane
#      over the field of  l  elements,  l  prime, following probabilistic
#      matching of text with the points of an elliptic curve having
#      failure rate  1/2^k
`help/text/decell` :=
TEXT(`HELP FOR:  decell    [elliptic curve routine]`,` `,
`CALLING SEQUENCE:  decell(p,l,k)`,
`       -- p  point in the plane over the field of  l  elements, or`,
`             the first affine coordinate of such a point`,
`       -- l  prime corresponding to a given field`,
`       -- k  negative base 2 logarithm of the failure rate associated`,
`             with the process of assigning a text unit to a point`,
`             on an elliptic curve`,` `,
`SYNOPSIS:  Returns a non-negative integer representing a text unit`,
`           that was previously associated, using encell, with a point`,
`           on an elliptic curve.`,` `,
`SEE ALSO:  encell`,` `):
decell := proc(p,l,k) local d,t,x;
if type(p,integer) then
  x := modp(p,l);
else
  d := linalg[vectdim](p);
  if d < 3 then
    x := modp(p[1],l);
  elif d = 3 then
    if not modp(p[3],l) = 0 then x := modp(p[1]/p[3],l); fi;
  else 
    print(`decell: 1st arg wrong no. of coords.`);
    RETURN()
  fi;
fi;
t := trunc((x-1)/k);
RETURN(t)
end:
print(`Pre-defining global names 'e0', 'e1', 'ec'`);
winitell([1,0,0,1,-1,0],e0);
winitell([1,0,1,1,0,0],e1);
winitell([1,0,0,0,0,1],ec);