1 # Copyright (c) 2022-present The Bitcoin Core developers
2 # Distributed under the MIT software license, see the accompanying
3 # file COPYING or http://www.opensource.org/licenses/mit-license.php.
4 5 """Test-only implementation of low-level secp256k1 field and group arithmetic
6 7 It is designed for ease of understanding, not performance.
8 9 WARNING: This code is slow and trivially vulnerable to side channel attacks. Do not use for
10 anything but tests.
11 12 Exports:
13 * FE: class for secp256k1 field elements
14 * GE: class for secp256k1 group elements
15 * G: the secp256k1 generator point
16 """
17 18 import unittest
19 from hashlib import sha256
20 from test_framework.util import assert_equal, assert_not_equal
21 22 class FE:
23 """Objects of this class represent elements of the field GF(2**256 - 2**32 - 977).
24 25 They are represented internally in numerator / denominator form, in order to delay inversions.
26 """
27 28 # The size of the field (also its modulus and characteristic).
29 SIZE = 2**256 - 2**32 - 977
30 31 def __init__(self, a=0, b=1):
32 """Initialize a field element a/b; both a and b can be ints or field elements."""
33 if isinstance(a, FE):
34 num = a._num
35 den = a._den
36 else:
37 num = a % FE.SIZE
38 den = 1
39 if isinstance(b, FE):
40 den = (den * b._num) % FE.SIZE
41 num = (num * b._den) % FE.SIZE
42 else:
43 den = (den * b) % FE.SIZE
44 assert_not_equal(den, 0)
45 if num == 0:
46 den = 1
47 self._num = num
48 self._den = den
49 50 def __add__(self, a):
51 """Compute the sum of two field elements (second may be int)."""
52 if isinstance(a, FE):
53 return FE(self._num * a._den + self._den * a._num, self._den * a._den)
54 return FE(self._num + self._den * a, self._den)
55 56 def __radd__(self, a):
57 """Compute the sum of an integer and a field element."""
58 return FE(a) + self
59 60 def __sub__(self, a):
61 """Compute the difference of two field elements (second may be int)."""
62 if isinstance(a, FE):
63 return FE(self._num * a._den - self._den * a._num, self._den * a._den)
64 return FE(self._num - self._den * a, self._den)
65 66 def __rsub__(self, a):
67 """Compute the difference of an integer and a field element."""
68 return FE(a) - self
69 70 def __mul__(self, a):
71 """Compute the product of two field elements (second may be int)."""
72 if isinstance(a, FE):
73 return FE(self._num * a._num, self._den * a._den)
74 return FE(self._num * a, self._den)
75 76 def __rmul__(self, a):
77 """Compute the product of an integer with a field element."""
78 return FE(a) * self
79 80 def __truediv__(self, a):
81 """Compute the ratio of two field elements (second may be int)."""
82 return FE(self, a)
83 84 def __pow__(self, a):
85 """Raise a field element to an integer power."""
86 return FE(pow(self._num, a, FE.SIZE), pow(self._den, a, FE.SIZE))
87 88 def __neg__(self):
89 """Negate a field element."""
90 return FE(-self._num, self._den)
91 92 def __int__(self):
93 """Convert a field element to an integer in range 0..p-1. The result is cached."""
94 if self._den != 1:
95 self._num = (self._num * pow(self._den, -1, FE.SIZE)) % FE.SIZE
96 self._den = 1
97 return self._num
98 99 def sqrt(self):
100 """Compute the square root of a field element if it exists (None otherwise).
101 102 Due to the fact that our modulus is of the form (p % 4) == 3, the Tonelli-Shanks
103 algorithm (https://en.wikipedia.org/wiki/Tonelli-Shanks_algorithm) is simply
104 raising the argument to the power (p + 1) / 4.
105 106 To see why: (p-1) % 2 = 0, so 2 divides the order of the multiplicative group,
107 and thus only half of the non-zero field elements are squares. An element a is
108 a (nonzero) square when Euler's criterion, a^((p-1)/2) = 1 (mod p), holds. We're
109 looking for x such that x^2 = a (mod p). Given a^((p-1)/2) = 1, that is equivalent
110 to x^2 = a^(1 + (p-1)/2) mod p. As (1 + (p-1)/2) is even, this is equivalent to
111 x = a^((1 + (p-1)/2)/2) mod p, or x = a^((p+1)/4) mod p."""
112 v = int(self)
113 s = pow(v, (FE.SIZE + 1) // 4, FE.SIZE)
114 if s**2 % FE.SIZE == v:
115 return FE(s)
116 return None
117 118 def is_square(self):
119 """Determine if this field element has a square root."""
120 # A more efficient algorithm is possible here (Jacobi symbol).
121 return self.sqrt() is not None
122 123 def is_even(self):
124 """Determine whether this field element, represented as integer in 0..p-1, is even."""
125 return int(self) & 1 == 0
126 127 def __eq__(self, a):
128 """Check whether two field elements are equal (second may be an int)."""
129 if isinstance(a, FE):
130 return (self._num * a._den - self._den * a._num) % FE.SIZE == 0
131 return (self._num - self._den * a) % FE.SIZE == 0
132 133 def to_bytes(self):
134 """Convert a field element to a 32-byte array (BE byte order)."""
135 return int(self).to_bytes(32, 'big')
136 137 @staticmethod
138 def from_bytes(b):
139 """Convert a 32-byte array to a field element (BE byte order, no overflow allowed)."""
140 v = int.from_bytes(b, 'big')
141 if v >= FE.SIZE:
142 return None
143 return FE(v)
144 145 def __str__(self):
146 """Convert this field element to a 64 character hex string."""
147 return f"{int(self):064x}"
148 149 def __repr__(self):
150 """Get a string representation of this field element."""
151 return f"FE(0x{int(self):x})"
152 153 154 class GE:
155 """Objects of this class represent secp256k1 group elements (curve points or infinity)
156 157 Normal points on the curve have fields:
158 * x: the x coordinate (a field element)
159 * y: the y coordinate (a field element, satisfying y^2 = x^3 + 7)
160 * infinity: False
161 162 The point at infinity has field:
163 * infinity: True
164 """
165 166 # Order of the group (number of points on the curve, plus 1 for infinity)
167 ORDER = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141
168 169 # Number of valid distinct x coordinates on the curve.
170 ORDER_HALF = ORDER // 2
171 172 def __init__(self, x=None, y=None):
173 """Initialize a group element with specified x and y coordinates, or infinity."""
174 if x is None:
175 # Initialize as infinity.
176 assert y is None
177 self.infinity = True
178 else:
179 # Initialize as point on the curve (and check that it is).
180 fx = FE(x)
181 fy = FE(y)
182 assert_equal(fy**2, fx**3 + 7)
183 self.infinity = False
184 self.x = fx
185 self.y = fy
186 187 def __add__(self, a):
188 """Add two group elements together."""
189 # Deal with infinity: a + infinity == infinity + a == a.
190 if self.infinity:
191 return a
192 if a.infinity:
193 return self
194 if self.x == a.x:
195 if self.y != a.y:
196 # A point added to its own negation is infinity.
197 assert_equal(self.y + a.y, 0)
198 return GE()
199 else:
200 # For identical inputs, use the tangent (doubling formula).
201 lam = (3 * self.x**2) / (2 * self.y)
202 else:
203 # For distinct inputs, use the line through both points (adding formula).
204 lam = (self.y - a.y) / (self.x - a.x)
205 # Determine point opposite to the intersection of that line with the curve.
206 x = lam**2 - (self.x + a.x)
207 y = lam * (self.x - x) - self.y
208 return GE(x, y)
209 210 @staticmethod
211 def mul(*aps):
212 """Compute a (batch) scalar group element multiplication.
213 214 GE.mul((a1, p1), (a2, p2), (a3, p3)) is identical to a1*p1 + a2*p2 + a3*p3,
215 but more efficient."""
216 # Reduce all the scalars modulo order first (so we can deal with negatives etc).
217 naps = [(a % GE.ORDER, p) for a, p in aps]
218 # Start with point at infinity.
219 r = GE()
220 # Iterate over all bit positions, from high to low.
221 for i in range(255, -1, -1):
222 # Double what we have so far.
223 r = r + r
224 # Add then add the points for which the corresponding scalar bit is set.
225 for (a, p) in naps:
226 if (a >> i) & 1:
227 r += p
228 return r
229 230 def __rmul__(self, a):
231 """Multiply an integer with a group element."""
232 if self == G:
233 return FAST_G.mul(a)
234 return GE.mul((a, self))
235 236 def __neg__(self):
237 """Compute the negation of a group element."""
238 if self.infinity:
239 return self
240 return GE(self.x, -self.y)
241 242 def to_bytes_compressed(self):
243 """Convert a non-infinite group element to 33-byte compressed encoding."""
244 assert not self.infinity
245 return bytes([3 - self.y.is_even()]) + self.x.to_bytes()
246 247 def to_bytes_uncompressed(self):
248 """Convert a non-infinite group element to 65-byte uncompressed encoding."""
249 assert not self.infinity
250 return b'\x04' + self.x.to_bytes() + self.y.to_bytes()
251 252 def to_bytes_xonly(self):
253 """Convert (the x coordinate of) a non-infinite group element to 32-byte xonly encoding."""
254 assert not self.infinity
255 return self.x.to_bytes()
256 257 @staticmethod
258 def lift_x(x):
259 """Return group element with specified field element as x coordinate (and even y)."""
260 y = (FE(x)**3 + 7).sqrt()
261 if y is None:
262 return None
263 if not y.is_even():
264 y = -y
265 return GE(x, y)
266 267 @staticmethod
268 def from_bytes(b):
269 """Convert a compressed or uncompressed encoding to a group element."""
270 assert len(b) in (33, 65)
271 if len(b) == 33:
272 if b[0] != 2 and b[0] != 3:
273 return None
274 x = FE.from_bytes(b[1:])
275 if x is None:
276 return None
277 r = GE.lift_x(x)
278 if r is None:
279 return None
280 if b[0] == 3:
281 r = -r
282 return r
283 else:
284 if b[0] != 4:
285 return None
286 x = FE.from_bytes(b[1:33])
287 y = FE.from_bytes(b[33:])
288 if y**2 != x**3 + 7:
289 return None
290 return GE(x, y)
291 292 @staticmethod
293 def from_bytes_xonly(b):
294 """Convert a point given in xonly encoding to a group element."""
295 assert_equal(len(b), 32)
296 x = FE.from_bytes(b)
297 if x is None:
298 return None
299 return GE.lift_x(x)
300 301 @staticmethod
302 def is_valid_x(x):
303 """Determine whether the provided field element is a valid X coordinate."""
304 return (FE(x)**3 + 7).is_square()
305 306 def __str__(self):
307 """Convert this group element to a string."""
308 if self.infinity:
309 return "(inf)"
310 return f"({self.x},{self.y})"
311 312 def __repr__(self):
313 """Get a string representation for this group element."""
314 if self.infinity:
315 return "GE()"
316 return f"GE(0x{int(self.x):x},0x{int(self.y):x})"
317 318 # The secp256k1 generator point
319 G = GE.lift_x(0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798)
320 321 322 class FastGEMul:
323 """Table for fast multiplication with a constant group element.
324 325 Speed up scalar multiplication with a fixed point P by using a precomputed lookup table with
326 its powers of 2:
327 328 table = [P, 2*P, 4*P, (2^3)*P, (2^4)*P, ..., (2^255)*P]
329 330 During multiplication, the points corresponding to each bit set in the scalar are added up,
331 i.e. on average ~128 point additions take place.
332 """
333 334 def __init__(self, p):
335 self.table = [p] # table[i] = (2^i) * p
336 for _ in range(255):
337 p = p + p
338 self.table.append(p)
339 340 def mul(self, a):
341 result = GE()
342 a = a % GE.ORDER
343 for bit in range(a.bit_length()):
344 if a & (1 << bit):
345 result += self.table[bit]
346 return result
347 348 # Precomputed table with multiples of G for fast multiplication
349 FAST_G = FastGEMul(G)
350 351 class TestFrameworkSecp256k1(unittest.TestCase):
352 def test_H(self):
353 H = sha256(G.to_bytes_uncompressed()).digest()
354 assert GE.lift_x(FE.from_bytes(H)) is not None
355 self.assertEqual(H.hex(), "50929b74c1a04954b78b4b6035e97a5e078a5a0f28ec96d547bfee9ace803ac0")
356