Binary Curve Affine Coordinates

Point Doubling (1D + 1M)

Let (x,y) be a point (unequal to the 'point at infinity') on the elliptic (binary) curve given by the equation y^2 + xy = x^3 + ax^2 + b. Then the point (x',y') := 2*(x,y) can be computed by
 if (x == 0)
   return POINT_AT_INFINITY
 else
   l = x + (y / x)
   x' = l^2 + l + a = x^2 + b / x^2
   y' = x^2 + l*x' + x'
   return (x', y')
Note that there are two different ways to compute x'.

Point Halving

Let G be an odd order subgroup of the elliptic curve and P=(x,y) be a point (unequal to the 'point at infinity') in G. Then there exists a unique point Q = (x',y') in G with 2Q = P. If Trace(a) is 1, then Q can be computed as follows:
 find a solution l for the quadratic equation l^2 + l = x + a
 t = y + x*l
 if Trace(t) == 1
   l = l + 1
 else
   t = t + x
 x' = sqrt(t)
 y' = l*x' + t
 return (x', y')
A routine for m-fold repeated halving (of points unequal to the 'point at infinity') is given below.
 if (m == 0)
   return (x, y)
 first = 1
 x' = x
 while(m--)
   find a solution l for the quadratic equation l^2 + l = x' + a
   if (first)
     t = y + x*l
     first = 0
   else
     t = x'*(x' + lq + l)
   if Trace(t) == 1
     lq = l + 1
   else
     lq = l 
     t = t + x'
   x' = sqrt(t)
 y' = lq*x' + t
 return (x', y')

Point addition (1D + 1M)

Let (x1,y1) and (x2,y2) be two points (both unequal to the 'point at infinity'). Then the point (x3,y3) := (x1,y1) + (x2,y2) can be computed by
 if (x1 == x2)
   if (y1 != y2)
     return POINT_AT_INFINITY
   else
     return POINT_DOUBLE(x1, y1)
 l = (y2 + y1) / (x2 + x1)
 x3 = l^2 + l + x1 + x2 + a
 y3 = l(x1 + x3) + x3 + y1 = l(x2 + x3) + x3 + y2
 return (x3, y3)

Point Decompression

The following algorithm calculates for a given x a value y, such that (x,y) is a point on the elliptic curve.
 if x == 0
   return y = sqrt(b)
 else
   t = x + a + b/(x^2)
   find a solution l for the quadratic equation l^2 + l = t
   if no such solution exists 
     return POINT_NOT_EXPANDABLE
   else
     return y = l*x  (the result (l + 1)*x would be correct, too)