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