diff --git a/src/secp256k1/sage/gen_exhaustive_groups.sage b/src/secp256k1/sage/gen_exhaustive_groups.sage index 3c3c98481..01d15dcde 100644 --- a/src/secp256k1/sage/gen_exhaustive_groups.sage +++ b/src/secp256k1/sage/gen_exhaustive_groups.sage @@ -1,129 +1,124 @@ -# Define field size and field -P = 2^256 - 2^32 - 977 -F = GF(P) -BETA = F(0x7ae96a2b657c07106e64479eac3434e99cf0497512f58995c1396c28719501ee) - -assert(BETA != F(1) and BETA^3 == F(1)) +load("secp256k1_params.sage") orders_done = set() results = {} first = True for b in range(1, P): # There are only 6 curves (up to isomorphism) of the form y^2=x^3+B. Stop once we have tried all. if len(orders_done) == 6: break E = EllipticCurve(F, [0, b]) print("Analyzing curve y^2 = x^3 + %i" % b) n = E.order() # Skip curves with an order we've already tried if n in orders_done: print("- Isomorphic to earlier curve") continue orders_done.add(n) # Skip curves isomorphic to the real secp256k1 if n.is_pseudoprime(): print(" - Isomorphic to secp256k1") continue print("- Finding subgroups") # Find what prime subgroups exist for f, _ in n.factor(): print("- Analyzing subgroup of order %i" % f) # Skip subgroups of order >1000 if f < 4 or f > 1000: print(" - Bad size") continue # Iterate over X coordinates until we find one that is on the curve, has order f, # and for which curve isomorphism exists that maps it to X coordinate 1. for x in range(1, P): # Skip X coordinates not on the curve, and construct the full point otherwise. if not E.is_x_coord(x): continue G = E.lift_x(F(x)) print(" - Analyzing (multiples of) point with X=%i" % x) # Skip points whose order is not a multiple of f. Project the point to have # order f otherwise. if (G.order() % f): print(" - Bad order") continue G = G * (G.order() // f) # Find lambda for endomorphism. Skip if none can be found. lam = None for l in Integers(f)(1).nth_root(3, all=True): if int(l)*G == E(BETA*G[0], G[1]): lam = int(l) break if lam is None: print(" - No endomorphism for this subgroup") break # Now look for an isomorphism of the curve that gives this point an X # coordinate equal to 1. # If (x,y) is on y^2 = x^3 + b, then (a^2*x, a^3*y) is on y^2 = x^3 + a^6*b. # So look for m=a^2=1/x. m = F(1)/G[0] if not m.is_square(): print(" - No curve isomorphism maps it to a point with X=1") continue a = m.sqrt() rb = a^6*b RE = EllipticCurve(F, [0, rb]) # Use as generator twice the image of G under the above isormorphism. # This means that generator*(1/2 mod f) will have X coordinate 1. RG = RE(1, a^3*G[1]) * 2 # And even Y coordinate. if int(RG[1]) % 2: RG = -RG assert(RG.order() == f) assert(lam*RG == RE(BETA*RG[0], RG[1])) # We have found curve RE:y^2=x^3+rb with generator RG of order f. Remember it results[f] = {"b": rb, "G": RG, "lambda": lam} print(" - Found solution") break print("") print("") print("") print("/* To be put in src/group_impl.h: */") first = True for f in sorted(results.keys()): b = results[f]["b"] G = results[f]["G"] print("# %s EXHAUSTIVE_TEST_ORDER == %i" % ("if" if first else "elif", f)) first = False print("static const secp256k1_ge secp256k1_ge_const_g = SECP256K1_GE_CONST(") print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[0]) >> (32 * (7 - i))) & 0xffffffff for i in range(4))) print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[0]) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8))) print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(G[1]) >> (32 * (7 - i))) & 0xffffffff for i in range(4))) print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x" % tuple((int(G[1]) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8))) print(");") print("static const secp256k1_fe secp256k1_fe_const_b = SECP256K1_FE_CONST(") print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x," % tuple((int(b) >> (32 * (7 - i))) & 0xffffffff for i in range(4))) print(" 0x%08x, 0x%08x, 0x%08x, 0x%08x" % tuple((int(b) >> (32 * (7 - i))) & 0xffffffff for i in range(4, 8))) print(");") print("# else") print("# error No known generator for the specified exhaustive test group order.") print("# endif") print("") print("") print("/* To be put in src/scalar_impl.h: */") first = True for f in sorted(results.keys()): lam = results[f]["lambda"] print("# %s EXHAUSTIVE_TEST_ORDER == %i" % ("if" if first else "elif", f)) first = False print("# define EXHAUSTIVE_TEST_LAMBDA %i" % lam) print("# else") print("# error No known lambda for the specified exhaustive test group order.") print("# endif") print("") diff --git a/src/secp256k1/sage/gen_split_lambda_constants.sage b/src/secp256k1/sage/gen_split_lambda_constants.sage new file mode 100644 index 000000000..7d4359e0f --- /dev/null +++ b/src/secp256k1/sage/gen_split_lambda_constants.sage @@ -0,0 +1,114 @@ +""" Generates the constants used in secp256k1_scalar_split_lambda. + +See the comments for secp256k1_scalar_split_lambda in src/scalar_impl.h for detailed explanations. +""" + +load("secp256k1_params.sage") + +def inf_norm(v): + """Returns the infinity norm of a vector.""" + return max(map(abs, v)) + +def gauss_reduction(i1, i2): + v1, v2 = i1.copy(), i2.copy() + while True: + if inf_norm(v2) < inf_norm(v1): + v1, v2 = v2, v1 + # This is essentially + # m = round((v1[0]*v2[0] + v1[1]*v2[1]) / (inf_norm(v1)**2)) + # (rounding to the nearest integer) without relying on floating point arithmetic. + m = ((v1[0]*v2[0] + v1[1]*v2[1]) + (inf_norm(v1)**2) // 2) // (inf_norm(v1)**2) + if m == 0: + return v1, v2 + v2[0] -= m*v1[0] + v2[1] -= m*v1[1] + +def find_split_constants_gauss(): + """Find constants for secp256k1_scalar_split_lamdba using gauss reduction.""" + (v11, v12), (v21, v22) = gauss_reduction([0, N], [1, int(LAMBDA)]) + + # We use related vectors in secp256k1_scalar_split_lambda. + A1, B1 = -v21, -v11 + A2, B2 = v22, -v21 + + return A1, B1, A2, B2 + +def find_split_constants_explicit_tof(): + """Find constants for secp256k1_scalar_split_lamdba using the trace of Frobenius. + + See Benjamin Smith: "Easy scalar decompositions for efficient scalar multiplication on + elliptic curves and genus 2 Jacobians" (https://eprint.iacr.org/2013/672), Example 2 + """ + assert P % 3 == 1 # The paper says P % 3 == 2 but that appears to be a mistake, see [10]. + assert C.j_invariant() == 0 + + t = C.trace_of_frobenius() + + c = Integer(sqrt((4*P - t**2)/3)) + A1 = Integer((t - c)/2 - 1) + B1 = c + + A2 = Integer((t + c)/2 - 1) + B2 = Integer(1 - (t - c)/2) + + # We use a negated b values in secp256k1_scalar_split_lambda. + B1, B2 = -B1, -B2 + + return A1, B1, A2, B2 + +A1, B1, A2, B2 = find_split_constants_explicit_tof() + +# For extra fun, use an independent method to recompute the constants. +assert (A1, B1, A2, B2) == find_split_constants_gauss() + +# PHI : Z[l] -> Z_n where phi(a + b*l) == a + b*lambda mod n. +def PHI(a,b): + return Z(a + LAMBDA*b) + +# Check that (A1, B1) and (A2, B2) are in the kernel of PHI. +assert PHI(A1, B1) == Z(0) +assert PHI(A2, B2) == Z(0) + +# Check that the parallelogram generated by (A1, A2) and (B1, B2) +# is a fundamental domain by containing exactly N points. +# Since the LHS is the determinant and N != 0, this also checks that +# (A1, A2) and (B1, B2) are linearly independent. By the previous +# assertions, (A1, A2) and (B1, B2) are a basis of the kernel. +assert A1*B2 - B1*A2 == N + +# Check that their components are short enough. +assert (A1 + A2)/2 < sqrt(N) +assert B1 < sqrt(N) +assert B2 < sqrt(N) + +G1 = round((2**384)*B2/N) +G2 = round((2**384)*(-B1)/N) + +def rnddiv2(v): + if v & 1: + v += 1 + return v >> 1 + +def scalar_lambda_split(k): + """Equivalent to secp256k1_scalar_lambda_split().""" + c1 = rnddiv2((k * G1) >> 383) + c2 = rnddiv2((k * G2) >> 383) + c1 = (c1 * -B1) % N + c2 = (c2 * -B2) % N + r2 = (c1 + c2) % N + r1 = (k + r2 * -LAMBDA) % N + return (r1, r2) + +# The result of scalar_lambda_split can depend on the representation of k (mod n). +SPECIAL = (2**383) // G2 + 1 +assert scalar_lambda_split(SPECIAL) != scalar_lambda_split(SPECIAL + N) + +print(' A1 =', hex(A1)) +print(' -B1 =', hex(-B1)) +print(' A2 =', hex(A2)) +print(' -B2 =', hex(-B2)) +print(' =', hex(Z(-B2))) +print(' -LAMBDA =', hex(-LAMBDA)) + +print(' G1 =', hex(G1)) +print(' G2 =', hex(G2)) diff --git a/src/secp256k1/sage/secp256k1.sage b/src/secp256k1/sage/prove_group_implementations.sage similarity index 100% rename from src/secp256k1/sage/secp256k1.sage rename to src/secp256k1/sage/prove_group_implementations.sage diff --git a/src/secp256k1/sage/secp256k1_params.sage b/src/secp256k1/sage/secp256k1_params.sage new file mode 100644 index 000000000..4e000726e --- /dev/null +++ b/src/secp256k1/sage/secp256k1_params.sage @@ -0,0 +1,36 @@ +"""Prime order of finite field underlying secp256k1 (2^256 - 2^32 - 977)""" +P = 0xFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F + +"""Finite field underlying secp256k1""" +F = FiniteField(P) + +"""Elliptic curve secp256k1: y^2 = x^3 + 7""" +C = EllipticCurve([F(0), F(7)]) + +"""Base point of secp256k1""" +G = C.lift_x(0x79BE667EF9DCBBAC55A06295CE870B07029BFCDB2DCE28D959F2815B16F81798) + +"""Prime order of secp256k1""" +N = C.order() + +"""Finite field of scalars of secp256k1""" +Z = FiniteField(N) + +""" Beta value of secp256k1 non-trivial endomorphism: lambda * (x, y) = (beta * x, y)""" +BETA = F(2)^((P-1)/3) + +""" Lambda value of secp256k1 non-trivial endomorphism: lambda * (x, y) = (beta * x, y)""" +LAMBDA = Z(3)^((N-1)/3) + +assert is_prime(P) +assert is_prime(N) + +assert BETA != F(1) +assert BETA^3 == F(1) +assert BETA^2 + BETA + 1 == 0 + +assert LAMBDA != Z(1) +assert LAMBDA^3 == Z(1) +assert LAMBDA^2 + LAMBDA + 1 == 0 + +assert Integer(LAMBDA)*G == C(BETA*G[0], G[1])