How to calculate point addition using Jacobian coordinate-system over elliptic curves

Question

I'm writing a small project of elliptic curve cryptography, and the program works well when I use affine coordinate system, which means each point is represented by 2 coordinates (x',y').

Now I'm trying to replace affine coordinate system by jacobian coordinate system in which each point is represented by 3 coordinates (x,y,z), x' = x/z² and y' = y/z³.

I'd like to know how to transform affine coordinates to jacobian coordinates**. In some tutorials, people uses the formula: (x,y) = (x,y,1) which means the z-coordinate is always set to one. But I'm not sure if it's correct.

Then for points additions over elliptic curve, to calculate P(x1,y1,z1) + Q(x2,y2,z2) = R(x3,y3,z3). I've used the following formulas in my program:

u1 = x1.z2²
u2 = x2.z1²
s1 = y1.z2³
s2 = y2.z1³
h = u2 - u1
r = s2 - s1
x3 = r² - h³ - 2.u1.h²
Y3 = r. (U1.h² - x3) - s1.h³
z3 = z1.z2.h

But when I test my program, I get some negative coordinates e.g. (-2854978200,-5344897546224,578). And when I try to convert the result back to affine coordinate system with the formular (x'=x/z²,y'=y/z³), I get (-8545, -27679), actually the x coordinate is -8545.689.... The jacobian x coordinate isn't divisible by z².

What should I do if the coordinates are not integers? And if they're negative? I've tried to MOD with the field size of my curve, but the result isn't correct either.

So the point using jacobian coordinates (x,y,1) is correct, but not unique. All points satisfying (a^2.x,a^3.y,a) are equivalent. And in my program the curve is defined in a prime field, so when I calculate u1, u2, s1, s2 ... I should apply MOD p to each variable?

And for transforming the final result back to affine coordinates, e.g. The x coordinate, in fact it is not a division, it's a modular inverse? For example, my curve is defined in a finite prime field p=11, and I have a point using jacobian coordinates (15,3,2), to transform jacobian x coordinate to affine x coordinate, I have to calculate 2^2 = 4 => x = 4^-1 mod p => x = 3, and 15.3 mod p = 1, so the affine x coordinate is 1, is that right?

The objective of jacobian coordinates is to avoid the division during addition. But as Thomas Pornin said, when we calculate P1 + P2 = P3, there are some special cases to handle.

1. P1 and P2 are both infinite: P3=infinite.
2. P1 is infinite: P3=P2.
3. P2 is infinite: P3=P1.
4. P1 and P2 have the same x coordinate, but different y coordinates or both y coordinates equal 0: P3=infinite.
5. P1 and P2 have different x coordinate: Addition formula.
6. P1 and P2 have same coordinates: Doubling formula.

And here's prototypes of my C functions:

jac_addition(jacobian *, point *, jacobian *);
jac_doubling(jacobian *, jacobian *);

point is a structure representing a point defined in affine coordinate system, and jacobian for jacobian system.

The problem is when I handle those special cases, especially the 4th one, I still have convert both points back to affine coordinates, or I can't compare their coordinates, which means I still have to calculate the division.

Answer 1

The Jacobian form of the projective coordinates (as any other form) is not unique - for every value of Z (other than 0) you get other X and Y without the actual point changing.

Thus, if you have a point in affine coordinates (X', Y'), the pair (X', Y', 1) is a projective representative of this point, as well as (4·X', 8·Y', 2), (9·X', 27·Y', 3), etc. The one with 1 is the easiest to create, so usually you would use this one.

While one can define (and study) elliptic curves over any field, and many mathematicians study curves defined over the complex numbers, for cryptographic uses, the coordinates are elements of some finite field. In the simplest case, we have a prime field (i.e. integers modulo a prime number), and you'll have to do addition, subtraction, multiplication and division (and probably exponentiation) in this field.

As long as Z is not zero, you should be able to divide by Z² - this is equivalent to multiplying by the inverse of Z², and such an element exists, and can be calculated efficiently using the extended euclidean algorithm.

This is easiest if your programming language comes with some big number library which has the necessary operations predefined, like Java's BigInteger class (with its mod, modPow and modInverse methods).

The field in question (i.e. the modulus) is part of the definition of the elliptic curve, and operations over one field give totally different results than operations in another one. So make sure you are using the right field.

Answer 2

When dealing with elliptic curves, the coordinates are in a field. For cryptography, this is a finite field; in your case, the "integers modulo a prime p". All operations are meant in that field, i.e. you should do every single addition, multiplication or inversion modulo p.

When doing addition of points, there are a few special cases which you must handle specially:

• There is a special "point at infinity" which does not have coordinates x and y. It is the "zero" of the curve addition. In a generic point addition routine, you must have a way to encode a "point at infinity" and check it specially.

• When adding (x,y) to (x',y'), it may happen that x = x'. In that case, either y =y', and then it is a point doubling which has its specific formula (if you apply the generic formula, you end up dividing by zero, which will not work); or, y = -y', in which case the sum is the "point at infinity".

So the generic formula is to be applied only once you have handled the special case. In all generality, in a curve y² = x³+a·x+b, the sum of (x₁,y₁) and (x₂,y₂) is (x₃,y₃) such that x₃ = f²-x₁-x₂ and y₃ = f·(x₁-x₃)-y₁, where f = (y₂-y₁)/(x₂-x₁). This implies computing one division and two multiplications, and a few subtractions (all operations being done on integers modulo p, as explained above).

Division and inversion modulo p are relatively expensive (a modular division has typically the same cost as about 80 multiplications) so in some situations, we use projective or Jacobian coordinate systems. Jacobian coordinates are about representing a point as three values (X,Y,Z) (all of them in the finite field, i.e. integers modulo p) such that x = X/Z² and y = Y/Z³.

Each point (x,y) has many possible representations as (X,Y,Z). Conversion to Jacobian coordinates is easy by setting X = x, Y = y and Z = 1: (x,y,1) is a perfectly valid Jacobian representation of the (x,y) point. Conversion from Jacobian coordinates is computationally harder: you have to do a modular inversion, and a few multiplications (you compute U = 1/Z, then x = X·U² and y = Y·U³).

With Jacobian coordinates, the addition of two points is done in a dozen field multiplications, and no division at all. You get only a Jacobian representation of the result, so you still have to do a modular inversion or division at some point; however (and that's why you bother using Jacobian coordinates), that division can be mutualized. If you have to do a hundred or so successive point additions (as is typical in a cryptographic context, when "multiplying" a point with an integer), then you can use Jacobian representations throughout, and do a single conversion back to cartesian coordinates (x,y) at the end. So instead of doing 200 multiplications and 100 divisions, you do 1000 multiplications and 1 inversion; since an inversion costs the same as 80 multiplications, the gain is appreciable.

Try to avail yourself to this book; any good college library should have one. It explains all of that very clearly.

Answer 3

Here is a sample in python:

def to_jacobian((Xp, Yp)):
    """
    Convert point to Jacobian coordinates

    :param (Xp,Yp,Zp): First Point you want to add
    :return: Point in Jacobian coordinates
    """
    return (Xp, Yp, 1)

def from_jacobian((Xp, Yp, Zp), P):
    """
    Convert point back from Jacobian coordinates

    :param (Xp,Yp,Zp): First Point you want to add
    :param P: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    :return: Point in default coordinates
    """
    z = inv(Zp, P)
    return ((Xp * z**2) % P, (Yp * z**3) % P)

def jacobian_add((Xp, Yp, Zp), (Xq, Yq, Zq), A, P):
    """
    Add two points in elliptic curves

    :param (Xp,Yp,Zp): First Point you want to add
    :param (Xq,Yq,Zq): Second Point you want to add
    :param P: Prime number in the module of the equation Y^2 = X^3 + A*X + B (mod p)
    :param A: Coefficient of the first-order term of the equation Y^2 = X^3 + A*X + B (mod p)
    :return: Point that represents the sum of First and Second Point
    """
    if not Yp:
        return (Xq, Yq, Zq)
    if not Yq:
        return (Xp, Yp, Zp)
    U1 = (Xp * Zq ** 2) % P
    U2 = (Xq * Zp ** 2) % P
    S1 = (Yp * Zq ** 3) % P
    S2 = (Yq * Zp ** 3) % P
    if U1 == U2:
        if S1 != S2:
            return (0, 0, 1)
        return jacobian_double((Xp, Yp, Zp), A, P)
    H = U2 - U1
    R = S2 - S1
    H2 = (H * H) % P
    H3 = (H * H2) % P
    U1H2 = (U1 * H2) % P
    nx = (R ** 2 - H3 - 2 * U1H2) % P
    ny = (R * (U1H2 - nx) - S1 * H3) % P
    nz = (H * Zp * Zq) % P
    return (nx, ny, nz)

So you can sum with:

def fast_add(a, b, A, P):
    return from_jacobian(jacobian_add(to_jacobian(a), to_jacobian(b), A, P), P)