==============================================================================
This is a (preliminary version of) a short description of an efficient
algorithm for computing the group law in the Jacobian of a genus-2
curve.

The first section gives definitions and notation and describes some
well-known properties.

The second gives pseudo-code for the algorithm.  The details of some
complex formulae are omitted.

The third discusses in detail the complex formulae for the frequent
case when divisors are "in general position" and shows why the
algorithm requires 2 inversion and 27 multiplications in such a case.

Rob Harley, July 2000.
==============================================================================
1. Definitions and notation.

Let F_q be the finite field with q = p^d elements.

An elliptic or hyperelliptic curve over F_q is given by an equation:
  y^2 = f(x)
where f is a monic polynomial over F_q of degree 2g+1 and with distinct
roots.  The integer g is the GENUS.  For elliptic curves, g = 1.  For
hyperelliptic curves, g >= 2.

We will use letters P, Q, R for arbitrary (finite) points on the
curve.  When P = (x_P, y_P) is such a point, then the point -P = (x_P,
-y_P) is too and is called the OPPOSITE of P.  Note that x_P and y_P
may be in any extension field of F_q.  We say that a point is DEFINED
OVER a field F when its coordinates are in F.


A DIVISOR is a formal sum of at most g points with the property that
no two of them are opposite.  We write extra parentheses around such a
formal sum to turn it into a divisor.  The number of points is called
the WEIGHT of the divisor.

In practice we employ MUMFORD's representation of divisors.  Each
divisor (Sum_i P_i) is represented as a pair of polynomials (u, v)
where u has its roots at the abscissae x_{P_i} and v interpolates the
points P_i.  Here u may have repeated roots and v is tangent to the
curve, both according to multiplicity.

More precisely the divisors are (u, v) where:
  u is monic of degree at most g,
  v is reduced modulo u,
  u | f-v^2.

For instance the unique divisor of weight 0 is O = (1, 0).  The
divisor of weight 1, with a single point P = (x_P, y_P) on the curve,
is (x-x_P, y_P).

The set of divisors can be given a group law.  The resulting group is
called the JACOBIAN.  The divisor O is the neutral element of the
group.  We will use letters D, E for arbitrary divisors.  The opposite
of a divisor D = (u, v) is -D = (u, -v).  A divisor D = (u, v) is
2-torsion i.e., D+D = O, if and only if D = -D in other words iff v = 0.

We say that a divisor is defined over a field F when the coefficients
of u and v are in F.  Beware that the points composing such a divisor
are not necessarily defined over F although their coordinates are
algebraic over F.

The sum of two divisors defined over F is again defined over F.  Thus
the set of divisors defined over F forms a subgroup of the Jacobian.
It is possible to compute in this subgroup by using Mumford's
representation and doing arithmetic over F without going into
extension fields.  Two algorithms for doing so in the general case
were given by CANTOR.

We will concentrate on genus 2 and give a very efficient algorithm for
adding divisors in Jacobians of genus-2 curves.

The divisors of weight 0 or 1 were described above.  The remaining
ones contain two points P = (x_P, y_P) and Q = (x_Q, y_Q) on the curve
and are of the form D = (u, v) with u = (x-x_P)(x-x_Q).  When D
is defined over F then the coordinates x_P and x_Q are either in F or
are conjugate elements in a quadratic extension, and likewise for y_P
and y_Q.

In the following we describe the genus-2 group law using pseudo-code.

Note that the symbol + between points denotes the formal sum whereas
between divisors it denotes the group law.  For instance (P+Q) is a
divisor of weight 2, with necessarily P != -Q, whereas (P)+(Q) is the
sum of two divisors of weight 1, with possibly P = -Q.  In fact
(P)+(Q) = O when P = -Q and (P)+(Q) = (P+Q) otherwise.

The notation [n]D denotes the sum D+D+...+D repeated n times.

We write D.u and D.v for the u and v polynomials of divisor D in
Mumford's representation.  Their coefficients are given by u(x) =
u_0*x^2+u_1*x+u_2 and v(x) = v_0*x+v_1.

==============================================================================
2. The algorithm in pseudo-code.

Here is a description of the method for adding divisors.  It is given
as a list of functions, essentially a large case-by-case analysis
which causes no particular difficulty.

Some complex formulae are omitted for simplicity.  However the detail
of the most frequently occurring one is shown in section 3.

An old version of the doubling algorithm is available in doubling.c.
Complete C code will be available soon (it is too messy for public
consumption at the moment!)

------------------------------------------------------------------------------
Function: jacDoubleOneGivingTwo()
Input:    P with y_P non-zero.
Output:   [2](P)

/* Use simple formula. */
Set s = f'(x_P)/(2*y_P) and t = y_P-s*x_P.
Set u(x) = (x-x_P)^2 and v(x) = s*x+t.
Result is (P+P) = (u, v).

------------------------------------------------------------------------------
Function: jacDoubleOne()
Input:    P
Output:   [2](P)

IF y_P = 0
THEN
  /* Divisor (P) has order 2. */
  Result is O.
ELSE
  Call jacDoubleOneGivingTwo().

------------------------------------------------------------------------------
Function: jacDoubleTwo()
Input:    D of weight 2, say D = (u, v)
Output:   [2]D

IF v = 0
THEN
  /* Divisor D has order 2. */
  Result is O.
ELSE
  BEGIN
    /* Get r = resultant(u, v). */
    Set r = u_2*v_0^2-u_1*v_0*v_1+v_1^2.
    IF r = 0
    THEN
      /* One point has zero ordinate, get the other. */
      Set x_P = v1/v0-u1 and y_P = v(x_P).
      Call jacDoubleOneGivingTwo() to get result (P+P).
    ELSE
      Use complex formula for result of weight 1 or 2.
  END

------------------------------------------------------------------------------
Function: jacDouble()
Input:    D
Output:   [2]D

IF D has weight 2
THEN
  Call jacDoubleTwo().
ELSE IF D has weight 1
THEN
  Call jacDoubleOne().
ELSE /* D = O. */
  Result is O.

------------------------------------------------------------------------------
Function: jacSumOneOne()
Inputs:   P, Q
Output:   (P)+(Q)

IF P = -Q
THEN
  Result is O.
ELSE IF P = Q
THEN
  Call jacDoubleOneGivingTwo() to get result (P+P).
ELSE
  /* Use simple formula. */
  Set s = (y_P-y_Q)/(x_P-x_Q) and t = y_P-s*x_P.
  Set u(x) = (x-x_P)(x-x_Q) and v(x) = s*x+t.
  Result is (P+Q) = (u, v).

------------------------------------------------------------------------------
Function: jacSumOneTwoCoprime()
Inputs:   P, E with E of weight 2 and neither P nor -P occurs in E
Output:   (P)+E

Use complex formula. /* The result has weight 2. */

------------------------------------------------------------------------------
Function: jacSumOneTwoCommon()
Inputs:   P, E with E of weight 2 and either P or -P occurs in E.
Output:   (P)+E

Let E be (u, v).
Set y = v(x_P).

IF y = -y_P
THEN
  /* Point -P occurs in E, so return other point. */
  Set x_Q = -u_1-x_P and y_Q = v(x_Q).
  Result is (Q).
ELSE IF u_1 = -2*x_P
THEN
  /* Here E = (P+P). */
  Use complex formula to get result (P+P)+(P). /* Note it has weight 2. */
ELSE
  /* Here E is of the form (P+Q) with P != -P and P != +/-Q. */
  Extract Q by x_Q = -u_1-x_P and y_Q = v(x_Q).
  Call jacDoubleOneGivingTwo() to get (P+P).
  Call jacSumOneTwoCoprime() to get result (P+P)+(Q).

------------------------------------------------------------------------------
Function: jacSumOneTwo()
Inputs:   P, E with E of weight 2.
Output:   (P)+E

IF E.u(x_P) = 0
THEN
  Call jacSumOneTwoCommon().
ELSE
  Call jacSumOneTwoCoprime().

------------------------------------------------------------------------------
Function: jacSumOneN()
Inputs:   P, E
Output:   (P)+E

IF E has weight 2
THEN
  Call jacSumOneTwo()
ELSE IF E has weight 1
THEN
  Call jacSumOneOne()
ELSE
  Result is (P).

------------------------------------------------------------------------------
Function: jacSumTwoTwo()
Inputs:   D, E with D and E of weight 2
Output:   D+E

IF D.u = E.u
THEN
  BEGIN
    IF D.v = -E.v
    THEN
      /* Opposite divisors. */
      Result is O.
    ELSE IF D.v = E.v
    THEN
      /* Equal divisors. */
      Call jacDouble().
    ELSE
      /* Have divisors of the form D = (P+Q) and E = (P-Q). */
      /* Here Q != -Q so D.v0-E.v0 is non-zero. */
      Extract P by x_P = (E.v1-D.v1)/(D.v0-E.v0) and yP = D.v(x_P).
      Call jacDoubleOneGivingTwo() to get result (P+P).
  END
ELSE
  BEGIN
    /* Get r = resultant(D.u, E.u). */
    Set t = D.u2 - D.u1*E.u1 + E.u1^2 - E.u2.
    Set r = D.u2*(t - E.u2) + E.u2*(D.u1^2 - D.u1*E.u1 + E.u2).
    IF r != 0
    THEN
      Use complex formula (frequent case).
    ELSE
      BEGIN
        /* Here D.u and E.u have a linear factor in common. */
        /* I.e., we have D = (P+Q) and E = (+/-P+R) with Q != +/-R. */
        Set x_P = (E.u2-D.u2)/(D.u1-E.u1).
        Set yy = D.v(x_P)+E.v(x_P).
        Set x_Q = -D.u1-x_P and y_Q = D.v(x_Q).
        Set x_R = -E.u1-x_P and y_R = E.v(x_R).
        IF yy = 0
        THEN
          /* Have D = (P+Q) and E = (-P+R). */
          Call jacSumOneOne() to get result (Q+R).
        ELSE
          /* Have D = (P+Q) and E = (P+R). */
          Set y_P = yy/2.
          Call jacDoubleOneGivingTwo() to get (P+P).
          Call jacSumOneTwo() to get (P+P)+(Q). /* Note Q != -P. */
          Call jacSumOneTwo() to get result (P+P)+(Q)+(R).
      END
  END
------------------------------------------------------------------------------
Function: jacSumTwoN()
Inputs:   D, E with D of weight 2
Output:   D+E

IF E has weight 2
THEN
  Call jacSumTwoTwo()
ELSE IF E has weight 1
THEN
  Call jacSumOneTwo()
ELSE
  Result is D.

------------------------------------------------------------------------------
Function: jacSum()
Inputs:   D, E
Output:   D+E

IF E has weight 2
THEN
  Call jacSumTwoN()
ELSE IF E has weight 1
THEN
  Call jacSumOneN()
ELSE
  Result is D.

==============================================================================
3. Detail of the frequent case.

Initially one calls jacSum() to compute D+E.
In the large majority of cases, the algorithm runs through as follows.

First it finds that D has weight 2 and calls jacSumTwoN().
Then it finds that E has weight 2 and calls jacSumTwoTwo().
Then it finds that D.u != E.u and computes the resultant r.
Then it finds that r is non-zero and applies the general formula to finish.

In outline, the general formula is:

  Apply Chinese Remainder Theorem to interpolate:
    Solve nv = D.v mod D.u and nv = E.v mod E.u

  Do exact division to get other roots:
    nu = (nv^2-f)/D.u/E.u

  Finally make nu monic, reduce nv modulo nu and negate it.


In practice the outlined form is optimised into:

  s = (E.v-D.v)/D.u mod E.u
  k = (f-D.v^2)/D.u                                     (exact division)
  su = s*D.u
  nu = (s*(su+2*D.v)-k)/E.u                             (exact division)
  make nu monic
  nv = (su+D.v) mod nu
  nv = -nv

The degrees of some of the quantities are given here along with the
coefficients that are actually required.  Note that k is monic.

                        degree  coeffs
                        bound   needed
        s               1       c0,c1
        k               3          c1
        su              3       c0,c1,c2,c3
        su+2*D.v        3       c0,c1,c2
        s*(that)        4       c0,c1,c2
        (that)-k        4       c0,c1,c2
        (that)/E.u      2       c0,c1,c2
        su+D.v          3       c0,c1,c2,c3
        (that) mod nu   <= 2


The final algorithm is as follows, with the costs as shown (P =
product, S = squaring, I = inversion, additions and such are ignored).

                                                          Cost
  Common subexpressions:
    de = D.u_1*E.u_1                                      P
    t = D.u_2 - de + E.u_1^2 - E.u_2                      S

  Resultant:
    r = D.u_2*(t-E.u_2) + E.u_2*(D.u_1^2-de+E.u_2)        S + 2 P

  Inverse of D.u modulo E.u:
     idu = ((E.u_1 - D.u_1)*x + t)/r                      I + 2 P

  Use modular inverse:
    s = (E.v-D.v)*idu mod E.u                             5 P

  Exact division:
    k = (f-D.v^2)/D.u                                     Free!

  Common subexpression:
    su = s*D.u                                            3 P

  Exact division:
    nu = (s*(su+2*D.v)-k)/E.u                             S + 6 P

  Make nu monic
    ...                                                   I + 2 P

  Reduce nv:
    nv = (su+D.v) mod nu                                  3 P

  Negate nv.
    nv = -nv                                              Free

                                                          --------

                                                  Total:  2 I + 3 S + 24 M
==============================================================================