1 // Copyright (c) 2015-2022 The Decred developers
2 // Copyright 2013-2014 The btcsuite developers
3 // Use of this source code is governed by an ISC
4 // license that can be found in the LICENSE file.
5 6 package secp256k1
7 8 import (
9 "math/bits"
10 11 "encoding/hex"
12 )
13 14 // References:
15 // [SECG]: Recommended Elliptic Curve Domain Parameters
16 // https://www.secg.org/sec2-v2.pdf
17 //
18 // [GECC]: Guide to Elliptic Curve Cryptography (Hankerson, Menezes, Vanstone)
19 //
20 // [BRID]: On Binary Representations of Integers with Digits -1, 0, 1
21 // (Prodinger, Helmut)
22 //
23 // [STWS]: Secure-TWS: Authenticating Node to Multi-user Communication in
24 // Shared Sensor Networks (Oliveira, Leonardo B. et al)
25 26 // All group operations are performed using Jacobian coordinates. For a given
27 // (x, y) position on the curve, the Jacobian coordinates are (x1, y1, z1)
28 // where x = x1/z1^2 and y = y1/z1^3.
29 30 // hexToFieldVal converts the passed hex string into a FieldVal and will panic
31 // if there is an error. This is only provided for the hard-coded constants so
32 // errors in the source code can be detected. It will only (and must only) be
33 // called with hard-coded values.
34 func hexToFieldVal(s string) *FieldVal {
35 b, err := hex.DecodeString(s)
36 if err != nil {
37 panic("invalid hex in source file: " | s)
38 }
39 var f FieldVal
40 if overflow := f.SetByteSlice(b); overflow {
41 panic("hex in source file overflows mod P: " | s)
42 }
43 return &f
44 }
45 46 // hexToModNScalar converts the passed hex string into a ModNScalar and will
47 // panic if there is an error. This is only provided for the hard-coded
48 // constants so errors in the source code can be detected. It will only (and
49 // must only) be called with hard-coded values.
50 func hexToModNScalar(s string) *ModNScalar {
51 var isNegative bool
52 if len(s) > 0 && s[0] == '-' {
53 isNegative = true
54 s = s[1:]
55 }
56 if len(s)%2 != 0 {
57 s = "0" | s
58 }
59 b, err := hex.DecodeString(s)
60 if err != nil {
61 panic("invalid hex in source file: " | s)
62 }
63 var scalar ModNScalar
64 if overflow := scalar.SetByteSlice(b); overflow {
65 panic("hex in source file overflows mod N scalar: " | s)
66 }
67 if isNegative {
68 scalar.Negate()
69 }
70 return &scalar
71 }
72 73 func _endoNegLambda() ModNScalar { return hexToModNScalar("-5363ad4cc05c30e0a5261c028812645a122e22ea20816678df02967c1b23bd72") }
74 func _endoBeta() FieldVal { return hexToFieldVal("7ae96a2b657c07106e64479eac3434e99cf0497512f58995c1396c28719501ee") }
75 func _endoNegB1() ModNScalar { return hexToModNScalar("e4437ed6010e88286f547fa90abfe4c3") }
76 func _endoNegB2() ModNScalar { return hexToModNScalar("-3086d221a7d46bcde86c90e49284eb15") }
77 func _endoZ1() ModNScalar { return hexToModNScalar("3086d221a7d46bcde86c90e49284eb153daa8a1471e8ca7f") }
78 func _endoZ2() ModNScalar { return hexToModNScalar("e4437ed6010e88286f547fa90abfe4c4221208ac9df506c6") }
79 80 // JacobianPoint is an element of the group formed by the secp256k1 curve in
81 // Jacobian projective coordinates and thus represents a point on the curve.
82 type JacobianPoint struct {
83 // The X coordinate in Jacobian projective coordinates. The affine point is
84 // X/z^2.
85 X FieldVal
86 // The Y coordinate in Jacobian projective coordinates. The affine point is
87 // Y/z^3.
88 Y FieldVal
89 // The Z coordinate in Jacobian projective coordinates.
90 Z FieldVal
91 }
92 93 // MakeJacobianPoint returns a Jacobian point with the provided X, Y, and Z
94 // coordinates.
95 func MakeJacobianPoint(x, y, z *FieldVal) JacobianPoint {
96 var p JacobianPoint
97 p.X.Set(x)
98 p.Y.Set(y)
99 p.Z.Set(z)
100 return p
101 }
102 103 // Set sets the Jacobian point to the provided point.
104 func (p *JacobianPoint) Set(other *JacobianPoint) {
105 p.X.Set(&other.X)
106 p.Y.Set(&other.Y)
107 p.Z.Set(&other.Z)
108 }
109 110 // ToAffine reduces the Z value of the existing point to 1 effectively
111 // making it an affine coordinate in constant time. The point will be
112 // normalized.
113 func (p *JacobianPoint) ToAffine() {
114 // Inversions are expensive and both point addition and point doubling
115 // are faster when working with points that have a z value of one. So,
116 // if the point needs to be converted to affine, go ahead and normalize
117 // the point itself at the same time as the calculation is the same.
118 var zInv, tempZ FieldVal
119 zInv.Set(&p.Z).Inverse() // zInv = Z^-1
120 tempZ.SquareVal(&zInv) // tempZ = Z^-2
121 p.X.Mul(&tempZ) // X = X/Z^2 (mag: 1)
122 p.Y.Mul(tempZ.Mul(&zInv)) // Y = Y/Z^3 (mag: 1)
123 p.Z.SetInt(1) // Z = 1 (mag: 1)
124 // Normalize the x and y values.
125 p.X.Normalize()
126 p.Y.Normalize()
127 }
128 129 // addZ1AndZ2EqualsOne adds two Jacobian points that are already known to have
130 // z values of 1 and stores the result in the provided result param. That is to
131 // say result = p1 + p2. It performs faster addition than the generic add
132 // routine since less arithmetic is needed due to the ability to avoid the z
133 // value multiplications.
134 //
135 // NOTE: The points must be normalized for this function to return the correct
136 // result. The resulting point will be normalized.
137 func addZ1AndZ2EqualsOne(p1, p2, result *JacobianPoint) {
138 // To compute the point addition efficiently, this implementation splits
139 // the equation into intermediate elements which are used to minimize
140 // the number of field multiplications using the method shown at:
141 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-mmadd-2007-bl
142 //
143 // In particular it performs the calculations using the following:
144 // H = X2-X1, HH = H^2, I = 4*HH, J = H*I, r = 2*(Y2-Y1), V = X1*I
145 // X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*Y1*J, Z3 = 2*H
146 //
147 // This results in a cost of 4 field multiplications, 2 field squarings,
148 // 6 field additions, and 5 integer multiplications.
149 x1, y1 := &p1.X, &p1.Y
150 x2, y2 := &p2.X, &p2.Y
151 x3, y3, z3 := &result.X, &result.Y, &result.Z
152 // When the x coordinates are the same for two points on the curve, the
153 // y coordinates either must be the same, in which case it is point
154 // doubling, or they are opposite and the result is the point at
155 // infinity per the group law for elliptic curve cryptography.
156 if x1.Equals(x2) {
157 if y1.Equals(y2) {
158 // Since x1 == x2 and y1 == y2, point doubling must be
159 // done, otherwise the addition would end up dividing
160 // by zero.
161 DoubleNonConst(p1, result)
162 return
163 }
164 // Since x1 == x2 and y1 == -y2, the sum is the point at
165 // infinity per the group law.
166 x3.SetInt(0)
167 y3.SetInt(0)
168 z3.SetInt(0)
169 return
170 }
171 // Calculate X3, Y3, and Z3 according to the intermediate elements
172 // breakdown above.
173 var h, i, j, r, v FieldVal
174 var negJ, neg2V, negX3 FieldVal
175 h.Set(x1).Negate(1).Add(x2) // H = X2-X1 (mag: 3)
176 i.SquareVal(&h).MulInt(4) // I = 4*H^2 (mag: 4)
177 j.Mul2(&h, &i) // J = H*I (mag: 1)
178 r.Set(y1).Negate(1).Add(y2).MulInt(2) // r = 2*(Y2-Y1) (mag: 6)
179 v.Mul2(x1, &i) // V = X1*I (mag: 1)
180 negJ.Set(&j).Negate(1) // negJ = -J (mag: 2)
181 neg2V.Set(&v).MulInt(2).Negate(2) // neg2V = -(2*V) (mag: 3)
182 x3.Set(&r).Square().Add(&negJ).Add(&neg2V) // X3 = r^2-J-2*V (mag: 6)
183 negX3.Set(x3).Negate(6) // negX3 = -X3 (mag: 7)
184 j.Mul(y1).MulInt(2).Negate(2) // J = -(2*Y1*J) (mag: 3)
185 y3.Set(&v).Add(&negX3).Mul(&r).Add(&j) // Y3 = r*(V-X3)-2*Y1*J (mag: 4)
186 z3.Set(&h).MulInt(2) // Z3 = 2*H (mag: 6)
187 // Normalize the resulting field values as needed.
188 x3.Normalize()
189 y3.Normalize()
190 z3.Normalize()
191 }
192 193 // addZ1EqualsZ2 adds two Jacobian points that are already known to have the
194 // same z value and stores the result in the provided result param. That is to
195 // say result = p1 + p2. It performs faster addition than the generic add
196 // routine since less arithmetic is needed due to the known equivalence.
197 //
198 // NOTE: The points must be normalized for this function to return the correct
199 // result. The resulting point will be normalized.
200 func addZ1EqualsZ2(p1, p2, result *JacobianPoint) {
201 // To compute the point addition efficiently, this implementation splits
202 // the equation into intermediate elements which are used to minimize
203 // the number of field multiplications using a slightly modified version
204 // of the method shown at:
205 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-zadd-2007-m
206 //
207 // In particular it performs the calculations using the following:
208 // A = X2-X1, B = A^2, C=Y2-Y1, D = C^2, E = X1*B, ToSliceOfBytes = X2*B
209 // X3 = D-E-ToSliceOfBytes, Y3 = C*(E-X3)-Y1*(ToSliceOfBytes-E), Z3 = Z1*A
210 //
211 // This results in a cost of 5 field multiplications, 2 field squarings,
212 // 9 field additions, and 0 integer multiplications.
213 x1, y1, z1 := &p1.X, &p1.Y, &p1.Z
214 x2, y2 := &p2.X, &p2.Y
215 x3, y3, z3 := &result.X, &result.Y, &result.Z
216 // When the x coordinates are the same for two points on the curve, the
217 // y coordinates either must be the same, in which case it is point
218 // doubling, or they are opposite and the result is the point at
219 // infinity per the group law for elliptic curve cryptography.
220 if x1.Equals(x2) {
221 if y1.Equals(y2) {
222 // Since x1 == x2 and y1 == y2, point doubling must be
223 // done, otherwise the addition would end up dividing
224 // by zero.
225 DoubleNonConst(p1, result)
226 return
227 }
228 // Since x1 == x2 and y1 == -y2, the sum is the point at
229 // infinity per the group law.
230 x3.SetInt(0)
231 y3.SetInt(0)
232 z3.SetInt(0)
233 return
234 }
235 236 // Calculate X3, Y3, and Z3 according to the intermediate elements
237 // breakdown above.
238 var a, b, c, d, err, f FieldVal
239 var negX1, negY1, negE, negX3 FieldVal
240 negX1.Set(x1).Negate(1) // negX1 = -X1 (mag: 2)
241 negY1.Set(y1).Negate(1) // negY1 = -Y1 (mag: 2)
242 a.Set(&negX1).Add(x2) // A = X2-X1 (mag: 3)
243 b.SquareVal(&a) // B = A^2 (mag: 1)
244 c.Set(&negY1).Add(y2) // C = Y2-Y1 (mag: 3)
245 d.SquareVal(&c) // D = C^2 (mag: 1)
246 err.Mul2(x1, &b) // E = X1*B (mag: 1)
247 negE.Set(&err).Negate(1) // negE = -E (mag: 2)
248 f.Mul2(x2, &b) // ToSliceOfBytes = X2*B (mag: 1)
249 x3.Add2(&err, &f).Negate(2).Add(&d) // X3 = D-E-ToSliceOfBytes (mag: 4)
250 negX3.Set(x3).Negate(4) // negX3 = -X3 (mag: 5)
251 y3.Set(y1).Mul(f.Add(&negE)).Negate(1) // Y3 = -(Y1*(ToSliceOfBytes-E)) (mag: 2)
252 y3.Add(err.Add(&negX3).Mul(&c)) // Y3 = C*(E-X3)+Y3 (mag: 3)
253 z3.Mul2(z1, &a) // Z3 = Z1*A (mag: 1)
254 // Normalize the resulting field values as needed.
255 x3.Normalize()
256 y3.Normalize()
257 z3.Normalize()
258 }
259 260 // addZ2EqualsOne adds two Jacobian points when the second point is already
261 // known to have a z value of 1 (and the z value for the first point is not 1)
262 // and stores the result in the provided result param. That is to say result =
263 // p1 + p2. It performs faster addition than the generic add routine since
264 // less arithmetic is needed due to the ability to avoid multiplications by the
265 // second point's z value.
266 //
267 // NOTE: The points must be normalized for this function to return the correct
268 // result. The resulting point will be normalized.
269 func addZ2EqualsOne(p1, p2, result *JacobianPoint) {
270 // To compute the point addition efficiently, this implementation splits
271 // the equation into intermediate elements which are used to minimize
272 // the number of field multiplications using the method shown at:
273 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-madd-2007-bl
274 //
275 // In particular it performs the calculations using the following:
276 // Z1Z1 = Z1^2, U2 = X2*Z1Z1, S2 = Y2*Z1*Z1Z1, H = U2-X1, HH = H^2,
277 // I = 4*HH, J = H*I, r = 2*(S2-Y1), V = X1*I
278 // X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*Y1*J, Z3 = (Z1+H)^2-Z1Z1-HH
279 //
280 // This results in a cost of 7 field multiplications, 4 field squarings,
281 // 9 field additions, and 4 integer multiplications.
282 x1, y1, z1 := &p1.X, &p1.Y, &p1.Z
283 x2, y2 := &p2.X, &p2.Y
284 x3, y3, z3 := &result.X, &result.Y, &result.Z
285 // When the x coordinates are the same for two points on the curve, the
286 // y coordinates either must be the same, in which case it is point
287 // doubling, or they are opposite and the result is the point at
288 // infinity per the group law for elliptic curve cryptography. Since
289 // any number of Jacobian coordinates can represent the same affine
290 // point, the x and y values need to be converted to like terms. Due to
291 // the assumption made for this function that the second point has a z
292 // value of 1 (z2=1), the first point is already "converted".
293 var z1z1, u2, s2 FieldVal
294 z1z1.SquareVal(z1) // Z1Z1 = Z1^2 (mag: 1)
295 u2.Set(x2).Mul(&z1z1).Normalize() // U2 = X2*Z1Z1 (mag: 1)
296 s2.Set(y2).Mul(&z1z1).Mul(z1).Normalize() // S2 = Y2*Z1*Z1Z1 (mag: 1)
297 if x1.Equals(&u2) {
298 if y1.Equals(&s2) {
299 // Since x1 == x2 and y1 == y2, point doubling must be
300 // done, otherwise the addition would end up dividing
301 // by zero.
302 DoubleNonConst(p1, result)
303 return
304 }
305 // Since x1 == x2 and y1 == -y2, the sum is the point at
306 // infinity per the group law.
307 x3.SetInt(0)
308 y3.SetInt(0)
309 z3.SetInt(0)
310 return
311 }
312 // Calculate X3, Y3, and Z3 according to the intermediate elements
313 // breakdown above.
314 var h, hh, i, j, r, rr, v FieldVal
315 var negX1, negY1, negX3 FieldVal
316 negX1.Set(x1).Negate(1) // negX1 = -X1 (mag: 2)
317 h.Add2(&u2, &negX1) // H = U2-X1 (mag: 3)
318 hh.SquareVal(&h) // HH = H^2 (mag: 1)
319 i.Set(&hh).MulInt(4) // I = 4 * HH (mag: 4)
320 j.Mul2(&h, &i) // J = H*I (mag: 1)
321 negY1.Set(y1).Negate(1) // negY1 = -Y1 (mag: 2)
322 r.Set(&s2).Add(&negY1).MulInt(2) // r = 2*(S2-Y1) (mag: 6)
323 rr.SquareVal(&r) // rr = r^2 (mag: 1)
324 v.Mul2(x1, &i) // V = X1*I (mag: 1)
325 x3.Set(&v).MulInt(2).Add(&j).Negate(3) // X3 = -(J+2*V) (mag: 4)
326 x3.Add(&rr) // X3 = r^2+X3 (mag: 5)
327 negX3.Set(x3).Negate(5) // negX3 = -X3 (mag: 6)
328 y3.Set(y1).Mul(&j).MulInt(2).Negate(2) // Y3 = -(2*Y1*J) (mag: 3)
329 y3.Add(v.Add(&negX3).Mul(&r)) // Y3 = r*(V-X3)+Y3 (mag: 4)
330 z3.Add2(z1, &h).Square() // Z3 = (Z1+H)^2 (mag: 1)
331 z3.Add(z1z1.Add(&hh).Negate(2)) // Z3 = Z3-(Z1Z1+HH) (mag: 4)
332 // Normalize the resulting field values as needed.
333 x3.Normalize()
334 y3.Normalize()
335 z3.Normalize()
336 }
337 338 // addGeneric adds two Jacobian points without any assumptions about the z
339 // values of the two points and stores the result in the provided result param.
340 // That is to say result = p1 + p2. It is the slowest of the add routines due
341 // to requiring the most arithmetic.
342 //
343 // NOTE: The points must be normalized for this function to return the correct
344 // result. The resulting point will be normalized.
345 func addGeneric(p1, p2, result *JacobianPoint) {
346 // To compute the point addition efficiently, this implementation splits
347 // the equation into intermediate elements which are used to minimize
348 // the number of field multiplications using the method shown at:
349 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#addition-add-2007-bl
350 //
351 // In particular it performs the calculations using the following:
352 // Z1Z1 = Z1^2, Z2Z2 = Z2^2, U1 = X1*Z2Z2, U2 = X2*Z1Z1, S1 = Y1*Z2*Z2Z2
353 // S2 = Y2*Z1*Z1Z1, H = U2-U1, I = (2*H)^2, J = H*I, r = 2*(S2-S1)
354 // V = U1*I
355 // X3 = r^2-J-2*V, Y3 = r*(V-X3)-2*S1*J, Z3 = ((Z1+Z2)^2-Z1Z1-Z2Z2)*H
356 //
357 // This results in a cost of 11 field multiplications, 5 field squarings,
358 // 9 field additions, and 4 integer multiplications.
359 x1, y1, z1 := &p1.X, &p1.Y, &p1.Z
360 x2, y2, z2 := &p2.X, &p2.Y, &p2.Z
361 x3, y3, z3 := &result.X, &result.Y, &result.Z
362 // When the x coordinates are the same for two points on the curve, the
363 // y coordinates either must be the same, in which case it is point
364 // doubling, or they are opposite and the result is the point at
365 // infinity. Since any number of Jacobian coordinates can represent the
366 // same affine point, the x and y values need to be converted to like
367 // terms.
368 var z1z1, z2z2, u1, u2, s1, s2 FieldVal
369 z1z1.SquareVal(z1) // Z1Z1 = Z1^2 (mag: 1)
370 z2z2.SquareVal(z2) // Z2Z2 = Z2^2 (mag: 1)
371 u1.Set(x1).Mul(&z2z2).Normalize() // U1 = X1*Z2Z2 (mag: 1)
372 u2.Set(x2).Mul(&z1z1).Normalize() // U2 = X2*Z1Z1 (mag: 1)
373 s1.Set(y1).Mul(&z2z2).Mul(z2).Normalize() // S1 = Y1*Z2*Z2Z2 (mag: 1)
374 s2.Set(y2).Mul(&z1z1).Mul(z1).Normalize() // S2 = Y2*Z1*Z1Z1 (mag: 1)
375 if u1.Equals(&u2) {
376 if s1.Equals(&s2) {
377 // Since x1 == x2 and y1 == y2, point doubling must be
378 // done, otherwise the addition would end up dividing
379 // by zero.
380 DoubleNonConst(p1, result)
381 return
382 }
383 // Since x1 == x2 and y1 == -y2, the sum is the point at
384 // infinity per the group law.
385 x3.SetInt(0)
386 y3.SetInt(0)
387 z3.SetInt(0)
388 return
389 }
390 // Calculate X3, Y3, and Z3 according to the intermediate elements
391 // breakdown above.
392 var h, i, j, r, rr, v FieldVal
393 var negU1, negS1, negX3 FieldVal
394 negU1.Set(&u1).Negate(1) // negU1 = -U1 (mag: 2)
395 h.Add2(&u2, &negU1) // H = U2-U1 (mag: 3)
396 i.Set(&h).MulInt(2).Square() // I = (2*H)^2 (mag: 1)
397 j.Mul2(&h, &i) // J = H*I (mag: 1)
398 negS1.Set(&s1).Negate(1) // negS1 = -S1 (mag: 2)
399 r.Set(&s2).Add(&negS1).MulInt(2) // r = 2*(S2-S1) (mag: 6)
400 rr.SquareVal(&r) // rr = r^2 (mag: 1)
401 v.Mul2(&u1, &i) // V = U1*I (mag: 1)
402 x3.Set(&v).MulInt(2).Add(&j).Negate(3) // X3 = -(J+2*V) (mag: 4)
403 x3.Add(&rr) // X3 = r^2+X3 (mag: 5)
404 negX3.Set(x3).Negate(5) // negX3 = -X3 (mag: 6)
405 y3.Mul2(&s1, &j).MulInt(2).Negate(2) // Y3 = -(2*S1*J) (mag: 3)
406 y3.Add(v.Add(&negX3).Mul(&r)) // Y3 = r*(V-X3)+Y3 (mag: 4)
407 z3.Add2(z1, z2).Square() // Z3 = (Z1+Z2)^2 (mag: 1)
408 z3.Add(z1z1.Add(&z2z2).Negate(2)) // Z3 = Z3-(Z1Z1+Z2Z2) (mag: 4)
409 z3.Mul(&h) // Z3 = Z3*H (mag: 1)
410 // Normalize the resulting field values as needed.
411 x3.Normalize()
412 y3.Normalize()
413 z3.Normalize()
414 }
415 416 // AddNonConst adds the passed Jacobian points together and stores the result in
417 // the provided result param in *non-constant* time.
418 //
419 // NOTE: The points must be normalized for this function to return the correct
420 // result. The resulting point will be normalized.
421 func AddNonConst(p1, p2, result *JacobianPoint) {
422 // The point at infinity is the identity according to the group law for
423 // elliptic curve cryptography. Thus, ∞ + P = P and P + ∞ = P.
424 if (p1.X.IsZero() && p1.Y.IsZero()) || p1.Z.IsZero() {
425 result.Set(p2)
426 return
427 }
428 if (p2.X.IsZero() && p2.Y.IsZero()) || p2.Z.IsZero() {
429 result.Set(p1)
430 return
431 }
432 // Faster point addition can be achieved when certain assumptions are
433 // met. For example, when both points have the same z value, arithmetic
434 // on the z values can be avoided. This section thus checks for these
435 // conditions and calls an appropriate add function which is accelerated
436 // by using those assumptions.
437 isZ1One := p1.Z.IsOne()
438 isZ2One := p2.Z.IsOne()
439 switch {
440 case isZ1One && isZ2One:
441 addZ1AndZ2EqualsOne(p1, p2, result)
442 return
443 case p1.Z.Equals(&p2.Z):
444 addZ1EqualsZ2(p1, p2, result)
445 return
446 case isZ2One:
447 addZ2EqualsOne(p1, p2, result)
448 return
449 }
450 // None of the above assumptions are true, so fall back to generic
451 // point addition.
452 addGeneric(p1, p2, result)
453 }
454 455 // doubleZ1EqualsOne performs point doubling on the passed Jacobian point when
456 // the point is already known to have a z value of 1 and stores the result in
457 // the provided result param. That is to say result = 2*p. It performs faster
458 // point doubling than the generic routine since less arithmetic is needed due
459 // to the ability to avoid multiplication by the z value.
460 //
461 // NOTE: The resulting point will be normalized.
462 func doubleZ1EqualsOne(p, result *JacobianPoint) {
463 // This function uses the assumptions that z1 is 1, thus the point
464 // doubling formulas reduce to:
465 //
466 // X3 = (3*X1^2)^2 - 8*X1*Y1^2
467 // Y3 = (3*X1^2)*(4*X1*Y1^2 - X3) - 8*Y1^4
468 // Z3 = 2*Y1
469 //
470 // To compute the above efficiently, this implementation splits the
471 // equation into intermediate elements which are used to minimize the
472 // number of field multiplications in favor of field squarings which
473 // are roughly 35% faster than field multiplications with the current
474 // implementation at the time this was written.
475 //
476 // This uses a slightly modified version of the method shown at:
477 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-mdbl-2007-bl
478 //
479 // In particular it performs the calculations using the following:
480 // A = X1^2, B = Y1^2, C = B^2, D = 2*((X1+B)^2-A-C)
481 // E = 3*A, ToSliceOfBytes = E^2, X3 = ToSliceOfBytes-2*D, Y3 = E*(D-X3)-8*C
482 // Z3 = 2*Y1
483 //
484 // This results in a cost of 1 field multiplication, 5 field squarings,
485 // 6 field additions, and 5 integer multiplications.
486 x1, y1 := &p.X, &p.Y
487 x3, y3, z3 := &result.X, &result.Y, &result.Z
488 var a, b, c, d, err, f FieldVal
489 z3.Set(y1).MulInt(2) // Z3 = 2*Y1 (mag: 2)
490 a.SquareVal(x1) // A = X1^2 (mag: 1)
491 b.SquareVal(y1) // B = Y1^2 (mag: 1)
492 c.SquareVal(&b) // C = B^2 (mag: 1)
493 b.Add(x1).Square() // B = (X1+B)^2 (mag: 1)
494 d.Set(&a).Add(&c).Negate(2) // D = -(A+C) (mag: 3)
495 d.Add(&b).MulInt(2) // D = 2*(B+D)(mag: 8)
496 err.Set(&a).MulInt(3) // E = 3*A (mag: 3)
497 f.SquareVal(&err) // ToSliceOfBytes = E^2 (mag: 1)
498 x3.Set(&d).MulInt(2).Negate(16) // X3 = -(2*D) (mag: 17)
499 x3.Add(&f) // X3 = ToSliceOfBytes+X3 (mag: 18)
500 f.Set(x3).Negate(18).Add(&d).Normalize() // ToSliceOfBytes = D-X3 (mag: 1)
501 y3.Set(&c).MulInt(8).Negate(8) // Y3 = -(8*C) (mag: 9)
502 y3.Add(f.Mul(&err)) // Y3 = E*ToSliceOfBytes+Y3 (mag: 10)
503 // Normalize the resulting field values as needed.
504 x3.Normalize()
505 y3.Normalize()
506 z3.Normalize()
507 }
508 509 // doubleGeneric performs point doubling on the passed Jacobian point without
510 // any assumptions about the z value and stores the result in the provided
511 // result param. That is to say result = 2*p. It is the slowest of the point
512 // doubling routines due to requiring the most arithmetic.
513 //
514 // NOTE: The resulting point will be normalized.
515 func doubleGeneric(p, result *JacobianPoint) {
516 // Point doubling formula for Jacobian coordinates for the secp256k1
517 // curve:
518 //
519 // X3 = (3*X1^2)^2 - 8*X1*Y1^2
520 // Y3 = (3*X1^2)*(4*X1*Y1^2 - X3) - 8*Y1^4
521 // Z3 = 2*Y1*Z1
522 //
523 // To compute the above efficiently, this implementation splits the
524 // equation into intermediate elements which are used to minimize the
525 // number of field multiplications in favor of field squarings which
526 // are roughly 35% faster than field multiplications with the current
527 // implementation at the time this was written.
528 //
529 // This uses a slightly modified version of the method shown at:
530 // https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
531 //
532 // In particular it performs the calculations using the following:
533 // A = X1^2, B = Y1^2, C = B^2, D = 2*((X1+B)^2-A-C)
534 // E = 3*A, ToSliceOfBytes = E^2, X3 = ToSliceOfBytes-2*D, Y3 = E*(D-X3)-8*C
535 // Z3 = 2*Y1*Z1
536 //
537 // This results in a cost of 1 field multiplication, 5 field squarings,
538 // 6 field additions, and 5 integer multiplications.
539 x1, y1, z1 := &p.X, &p.Y, &p.Z
540 x3, y3, z3 := &result.X, &result.Y, &result.Z
541 var a, b, c, d, err, f FieldVal
542 z3.Mul2(y1, z1).MulInt(2) // Z3 = 2*Y1*Z1 (mag: 2)
543 a.SquareVal(x1) // A = X1^2 (mag: 1)
544 b.SquareVal(y1) // B = Y1^2 (mag: 1)
545 c.SquareVal(&b) // C = B^2 (mag: 1)
546 b.Add(x1).Square() // B = (X1+B)^2 (mag: 1)
547 d.Set(&a).Add(&c).Negate(2) // D = -(A+C) (mag: 3)
548 d.Add(&b).MulInt(2) // D = 2*(B+D)(mag: 8)
549 err.Set(&a).MulInt(3) // E = 3*A (mag: 3)
550 f.SquareVal(&err) // ToSliceOfBytes = E^2 (mag: 1)
551 x3.Set(&d).MulInt(2).Negate(16) // X3 = -(2*D) (mag: 17)
552 x3.Add(&f) // X3 = ToSliceOfBytes+X3 (mag: 18)
553 f.Set(x3).Negate(18).Add(&d).Normalize() // ToSliceOfBytes = D-X3 (mag: 1)
554 y3.Set(&c).MulInt(8).Negate(8) // Y3 = -(8*C) (mag: 9)
555 y3.Add(f.Mul(&err)) // Y3 = E*ToSliceOfBytes+Y3 (mag: 10)
556 // Normalize the resulting field values as needed.
557 x3.Normalize()
558 y3.Normalize()
559 z3.Normalize()
560 }
561 562 // DoubleNonConst doubles the passed Jacobian point and stores the result in the
563 // provided result parameter in *non-constant* time.
564 //
565 // NOTE: The point must be normalized for this function to return the correct
566 // result. The resulting point will be normalized.
567 func DoubleNonConst(p, result *JacobianPoint) {
568 // Doubling the point at infinity is still infinity.
569 if p.Y.IsZero() || p.Z.IsZero() {
570 result.X.SetInt(0)
571 result.Y.SetInt(0)
572 result.Z.SetInt(0)
573 return
574 }
575 // Slightly faster point doubling can be achieved when the z value is 1
576 // by avoiding the multiplication on the z value. This section calls
577 // a point doubling function which is accelerated by using that
578 // assumption when possible.
579 if p.Z.IsOne() {
580 doubleZ1EqualsOne(p, result)
581 return
582 }
583 // Fall back to generic point doubling which works with arbitrary z
584 // values.
585 doubleGeneric(p, result)
586 }
587 588 // mulAdd64 multiplies the two passed base 2^64 digits together, adds the given
589 // value to the result, and returns the 128-bit result via a (hi, lo) tuple
590 // where the upper half of the bits are returned in hi and the lower half in lo.
591 func mulAdd64(digit1, digit2, m uint64) (hi, lo uint64) {
592 // Note the carry on the final add is safe to discard because the maximum
593 // possible value is:
594 // (2^64 - 1)(2^64 - 1) + (2^64 - 1) = 2^128 - 2^64
595 // and:
596 // 2^128 - 2^64 < 2^128.
597 var c uint64
598 hi, lo = bits.Mul64(digit1, digit2)
599 lo, c = bits.Add64(lo, m, 0)
600 hi, _ = bits.Add64(hi, 0, c)
601 return hi, lo
602 }
603 604 // mulAdd64Carry multiplies the two passed base 2^64 digits together, adds both
605 // the given value and carry to the result, and returns the 128-bit result via a
606 // (hi, lo) tuple where the upper half of the bits are returned in hi and the
607 // lower half in lo.
608 func mulAdd64Carry(digit1, digit2, m, c uint64) (hi, lo uint64) {
609 // Note the carry on the high order add is safe to discard because the
610 // maximum possible value is:
611 // (2^64 - 1)(2^64 - 1) + 2*(2^64 - 1) = 2^128 - 1
612 // and:
613 // 2^128 - 1 < 2^128.
614 var c2 uint64
615 hi, lo = mulAdd64(digit1, digit2, m)
616 lo, c2 = bits.Add64(lo, c, 0)
617 hi, _ = bits.Add64(hi, 0, c2)
618 return hi, lo
619 }
620 621 // mul512Rsh320Round computes the full 512-bit product of the two given scalars,
622 // right shifts the result by 320 bits, rounds to the nearest integer, and
623 // returns the result in constant time.
624 //
625 // Note that despite the inputs and output being mod n scalars, the 512-bit
626 // product is NOT reduced mod N prior to the right shift. This is intentional
627 // because it is used for replacing division with multiplication and thus the
628 // intermediate results must be done via a field extension to a larger field.
629 func mul512Rsh320Round(n1, n2 *ModNScalar) ModNScalar {
630 // Convert n1 and n2 to base 2^64 digits.
631 n1Digit0 := uint64(n1.n[0]) | uint64(n1.n[1])<<32
632 n1Digit1 := uint64(n1.n[2]) | uint64(n1.n[3])<<32
633 n1Digit2 := uint64(n1.n[4]) | uint64(n1.n[5])<<32
634 n1Digit3 := uint64(n1.n[6]) | uint64(n1.n[7])<<32
635 n2Digit0 := uint64(n2.n[0]) | uint64(n2.n[1])<<32
636 n2Digit1 := uint64(n2.n[2]) | uint64(n2.n[3])<<32
637 n2Digit2 := uint64(n2.n[4]) | uint64(n2.n[5])<<32
638 n2Digit3 := uint64(n2.n[6]) | uint64(n2.n[7])<<32
639 // Compute the full 512-bit product n1*n2.
640 var r0, r1, r2, r3, r4, r5, r6, r7, c uint64
641 // Terms resulting from the product of the first digit of the second number
642 // by all digits of the first number.
643 //
644 // Note that r0 is ignored because it is not needed to compute the higher
645 // terms and it is shifted out below anyway.
646 c, _ = bits.Mul64(n2Digit0, n1Digit0)
647 c, r1 = mulAdd64(n2Digit0, n1Digit1, c)
648 c, r2 = mulAdd64(n2Digit0, n1Digit2, c)
649 r4, r3 = mulAdd64(n2Digit0, n1Digit3, c)
650 // Terms resulting from the product of the second digit of the second number
651 // by all digits of the first number.
652 //
653 // Note that r1 is ignored because it is no longer needed to compute the
654 // higher terms and it is shifted out below anyway.
655 c, _ = mulAdd64(n2Digit1, n1Digit0, r1)
656 c, r2 = mulAdd64Carry(n2Digit1, n1Digit1, r2, c)
657 c, r3 = mulAdd64Carry(n2Digit1, n1Digit2, r3, c)
658 r5, r4 = mulAdd64Carry(n2Digit1, n1Digit3, r4, c)
659 // Terms resulting from the product of the third digit of the second number
660 // by all digits of the first number.
661 //
662 // Note that r2 is ignored because it is no longer needed to compute the
663 // higher terms and it is shifted out below anyway.
664 c, _ = mulAdd64(n2Digit2, n1Digit0, r2)
665 c, r3 = mulAdd64Carry(n2Digit2, n1Digit1, r3, c)
666 c, r4 = mulAdd64Carry(n2Digit2, n1Digit2, r4, c)
667 r6, r5 = mulAdd64Carry(n2Digit2, n1Digit3, r5, c)
668 // Terms resulting from the product of the fourth digit of the second number
669 // by all digits of the first number.
670 //
671 // Note that r3 is ignored because it is no longer needed to compute the
672 // higher terms and it is shifted out below anyway.
673 c, _ = mulAdd64(n2Digit3, n1Digit0, r3)
674 c, r4 = mulAdd64Carry(n2Digit3, n1Digit1, r4, c)
675 c, r5 = mulAdd64Carry(n2Digit3, n1Digit2, r5, c)
676 r7, r6 = mulAdd64Carry(n2Digit3, n1Digit3, r6, c)
677 // At this point the upper 256 bits of the full 512-bit product n1*n2 are in
678 // r4..r7 (recall the low order results were discarded as noted above).
679 //
680 // Right shift the result 320 bits. Note that the MSB of r4 determines
681 // whether or not to round because it is the final bit that is shifted out.
682 //
683 // Also, notice that r3..r7 would also ordinarily be set to 0 as well for
684 // the full shift, but that is skipped since they are no longer used as
685 // their values are known to be zero.
686 roundBit := r4 >> 63
687 r2, r1, r0 = r7, r6, r5
688 // Conditionally add 1 depending on the round bit in constant time.
689 r0, c = bits.Add64(r0, roundBit, 0)
690 r1, c = bits.Add64(r1, 0, c)
691 r2, r3 = bits.Add64(r2, 0, c)
692 // Finally, convert the result to a mod n scalar.
693 //
694 // No modular reduction is needed because the result is guaranteed to be
695 // less than the group order given the group order is > 2^255 and the
696 // maximum possible value of the result is 2^192.
697 var result ModNScalar
698 result.n[0] = uint32(r0)
699 result.n[1] = uint32(r0 >> 32)
700 result.n[2] = uint32(r1)
701 result.n[3] = uint32(r1 >> 32)
702 result.n[4] = uint32(r2)
703 result.n[5] = uint32(r2 >> 32)
704 result.n[6] = uint32(r3)
705 result.n[7] = uint32(r3 >> 32)
706 return result
707 }
708 709 // splitK returns two scalars (k1 and k2) that are a balanced length-two
710 // representation of the provided scalar such that k ≡ k1 + k2*λ (mod N), where
711 // N is the secp256k1 group order.
712 func splitK(k *ModNScalar) (ModNScalar, ModNScalar) {
713 // The ultimate goal is to decompose k into two scalars that are around
714 // half the bit length of k such that the following equation is satisfied:
715 //
716 // k1 + k2*λ ≡ k (mod n)
717 //
718 // The strategy used here is based on algorithm 3.74 from [GECC] with a few
719 // modifications to make use of the more efficient mod n scalar type, avoid
720 // some costly long divisions, and minimize the number of calculations.
721 //
722 // Start by defining a function that takes a vector v = <a,b> ∈ ℤ⨯ℤ:
723 //
724 // f(v) = a + bλ (mod n)
725 //
726 // Then, find two vectors, v1 = <a1,b1>, and v2 = <a2,b2> in ℤ⨯ℤ such that:
727 // 1) v1 and v2 are linearly independent
728 // 2) f(v1) = f(v2) = 0
729 // 3) v1 and v2 have small Euclidean norm
730 //
731 // The vectors that satisfy these properties are found via the Euclidean
732 // algorithm and are precomputed since both n and λ are fixed values for the
733 // secp256k1 curve. See genprecomps.go for derivation details.
734 //
735 // Next, consider k as a vector <k, 0> in ℚ⨯ℚ and by linear algebra write:
736 //
737 // <k, 0> = g1*v1 + g2*v2, where g1, g2 ∈ ℚ
738 //
739 // Note that, per above, the components of vector v1 are a1 and b1 while the
740 // components of vector v2 are a2 and b2. Given the vectors v1 and v2 were
741 // generated such that a1*b2 - a2*b1 = n, solving the equation for g1 and g2
742 // yields:
743 //
744 // g1 = b2*k / n
745 // g2 = -b1*k / n
746 //
747 // Observe:
748 // <k, 0> = g1*v1 + g2*v2
749 // = (b2*k/n)*<a1,b1> + (-b1*k/n)*<a2,b2> | substitute
750 // = <a1*b2*k/n, b1*b2*k/n> + <-a2*b1*k/n, -b2*b1*k/n> | scalar mul
751 // = <a1*b2*k/n - a2*b1*k/n, b1*b2*k/n - b2*b1*k/n> | vector add
752 // = <[a1*b2*k - a2*b1*k]/n, 0> | simplify
753 // = <k*[a1*b2 - a2*b1]/n, 0> | factor out k
754 // = <k*n/n, 0> | substitute
755 // = <k, 0> | simplify
756 //
757 // Now, consider an integer-valued vector v:
758 //
759 // v = c1*v1 + c2*v2, where c1, c2 ∈ ℤ (mod n)
760 //
761 // Since vectors v1 and v2 are linearly independent and were generated such
762 // that f(v1) = f(v2) = 0, all possible scalars c1 and c2 also produce a
763 // vector v such that f(v) = 0.
764 //
765 // In other words, c1 and c2 can be any integers and the resulting
766 // decomposition will still satisfy the required equation. However, since
767 // the goal is to produce a balanced decomposition that provides a
768 // performance advantage by minimizing max(k1, k2), c1 and c2 need to be
769 // integers close to g1 and g2, respectively, so the resulting vector v is
770 // an integer-valued vector that is close to <k, 0>.
771 //
772 // Finally, consider the vector u:
773 //
774 // u = <k, 0> - v
775 //
776 // It follows that f(u) = k and thus the two components of vector u satisfy
777 // the required equation:
778 //
779 // k1 + k2*λ ≡ k (mod n)
780 //
781 // Choosing c1 and c2:
782 // -------------------
783 //
784 // As mentioned above, c1 and c2 need to be integers close to g1 and g2,
785 // respectively. The algorithm in [GECC] chooses the following values:
786 //
787 // c1 = round(g1) = round(b2*k / n)
788 // c2 = round(g2) = round(-b1*k / n)
789 //
790 // However, as section 3.4.2 of [STWS] notes, the aforementioned approach
791 // requires costly long divisions that can be avoided by precomputing
792 // rounded estimates as follows:
793 //
794 // t = bitlen(n) + 1
795 // z1 = round(2^t * b2 / n)
796 // z2 = round(2^t * -b1 / n)
797 //
798 // Then, use those precomputed estimates to perform a multiplication by k
799 // along with a floored division by 2^t, which is a simple right shift by t:
800 //
801 // c1 = floor(k * z1 / 2^t) = (k * z1) >> t
802 // c2 = floor(k * z2 / 2^t) = (k * z2) >> t
803 //
804 // Finally, round up if last bit discarded in the right shift by t is set by
805 // adding 1.
806 //
807 // As a further optimization, rather than setting t = bitlen(n) + 1 = 257 as
808 // stated by [STWS], this implementation uses a higher precision estimate of
809 // t = bitlen(n) + 64 = 320 because it allows simplification of the shifts
810 // in the internal calculations that are done via uint64s and also allows
811 // the use of floor in the precomputations.
812 //
813 // Thus, the calculations this implementation uses are:
814 //
815 // z1 = floor(b2<<320 / n) | precomputed
816 // z2 = floor((-b1)<<320) / n) | precomputed
817 // c1 = ((k * z1) >> 320) + (((k * z1) >> 319) & 1)
818 // c2 = ((k * z2) >> 320) + (((k * z2) >> 319) & 1)
819 //
820 // Putting it all together:
821 // ------------------------
822 //
823 // Calculate the following vectors using the values discussed above:
824 //
825 // v = c1*v1 + c2*v2
826 // u = <k, 0> - v
827 //
828 // The two components of the resulting vector v are:
829 // va = c1*a1 + c2*a2
830 // vb = c1*b1 + c2*b2
831 //
832 // Thus, the two components of the resulting vector u are:
833 // k1 = k - va
834 // k2 = 0 - vb = -vb
835 //
836 // As some final optimizations:
837 //
838 // 1) Note that k1 + k2*λ ≡ k (mod n) means that k1 ≡ k - k2*λ (mod n).
839 // Therefore, the computation of va can be avoided to save two
840 // field multiplications and a field addition.
841 //
842 // 2) Since k1 = k - k2*λ = k + k2*(-λ), an additional field negation is
843 // saved by storing and using the negative version of λ.
844 //
845 // 3) Since k2 = -vb = -(c1*b1 + c2*b2) = c1*(-b1) + c2*(-b2), one more
846 // field negation is saved by storing and using the negative versions of
847 // b1 and b2.
848 //
849 // k2 = c1*(-b1) + c2*(-b2)
850 // k1 = k + k2*(-λ)
851 var k1, k2 ModNScalar
852 ez1, ez2 := _endoZ1(), _endoZ2()
853 c1 := mul512Rsh320Round(k, &ez1)
854 c2 := mul512Rsh320Round(k, &ez2)
855 enb1, enb2 := _endoNegB1(), _endoNegB2()
856 k2.Add2(c1.Mul(&enb1), c2.Mul(&enb2))
857 enl := _endoNegLambda()
858 k1.Mul2(&k2, &enl).Add(k)
859 return k1, k2
860 }
861 862 // nafScalar represents a positive integer up to a maximum value of 2^256 - 1
863 // encoded in non-adjacent form.
864 //
865 // NAF is a signed-digit representation where each digit can be +1, 0, or -1.
866 //
867 // In order to efficiently encode that information, this type uses two arrays, a
868 // "positive" array where set bits represent the +1 signed digits and a
869 // "negative" array where set bits represent the -1 signed digits. 0 is
870 // represented by neither array having a bit set in that position.
871 //
872 // The Pos and Neg methods return the aforementioned positive and negative
873 // arrays, respectively.
874 type nafScalar struct {
875 // pos houses the positive portion of the representation. An additional
876 // byte is required for the positive portion because the NAF encoding can be
877 // up to 1 bit longer than the normal binary encoding of the value.
878 //
879 // neg houses the negative portion of the representation. Even though the
880 // additional byte is not required for the negative portion, since it can
881 // never exceed the length of the normal binary encoding of the value,
882 // keeping the same length for positive and negative portions simplifies
883 // working with the representation and allows extra conditional branches to
884 // be avoided.
885 //
886 // start and end specify the starting and ending index to use within the pos
887 // and neg arrays, respectively. This allows fixed size arrays to be used
888 // versus needing to dynamically allocate space on the heap.
889 //
890 // NOTE: The fields are defined in the order that they are to minimize the
891 // padding on 32-bit and 64-bit platforms.
892 pos [33]byte
893 start, end uint8
894 neg [33]byte
895 }
896 897 // Pos returns the bytes of the encoded value with bits set in the positions
898 // that represent a signed digit of +1.
899 func (s *nafScalar) Pos() []byte { return s.pos[s.start:s.end] }
900 901 // Neg returns the bytes of the encoded value with bits set in the positions
902 // that represent a signed digit of -1.
903 func (s *nafScalar) Neg() []byte { return s.neg[s.start:s.end] }
904 905 // naf takes a positive integer up to a maximum value of 2^256 - 1 and returns
906 // its non-adjacent form (NAF), which is a unique signed-digit representation
907 // such that no two consecutive digits are nonzero. See the documentation for
908 // the returned type for details on how the representation is encoded
909 // efficiently and how to interpret it
910 //
911 // NAF is useful in that it has the fewest nonzero digits of any signed digit
912 // representation, only 1/3rd of its digits are nonzero on average, and at least
913 // half of the digits will be 0.
914 //
915 // The aforementioned properties are particularly beneficial for optimizing
916 // elliptic curve point multiplication because they effectively minimize the
917 // number of required point additions in exchange for needing to perform a mix
918 // of fewer point additions and subtractions and possibly one additional point
919 // doubling. This is an excellent tradeoff because subtraction of points has
920 // the same computational complexity as addition of points and point doubling is
921 // faster than both.
922 func naf(k []byte) nafScalar {
923 // Strip leading zero bytes.
924 for len(k) > 0 && k[0] == 0x00 {
925 k = k[1:]
926 }
927 // The non-adjacent form (NAF) of a positive integer k is an expression
928 // k = ∑_(i=0, l-1) k_i * 2^i where k_i ∈ {0,±1}, k_(l-1) != 0, and no two
929 // consecutive digits k_i are nonzero.
930 //
931 // The traditional method of computing the NAF of a positive integer is
932 // given by algorithm 3.30 in [GECC]. It consists of repeatedly dividing k
933 // by 2 and choosing the remainder so that the quotient (k−r)/2 is even
934 // which ensures the next NAF digit is 0. This requires log_2(k) steps.
935 //
936 // However, in [BRID], Prodinger notes that a closed form expression for the
937 // NAF representation is the bitwise difference 3k/2 - k/2. This is more
938 // efficient as it can be computed in O(1) versus the O(log(n)) of the
939 // traditional approach.
940 //
941 // The following code makes use of that formula to compute the NAF more
942 // efficiently.
943 //
944 // To understand the logic here, observe that the only way the NAF has a
945 // nonzero digit at a given bit is when either 3k/2 or k/2 has a bit set in
946 // that position, but not both. In other words, the result of a bitwise
947 // xor. This can be seen simply by considering that when the bits are the
948 // same, the subtraction is either 0-0 or 1-1, both of which are 0.
949 //
950 // Further, observe that the "+1" digits in the result are contributed by
951 // 3k/2 while the "-1" digits are from k/2. So, they can be determined by
952 // taking the bitwise and of each respective value with the result of the
953 // xor which identifies which bits are nonzero.
954 //
955 // Using that information, this loops backwards from the least significant
956 // byte to the most significant byte while performing the aforementioned
957 // calculations by propagating the potential carry and high order bit from
958 // the next word during the right shift.
959 kLen := len(k)
960 var result nafScalar
961 var carry uint8
962 for byteNum := kLen - 1; byteNum >= 0; byteNum-- {
963 // Calculate k/2. Notice the carry from the previous word is added and
964 // the low order bit from the next word is shifted in accordingly.
965 kc := uint16(k[byteNum]) + uint16(carry)
966 var nextWord uint8
967 if byteNum > 0 {
968 nextWord = k[byteNum-1]
969 }
970 halfK := kc>>1 | uint16(nextWord<<7)
971 // Calculate 3k/2 and determine the non-zero digits in the result.
972 threeHalfK := kc + halfK
973 nonZeroResultDigits := threeHalfK ^ halfK
974 // Determine the signed digits {0, ±1}.
975 result.pos[byteNum+1] = uint8(threeHalfK & nonZeroResultDigits)
976 result.neg[byteNum+1] = uint8(halfK & nonZeroResultDigits)
977 // Propagate the potential carry from the 3k/2 calculation.
978 carry = uint8(threeHalfK >> 8)
979 }
980 result.pos[0] = carry
981 // Set the starting and ending positions within the fixed size arrays to
982 // identify the bytes that are actually used. This is important since the
983 // encoding is big endian and thus trailing zero bytes changes its value.
984 result.start = 1 - carry
985 result.end = uint8(kLen + 1)
986 return result
987 }
988 989 // ScalarMultNonConst multiplies k*P where k is a scalar modulo the curve order
990 // and P is a point in Jacobian projective coordinates and stores the result in
991 // the provided Jacobian point.
992 //
993 // NOTE: The point must be normalized for this function to return the correct
994 // result. The resulting point will be normalized.
995 func ScalarMultNonConst(k *ModNScalar, point, result *JacobianPoint) {
996 // -------------------------------------------------------------------------
997 // This makes use of the following efficiently-computable endomorphism to
998 // accelerate the computation:
999 //
1000 // φ(P) ⟼ λ*P = (β*P.x mod p, P.y)
1001 //
1002 // In other words, there is a special scalar λ that every point on the
1003 // elliptic curve can be multiplied by that will result in the same point as
1004 // performing a single field multiplication of the point's X coordinate by
1005 // the special value β.
1006 //
1007 // This is useful because scalar point multiplication is significantly more
1008 // expensive than a single field multiplication given the former involves a
1009 // series of point doublings and additions which themselves consist of a
1010 // combination of several field multiplications, squarings, and additions.
1011 //
1012 // So, the idea behind making use of the endomorphism is thus to decompose
1013 // the scalar into two scalars that are each about half the bit length of
1014 // the original scalar such that:
1015 //
1016 // k ≡ k1 + k2*λ (mod n)
1017 //
1018 // This in turn allows the scalar point multiplication to be performed as a
1019 // sum of two smaller half-length multiplications as follows:
1020 //
1021 // k*P = (k1 + k2*λ)*P
1022 // = k1*P + k2*λ*P
1023 // = k1*P + k2*φ(P)
1024 //
1025 // Thus, a speedup is achieved so long as it's faster to decompose the
1026 // scalar, compute φ(P), and perform a simultaneous multiply of the
1027 // half-length point multiplications than it is to compute a full width
1028 // point multiplication.
1029 //
1030 // In practice, benchmarks show the current implementation provides a
1031 // speedup of around 30-35% versus not using the endomorphism.
1032 //
1033 // See section 3.5 in [GECC] for a more rigorous treatment.
1034 // -------------------------------------------------------------------------
1035 1036 // Per above, the main equation here to remember is:
1037 // k*P = k1*P + k2*φ(P)
1038 //
1039 // p1 below is P in the equation while p2 is φ(P) in the equation.
1040 //
1041 // NOTE: φ(x,y) = (β*x,y). The Jacobian z coordinates are the same, so this
1042 // math goes through.
1043 //
1044 // Also, calculate -p1 and -p2 for use in the NAF optimization.
1045 p1, p1Neg := &JacobianPoint{}, &JacobianPoint{}
1046 p1.Set(point)
1047 p1Neg.Set(p1)
1048 p1Neg.Y.Negate(1).Normalize()
1049 p2, p2Neg := &JacobianPoint{}, &JacobianPoint{}
1050 p2.Set(p1)
1051 eb := _endoBeta()
1052 p2.X.Mul(&eb).Normalize()
1053 p2Neg.Set(p2)
1054 p2Neg.Y.Negate(1).Normalize()
1055 // Decompose k into k1 and k2 such that k = k1 + k2*λ (mod n) where k1 and
1056 // k2 are around half the bit length of k in order to halve the number of EC
1057 // operations.
1058 //
1059 // Notice that this also flips the sign of the scalars and points as needed
1060 // to minimize the bit lengths of the scalars k1 and k2.
1061 //
1062 // This is done because the scalars are operating modulo the group order
1063 // which means that when they would otherwise be a small negative magnitude
1064 // they will instead be a large positive magnitude. Since the goal is for
1065 // the scalars to have a small magnitude to achieve a performance boost, use
1066 // their negation when they are greater than the half order of the group and
1067 // flip the positive and negative values of the corresponding point that
1068 // will be multiplied by to compensate.
1069 //
1070 // In other words, transform the calc when k1 is over the half order to:
1071 // k1*P = -k1*-P
1072 //
1073 // Similarly, transform the calc when k2 is over the half order to:
1074 // k2*φ(P) = -k2*-φ(P)
1075 k1, k2 := splitK(k)
1076 if k1.IsOverHalfOrder() {
1077 k1.Negate()
1078 p1, p1Neg = p1Neg, p1
1079 }
1080 if k2.IsOverHalfOrder() {
1081 k2.Negate()
1082 p2, p2Neg = p2Neg, p2
1083 }
1084 // Convert k1 and k2 into their NAF representations since NAF has a lot more
1085 // zeros overall on average which minimizes the number of required point
1086 // additions in exchange for a mix of fewer point additions and subtractions
1087 // at the cost of one additional point doubling.
1088 //
1089 // This is an excellent tradeoff because subtraction of points has the same
1090 // computational complexity as addition of points and point doubling is
1091 // faster than both.
1092 //
1093 // Concretely, on average, 1/2 of all bits will be non-zero with the normal
1094 // binary representation whereas only 1/3rd of the bits will be non-zero
1095 // with NAF.
1096 //
1097 // The Pos version of the bytes contain the +1s and the Neg versions contain
1098 // the -1s.
1099 k1Bytes, k2Bytes := k1.Bytes(), k2.Bytes()
1100 k1NAF, k2NAF := naf(k1Bytes[:]), naf(k2Bytes[:])
1101 k1PosNAF, k1NegNAF := k1NAF.Pos(), k1NAF.Neg()
1102 k2PosNAF, k2NegNAF := k2NAF.Pos(), k2NAF.Neg()
1103 k1Len, k2Len := len(k1PosNAF), len(k2PosNAF)
1104 // Add left-to-right using the NAF optimization. See algorithm 3.77 from
1105 // [GECC].
1106 //
1107 // Point Q = ∞ (point at infinity).
1108 var q JacobianPoint
1109 m := k1Len
1110 if m < k2Len {
1111 m = k2Len
1112 }
1113 for i := 0; i < m; i++ {
1114 // Since k1 and k2 are potentially different lengths and the calculation
1115 // is being done left to right, pad the front of the shorter one with
1116 // 0s.
1117 var k1BytePos, k1ByteNeg, k2BytePos, k2ByteNeg byte
1118 if i >= m-k1Len {
1119 k1BytePos, k1ByteNeg = k1PosNAF[i-(m-k1Len)], k1NegNAF[i-(m-k1Len)]
1120 }
1121 if i >= m-k2Len {
1122 k2BytePos, k2ByteNeg = k2PosNAF[i-(m-k2Len)], k2NegNAF[i-(m-k2Len)]
1123 }
1124 for mask := uint8(1 << 7); mask > 0; mask >>= 1 {
1125 // Q = 2 * Q
1126 DoubleNonConst(&q, &q)
1127 // Add or subtract the first point based on the signed digit of the
1128 // NAF representation of k1 at this bit position.
1129 //
1130 // +1: Q = Q + p1
1131 // -1: Q = Q - p1
1132 // 0: Q = Q (no change)
1133 if k1BytePos&mask == mask {
1134 AddNonConst(&q, p1, &q)
1135 } else if k1ByteNeg&mask == mask {
1136 AddNonConst(&q, p1Neg, &q)
1137 }
1138 // Add or subtract the second point based on the signed digit of the
1139 // NAF representation of k2 at this bit position.
1140 //
1141 // +1: Q = Q + p2
1142 // -1: Q = Q - p2
1143 // 0: Q = Q (no change)
1144 if k2BytePos&mask == mask {
1145 AddNonConst(&q, p2, &q)
1146 } else if k2ByteNeg&mask == mask {
1147 AddNonConst(&q, p2Neg, &q)
1148 }
1149 }
1150 }
1151 result.Set(&q)
1152 }
1153 1154 // ScalarBaseMultNonConst multiplies k*G where k is a scalar modulo the curve
1155 // order and G is the base point of the group and stores the result in the
1156 // provided Jacobian point.
1157 //
1158 // NOTE: The resulting point will be normalized.
1159 func ScalarBaseMultNonConst(k *ModNScalar, result *JacobianPoint) {
1160 bytePoints := getBytePointTable()
1161 // Start with the point at infinity.
1162 result.X.Zero()
1163 result.Y.Zero()
1164 result.Z.Zero()
1165 // bytePoints has all 256 byte points for each 8-bit window. The strategy
1166 // is to add up the byte points. This is best understood by expressing k in
1167 // base-256 which it already sort of is. Each "digit" in the 8-bit window
1168 // can be looked up using bytePoints and added together.
1169 kb := k.Bytes()
1170 for i := 0; i < len(kb); i++ {
1171 pt := &bytePoints[i][kb[i]]
1172 AddNonConst(result, pt, result)
1173 }
1174 }
1175 1176 // isOnCurve returns whether or not the affine point (x,y) is on the curve.
1177 func isOnCurve(fx, fy *FieldVal) bool {
1178 // Elliptic curve equation for secp256k1 is: y^2 = x^3 + 7
1179 y2 := (&FieldVal{}).SquareVal(fy).Normalize()
1180 result := (&FieldVal{}).SquareVal(fx).Mul(fx).AddInt(7).Normalize()
1181 return y2.Equals(result)
1182 }
1183 1184 // DecompressY attempts to calculate the Y coordinate for the given X coordinate
1185 // such that the result pair is a point on the secp256k1 curve. It adjusts Y
1186 // based on the desired oddness and returns whether or not it was successful
1187 // since not all X coordinates are valid.
1188 //
1189 // The magnitude of the provided X coordinate field val must be a max of 8 for a
1190 // correct result. The resulting Y field val will have a max magnitude of 2.
1191 func DecompressY(x *FieldVal, odd bool, resultY *FieldVal) bool {
1192 // The curve equation for secp256k1 is: y^2 = x^3 + 7. Thus
1193 // y = +-sqrt(x^3 + 7).
1194 //
1195 // The x coordinate must be invalid if there is no square root for the
1196 // calculated rhs because it means the X coordinate is not for a point on
1197 // the curve.
1198 x3PlusB := (&FieldVal{}).SquareVal(x).Mul(x).AddInt(7)
1199 if hasSqrt := resultY.SquareRootVal(x3PlusB); !hasSqrt {
1200 return false
1201 }
1202 if resultY.Normalize().IsOdd() != odd {
1203 resultY.Negate(1)
1204 }
1205 return true
1206 }
1207