ec: p256

[mmcloughlin]

Implement the P-256 curve as the first end-to-end curve implementation. Operations:

• point on curve (or fallback to Params() implementation)
• add
• double 991484c
• scalar mult gen/ec: scalar multiplication #67
• scalar base mult gen/curve: optimized fixed-base scalar multiplication #88
• inverse modulo N (private crypto/ecdsa.invertable interface) gen/curve: scalar inversion #85
• combined mult (private crypto/ecdsa.combinedMult interface)

Sub-tasks:

• generate functions from EFD formulae gen/ec: code generation for efd formulae #57
• field element conversion to/from big.Int gen/fp: math/big conversion #61
• point conversion to/from big.Int
• conversion between affine and jacobian points
• handle montgomery encoding/decoding
• implementation doc comments
• point cmov gen/ec: conditional point operations #68
• point neg
• point cneg gen/ec: conditional point operations #68
• test for exceptional inputs: neutral point, addition used for doubling, ...
• review use of math/big testing: audit use of math/big #79
• field operations mul/sqr/inv modulo N gen/fp: fields modulo group order N #77
• decision on montgomery encode/decode gen: handling of montgomery encode/decode #86
• generate curve operations tests
• extend curve tests gen/curve: shortw curve tests #89
• test point operations

Performance:

• benchmark curve operations gen/curve: shortw curve benchmarks #90
• benchmark point operations
• generate field tests gen/fp: test against math/big #50
• benchmark field operations fp: field operation benchmarks #44
• ensure temporaries do not escape gen/fp: mark assembly functions go:noescape #91
• differences from stdlib gen/curve: p256 performance differences #92

[mmcloughlin]

The elliptic.Curve interface to implement is

type Curve interface {
    // Params returns the parameters for the curve.
    Params() *CurveParams
    // IsOnCurve reports whether the given (x,y) lies on the curve.
    IsOnCurve(x, y *big.Int) bool
    // Add returns the sum of (x1,y1) and (x2,y2)
    Add(x1, y1, x2, y2 *big.Int) (x, y *big.Int)
    // Double returns 2*(x,y)
    Double(x1, y1 *big.Int) (x, y *big.Int)
    // ScalarMult returns k*(Bx,By) where k is a number in big-endian form.
    ScalarMult(x1, y1 *big.Int, k []byte) (x, y *big.Int)
    // ScalarBaseMult returns k*G, where G is the base point of the group
    // and k is an integer in big-endian form.
    ScalarBaseMult(k []byte) (x, y *big.Int)
}

Coordinates

Use Jacobian coordinates under the hood.

x = X/Z^2
y = Y/Z^3

(Note converting affine to Jacobian we just take X = x, Y = y, Z = 1).

Add

The existing implementation uses add-2007-bl.

https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-2007-bl

id           g1p/shortw/jacobian-3/addition/add-2007-bl
tag          add-2007-bl
class        g1p
shape        shortw
repr         jacobian-3
operation    addition
cost         11M + 5S + 9add + 4*2
source       2007 Bernstein--Lange; note that the improvement from 12M+4S to 11M+5S was already mentioned in 2001 Bernstein http://cr.yp.to/talks.html#2001.10.29
compute      Z1Z1 = Z1^2
             Z2Z2 = Z2^2
             U1 = X1 Z2Z2
             U2 = X2 Z1Z1
             S1 = Y1 Z2 Z2Z2
             S2 = Y2 Z1 Z1Z1
             H = U2-U1
             I = (2 H)^2
             J = H I
             r = 2 (S2-S1)
             V = U1 I
             X3 = r^2-J-2 V
             Y3 = r (V-X3)-2 S1 J
             Z3 = ((Z1+Z2)^2-Z1Z1-Z2Z2) H
op3          Z1Z1 = Z1^2
             Z2Z2 = Z2^2
             U1 = X1*Z2Z2
             U2 = X2*Z1Z1
             t0 = Z2*Z2Z2
             S1 = Y1*t0
             t1 = Z1*Z1Z1
             S2 = Y2*t1
             H = U2-U1
             t2 = 2*H
             I = t2^2
             J = H*I
             t3 = S2-S1
             r = 2*t3
             V = U1*I
             t4 = r^2
             t5 = 2*V
             t6 = t4-J
             X3 = t6-t5
             t7 = V-X3
             t8 = S1*J
             t9 = 2*t8
             t10 = r*t7
             Y3 = t10-t9
             t11 = Z1+Z2
             t12 = t11^2
             t13 = t12-Z1Z1
             t14 = t13-Z2Z2
             Z3 = t14*H
inputs       X1 X2 Y1 Y2 Z1 Z2

The cloudflare p384 implementation uses add-1998-cmo

https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo

id           g1p/shortw/jacobian-3/addition/add-1998-cmo
tag          add-1998-cmo
class        g1p
shape        shortw
repr         jacobian-3
operation    addition
url          https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-add-1998-cmo
cost         10M + 5S + 4^3 + 6add + 1*2
source       1998 Cohen--Miyaji--Ono "Efficient elliptic curve exponentiation using mixed coordinates", formula (5)
compute      U1 = X1 Z2^2
             U2 = X2 Z1^2
             S1 = Y1 Z2^3
             S2 = Y2 Z1^3
             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 = Z1 Z2 H
op3          t0 = Z2^2
             U1 = X1*t0
             t1 = Z1^2
             U2 = X2*t1
             t2 = Z2^3
             S1 = Y1*t2
             t3 = Z1^3
             S2 = Y2*t3
             H = U2-U1
             r = S2-S1
             t4 = r^2
             t5 = H^3
             t6 = H^2
             t7 = U1*t6
             t8 = 2*t7
             t9 = t4-t5
             X3 = t9-t8
             t10 = H^2
             t11 = U1*t10
             t12 = t11-X3
             t13 = H^3
             t14 = S1*t13
             t15 = r*t12
             Y3 = t15-t14
             t16 = Z2*H
             Z3 = Z1*t16
inputs       X1 X2 Y1 Y2 Z1 Z2

Mixed Add

A mixed addition adds an affine point to a Jacobian point. In the EFD this is specified as an assumption that Z2 = 1, that is the second point's Z-coordinate is 1.

I cannot easily tell which formula the crypto/elliptic P-256 uses. I think cloudflare p384 uses madd-2007-bl

https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-madd-2007-bl

id           g1p/shortw/jacobian-3/addition/madd-2007-bl
tag          madd-2007-bl
class        g1p
shape        shortw
repr         jacobian-3
operation    addition
url          https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#addition-madd-2007-bl
cost         7M + 4S + 9add + 3*2 + 1*4
source       2007 Bernstein--Lange
assume       Z2=1
compute      Z1Z1 = Z1^2
             U2 = X2 Z1Z1
             S2 = Y2 Z1 Z1Z1
             H = U2-X1
             HH = H^2
             I = 4 HH
             J = H I
             r = 2 (S2-Y1)
             V = X1 I
             X3 = r^2-J-2 V
             Y3 = r (V-X3)-2 Y1 J
             Z3 = (Z1+H)^2-Z1Z1-HH
op3          Z1Z1 = Z1^2
             U2 = X2*Z1Z1
             t0 = Z1*Z1Z1
             S2 = Y2*t0
             H = U2-X1
             HH = H^2
             I = 4*HH
             J = H*I
             t1 = S2-Y1
             r = 2*t1
             V = X1*I
             t2 = r^2
             t3 = 2*V
             t4 = t2-J
             X3 = t4-t3
             t5 = V-X3
             t6 = Y1*J
             t7 = 2*t6
             t8 = r*t5
             Y3 = t8-t7
             t9 = Z1+H
             t10 = t9^2
             t11 = t10-Z1Z1
             Z3 = t11-HH
inputs       X1 X2 Y1 Y2 Z1

Doubling

I cannot easily tell which formula the crypto/elliptic P-256 uses.

I believe the cloudflare p384 implementation uses

https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b

id           g1p/shortw/jacobian-3/doubling/dbl-2001-b
tag          dbl-2001-b
class        g1p
shape        shortw
repr         jacobian-3
operation    doubling
url          https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-3.html#doubling-dbl-2001-b
cost         3M + 5S + 8add + 1*3 + 1*4 + 2*8
source       2001 Bernstein "A software implementation of NIST P-224"
appliesto    jacobian-3
compute      delta = Z1^2
             gamma = Y1^2
             beta = X1 gamma
             alpha = 3 (X1-delta) (X1+delta)
             X3 = alpha^2-8 beta
             Z3 = (Y1+Z1)^2-gamma-delta
             Y3 = alpha (4 beta-X3)-8 gamma^2
op3          delta = Z1^2
             gamma = Y1^2
             beta = X1*gamma
             t0 = X1-delta
             t1 = X1+delta
             t2 = t0*t1
             alpha = 3*t2
             t3 = alpha^2
             t4 = 8*beta
             X3 = t3-t4
             t5 = Y1+Z1
             t6 = t5^2
             t7 = t6-gamma
             Z3 = t7-delta
             t8 = 4*beta
             t9 = t8-X3
             t10 = gamma^2
             t11 = 8*t10
             t12 = alpha*t9
             Y3 = t12-t11
inputs       X1 Y1 Z1