Cryptography/Prime Curve/Jacobian Coordinates

&lt;Cryptography
Introduction.
Jacobian Coordinates are used to represent elliptic curve points on prime curves formula_1 They give a speed benefit over "Affine Coordinates" when the cost for field inversions is significantly higher than field multiplications. In "Jacobian Coordinates" the triple formula_2 represents the affine point formula_3.
Point Doubling (4M + 6S or 4M + 4S).
Let (X, Y, Z) be a point (unequal to the "point at infinity") represented in "Jacobian Coordinates". Then its double (X', Y', Z') can be calculated by
 if (Y == 0)
 return POINT_AT_INFINITY
 S = 4*X*Y^2
 M = 3*X^2 + a*Z^4
 X' = M^2 - 2*S
 Y' = M*(S - X') - 8*Y^4
 Z' = 2*Y*Z
 return (X', Y', Z')
Note: if a = -3, then M can also be calculated as M = 3*(X + Z^2)*(X - Z^2), saving 2 field squarings.
Repeated Doubling ((4m-1)M + (4m+2)S).
Let P=(X, Y, Z) be a point (unequal to the "point at infinity") represented in "Jacobian Coordinates". Then, in the case a = -3, its "m-fold double" (2^m)P can be calculated as follows:
 Y = 2*Y
 W = Z^4
 while(m--)
 if (Y == 0)
 return POINT_AT_INFINITY
 A = 3*(X^2 - W)
 B = X*Y^2
 X = A^2 - 2B
 Z = Z*Y
 if (m)
 W = W*Y^4
 Y = 2*A*(B - X) - Y^4
 Y = Y/2
 return (X, Y, Z)
An alternative repeated doubling routine with costs (4m)M + (4m+2)S for any value a can be derived from the Modified Jacobian doubling routine. For small values a (say 0 or -3) the costs reduce to (4m-1)M + (4m+2)S, competing nicely with the algorithm showed above.
Point Addition (12M + 4S).
Let (X1, Y1, Z1) and (X2, Y2, Z2) be two points (both unequal to the "point at infinity") represented in "Jacobian Coordinates". Then the sum (X3, Y3, Z3) can be calculated by
 U1 = X1*Z2^2
 U2 = X2*Z1^2
 S1 = Y1*Z2^3
 S2 = Y2*Z1^3
 if (U1 == U2)
 if (S1 != S2)
 return POINT_AT_INFINITY
 else 
 return POINT_DOUBLE(X1, Y1, Z1)
 H = U2 - U1
 R = S2 - S1
 X3 = R^2 - H^3 - 2*U1*H^2
 Y3 = R*(U1*H^2 - X3) - S1*H^3
 Z3 = H*Z1*Z2
 return (X3, Y3, Z3)
J + J -&gt; J ( 12M, 2S) (secp256k1 has 12M, 4S).
 U1 = Y2 * Z1
 U2 = Y1 * Z2
 V1 = X2 * Z1
 V2 = X1 * Z2
 if (V1 == V2)
 if (U1 != U2)
 return POINT_AT_INFINITY
 else
 return POINT_DOUBLE(X1, Y1, Z1)
 U = U1 - U2
 V = V1 - V2
 W = Z1 * Z2
 A = U ^ 2 * W - V ^ 3 - 2 * V ^ 2 * V2
 X3 = V * A
 Y3 = U * (V ^ 2 * V2 - A) - V ^ 3 * U2
 Z3 = V ^ 3 * W
 return (X3, Y3, Z3)
Mixed Addition (with affine point) (8M + 3S).
Let (X1, Y1, Z1) be a point represented in "Jacobian Coordinates" and (X2, Y2) a point in "Affine Coordinates" (both unequal to the "point at infinity"). A formula to add those points can be readily derived from the regular jacobian point addition by replacing each occurrence of "Z2" by "1" (and thereby dropping four field multiplications and one field squaring).
Mixed Addition (with chudnovsky point) (11M + 3S).
Let (X1, Y1, Z1) be a point represented in "Jacobian Coordinates" and (X2, Y2, Z2, Z2^2, Z2^3) a point in "Chudnovsky Coordinates" (both unequal to the "point at infinity"). Then the sum (X3, Y3, Z3) can be readily calculated using the addition formula given above (saving one field multiplication and one field squaring).