|
| 1 | +#!/usr/bin/env python3 |
| 2 | + |
| 3 | +import argparse |
| 4 | + |
| 5 | +def modinv(a, n): |
| 6 | + """Compute the modular inverse of a modulo n using the extended Euclidean |
| 7 | + Algorithm. See https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Modular_integers. |
| 8 | + """ |
| 9 | + # TODO: Change to pow(a, -1, n) available in Python 3.8 |
| 10 | + t1, t2 = 0, 1 |
| 11 | + r1, r2 = n, a |
| 12 | + while r2 != 0: |
| 13 | + q = r1 // r2 |
| 14 | + t1, t2 = t2, t1 - q * t2 |
| 15 | + r1, r2 = r2, r1 - q * r2 |
| 16 | + if r1 > 1: |
| 17 | + return None |
| 18 | + if t1 < 0: |
| 19 | + t1 += n |
| 20 | + return t1 |
| 21 | + |
| 22 | +def modsqrt(a, p): |
| 23 | + """Compute the square root of a modulo p when p % 4 = 3. |
| 24 | + The Tonelli-Shanks algorithm can be used. See https://en.wikipedia.org/wiki/Tonelli-Shanks_algorithm |
| 25 | + Limiting this function to only work for p % 4 = 3 means we don't need to |
| 26 | + iterate through the loop. The highest n such that p - 1 = 2^n Q with Q odd |
| 27 | + is n = 1. Therefore Q = (p-1)/2 and sqrt = a^((Q+1)/2) = a^((p+1)/4) |
| 28 | + secp256k1's is defined over field of size 2**256 - 2**32 - 977, which is 3 mod 4. |
| 29 | + """ |
| 30 | + if p % 4 != 3: |
| 31 | + raise NotImplementedError("modsqrt only implemented for p % 4 = 3") |
| 32 | + sqrt = pow(a, (p + 1)//4, p) |
| 33 | + if pow(sqrt, 2, p) == a % p: |
| 34 | + return sqrt |
| 35 | + return None |
| 36 | + |
| 37 | +class EllipticCurve: |
| 38 | + def __init__(self, p, a, b): |
| 39 | + """Initialize elliptic curve y^2 = x^3 + a*x + b over GF(p).""" |
| 40 | + self.p = p |
| 41 | + self.a = a % p |
| 42 | + self.b = b % p |
| 43 | + |
| 44 | + def affine(self, p1): |
| 45 | + """Convert a Jacobian point tuple p1 to affine form, or None if at infinity. |
| 46 | + An affine point is represented as the Jacobian (x, y, 1)""" |
| 47 | + x1, y1, z1 = p1 |
| 48 | + if z1 == 0: |
| 49 | + return None |
| 50 | + inv = modinv(z1, self.p) |
| 51 | + inv_2 = (inv**2) % self.p |
| 52 | + inv_3 = (inv_2 * inv) % self.p |
| 53 | + return ((inv_2 * x1) % self.p, (inv_3 * y1) % self.p, 1) |
| 54 | + |
| 55 | + def lift_x(self, x): |
| 56 | + """Given an X coordinate on the curve, return a corresponding affine point for which the Y coordinate is even.""" |
| 57 | + x_3 = pow(x, 3, self.p) |
| 58 | + v = x_3 + self.a * x + self.b |
| 59 | + y = modsqrt(v, self.p) |
| 60 | + if y is None: |
| 61 | + return None |
| 62 | + return (x, self.p - y if y & 1 else y, 1) |
| 63 | + |
| 64 | + def double(self, p1): |
| 65 | + """Double a Jacobian tuple p1 |
| 66 | + See https://en.wikibooks.org/wiki/Cryptography/Prime_Curve/Jacobian_Coordinates - Point Doubling""" |
| 67 | + x1, y1, z1 = p1 |
| 68 | + if z1 == 0: |
| 69 | + return (0, 1, 0) |
| 70 | + y1_2 = (y1**2) % self.p |
| 71 | + y1_4 = (y1_2**2) % self.p |
| 72 | + x1_2 = (x1**2) % self.p |
| 73 | + s = (4*x1*y1_2) % self.p |
| 74 | + m = 3*x1_2 |
| 75 | + if self.a: |
| 76 | + m += self.a * pow(z1, 4, self.p) |
| 77 | + m = m % self.p |
| 78 | + x2 = (m**2 - 2*s) % self.p |
| 79 | + y2 = (m*(s - x2) - 8*y1_4) % self.p |
| 80 | + z2 = (2*y1*z1) % self.p |
| 81 | + return (x2, y2, z2) |
| 82 | + |
| 83 | + def add(self, p1, p2): |
| 84 | + """Add two Jacobian tuples p1 and p2 |
| 85 | + See https://en.wikibooks.org/wiki/Cryptography/Prime_Curve/Jacobian_Coordinates - Point Addition""" |
| 86 | + x1, y1, z1 = p1 |
| 87 | + x2, y2, z2 = p2 |
| 88 | + # Adding the point at infinity is a no-op |
| 89 | + if z1 == 0: |
| 90 | + return p2 |
| 91 | + if z2 == 0: |
| 92 | + return p1 |
| 93 | + z1_2 = (z1**2) % self.p |
| 94 | + z1_3 = (z1_2 * z1) % self.p |
| 95 | + z2_2 = (z2**2) % self.p |
| 96 | + z2_3 = (z2_2 * z2) % self.p |
| 97 | + u1 = (x1 * z2_2) % self.p |
| 98 | + u2 = (x2 * z1_2) % self.p |
| 99 | + s1 = (y1 * z2_3) % self.p |
| 100 | + s2 = (y2 * z1_3) % self.p |
| 101 | + if u1 == u2: |
| 102 | + if (s1 != s2): |
| 103 | + # p1 and p2 are inverses. Return the point at infinity. |
| 104 | + return (0, 1, 0) |
| 105 | + # p1 == p2. The formulas below fail when the two points are equal. |
| 106 | + return self.double(p1) |
| 107 | + h = u2 - u1 |
| 108 | + r = s2 - s1 |
| 109 | + h_2 = (h**2) % self.p |
| 110 | + h_3 = (h_2 * h) % self.p |
| 111 | + u1_h_2 = (u1 * h_2) % self.p |
| 112 | + x3 = (r**2 - h_3 - 2*u1_h_2) % self.p |
| 113 | + y3 = (r*(u1_h_2 - x3) - s1*h_3) % self.p |
| 114 | + z3 = (h*z1*z2) % self.p |
| 115 | + return (x3, y3, z3) |
| 116 | + |
| 117 | + def mul(self, p, n): |
| 118 | + """Compute a point multiplication of Jacobian point p times n.""" |
| 119 | + r = (0, 1, 0) |
| 120 | + for i in range(255, -1, -1): |
| 121 | + r = self.double(r) |
| 122 | + if ((n >> i) & 1): |
| 123 | + r = self.add(r, p) |
| 124 | + return r |
| 125 | + |
| 126 | +# The secp256k1 field size |
| 127 | +SECP256K1_FIELD_SIZE = 2**256 - 2**32 - 977 |
| 128 | +# The secp256k1 curve itself |
| 129 | +SECP256K1 = EllipticCurve(SECP256K1_FIELD_SIZE, 0, 7) |
| 130 | +# The order of the secp256k1 curve |
| 131 | +SECP256K1_ORDER = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEBAAEDCE6AF48A03BBFD25E8CD0364141 |
| 132 | +# The standard generator |
| 133 | +SECP256K1_G = (0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798, 0x483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8, 1) |
| 134 | +# An alternative nothing-up-my-sleeve generator (with unknown DL w.r.t. G; only used for blinding) |
| 135 | +SECP256K1_U = SECP256K1.add(SECP256K1.lift_x(int.from_bytes(b"The scalar for this x is unknown", 'big')), SECP256K1_G) |
| 136 | + |
| 137 | +# Parse command-line arguments |
| 138 | +parser = argparse.ArgumentParser(description="Generate the precomputed context for libsecp256k1.") |
| 139 | +parser.add_argument('--ecmult-gen-precision', '-g', type=int, choices=[2,4,8], default=4, help="Precision bits to tune the precomputed table size for signing. Valid options: 2, 4, 8. The default is 4.") |
| 140 | +args = parser.parse_args() |
| 141 | + |
| 142 | +# Derive constants like ecmult_gen.h does |
| 143 | +ECMULT_GEN_PREC_B = args.ecmult_gen_precision |
| 144 | +ECMULT_GEN_PREC_G = 1 << ECMULT_GEN_PREC_B |
| 145 | +ECMULT_GEN_PREC_N = 256 // ECMULT_GEN_PREC_B |
| 146 | + |
| 147 | +# Compute precomputed points and output |
| 148 | +print("#ifndef SECP256K1_ECMULT_STATIC_CONTEXT_H") |
| 149 | +print("#define SECP256K1_ECMULT_STATIC_CONTEXT_H") |
| 150 | +print("#include \"src/group.h\"") |
| 151 | +print("#define SC SECP256K1_GE_STORAGE_CONST") |
| 152 | +print("#if ECMULT_GEN_PREC_N != %d || ECMULT_GEN_PREC_G != %d" % (ECMULT_GEN_PREC_N, ECMULT_GEN_PREC_G)) |
| 153 | +print(" #error configuration mismatch, invalid ECMULT_GEN_PREC_N, ECMULT_GEN_PREC_G. Try deleting ecmult_static_context.h before the build.") |
| 154 | +print("#endif"); |
| 155 | +print("static const secp256k1_ge_storage secp256k1_ecmult_static_context[ECMULT_GEN_PREC_N][ECMULT_GEN_PREC_G] = {") |
| 156 | +for outer in range(ECMULT_GEN_PREC_N): |
| 157 | + print("{") |
| 158 | + # All but the last bucket use SECP256K1_U * 2^outer as blinding term. The last one uses the negation of the |
| 159 | + # sum of all previous ones (so that they cancel out to 0). |
| 160 | + numsbase = SECP256K1.mul(SECP256K1_U, (1 << outer) if outer + 1 != ECMULT_GEN_PREC_N else (1 - (1 << outer)) % SECP256K1_ORDER) |
| 161 | + for inner in range(ECMULT_GEN_PREC_G): |
| 162 | + point = SECP256K1.affine(SECP256K1.add(SECP256K1.mul(SECP256K1_G, inner << (ECMULT_GEN_PREC_B * outer)), numsbase)) |
| 163 | + print(" SC(%uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu, %uu)%s" % tuple( |
| 164 | + [point[v] >> (32*(7-i)) & 0xFFFFFFFF for v in range(2) for i in range(8)] + |
| 165 | + ["," if inner + 1 != ECMULT_GEN_PREC_G else ""])) |
| 166 | + print("}," if outer + 1 != ECMULT_GEN_PREC_N else "}") |
| 167 | +print("};") |
| 168 | +print("#undef SC") |
| 169 | +print("#endif") |
0 commit comments