Inside Bitcoin

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
7It is designed for ease of understanding, not performance.
8
9WARNING: This code is slow and trivially vulnerable to side channel attacks. Do not use for
10anything but tests.
11
12Exports:
13* FE: class for secp256k1 field elements
14* GE: class for secp256k1 group elements
15* G: the secp256k1 generator point
16"""
17
18import unittest
19from hashlib import sha256
20from test_framework.util import assert_equal, assert_not_equal
21
22class 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
154class 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()