1 """ Generates the constants used in secp256k1_scalar_split_lambda.
2 3 See the comments for secp256k1_scalar_split_lambda in src/scalar_impl.h for detailed explanations.
4 """
5 6 load("secp256k1_params.sage")
7 8 def inf_norm(v):
9 """Returns the infinity norm of a vector."""
10 return max(map(abs, v))
11 12 def gauss_reduction(i1, i2):
13 v1, v2 = i1.copy(), i2.copy()
14 while True:
15 if inf_norm(v2) < inf_norm(v1):
16 v1, v2 = v2, v1
17 # This is essentially
18 # m = round((v1[0]*v2[0] + v1[1]*v2[1]) / (inf_norm(v1)**2))
19 # (rounding to the nearest integer) without relying on floating point arithmetic.
20 m = ((v1[0]*v2[0] + v1[1]*v2[1]) + (inf_norm(v1)**2) // 2) // (inf_norm(v1)**2)
21 if m == 0:
22 return v1, v2
23 v2[0] -= m*v1[0]
24 v2[1] -= m*v1[1]
25 26 def find_split_constants_gauss():
27 """Find constants for secp256k1_scalar_split_lamdba using gauss reduction."""
28 (v11, v12), (v21, v22) = gauss_reduction([0, N], [1, int(LAMBDA)])
29 30 # We use related vectors in secp256k1_scalar_split_lambda.
31 A1, B1 = -v21, -v11
32 A2, B2 = v22, -v21
33 34 return A1, B1, A2, B2
35 36 def find_split_constants_explicit_tof():
37 """Find constants for secp256k1_scalar_split_lamdba using the trace of Frobenius.
38 39 See Benjamin Smith: "Easy scalar decompositions for efficient scalar multiplication on
40 elliptic curves and genus 2 Jacobians" (https://eprint.iacr.org/2013/672), Example 2
41 """
42 assert P % 3 == 1 # The paper says P % 3 == 2 but that appears to be a mistake, see [10].
43 assert C.j_invariant() == 0
44 45 t = C.trace_of_frobenius()
46 47 c = Integer(sqrt((4*P - t**2)/3))
48 A1 = Integer((t - c)/2 - 1)
49 B1 = c
50 51 A2 = Integer((t + c)/2 - 1)
52 B2 = Integer(1 - (t - c)/2)
53 54 # We use a negated b values in secp256k1_scalar_split_lambda.
55 B1, B2 = -B1, -B2
56 57 return A1, B1, A2, B2
58 59 A1, B1, A2, B2 = find_split_constants_explicit_tof()
60 61 # For extra fun, use an independent method to recompute the constants.
62 assert (A1, B1, A2, B2) == find_split_constants_gauss()
63 64 # PHI : Z[l] -> Z_n where phi(a + b*l) == a + b*lambda mod n.
65 def PHI(a,b):
66 return Z(a + LAMBDA*b)
67 68 # Check that (A1, B1) and (A2, B2) are in the kernel of PHI.
69 assert PHI(A1, B1) == Z(0)
70 assert PHI(A2, B2) == Z(0)
71 72 # Check that the parallelogram generated by (A1, A2) and (B1, B2)
73 # is a fundamental domain by containing exactly N points.
74 # Since the LHS is the determinant and N != 0, this also checks that
75 # (A1, A2) and (B1, B2) are linearly independent. By the previous
76 # assertions, (A1, A2) and (B1, B2) are a basis of the kernel.
77 assert A1*B2 - B1*A2 == N
78 79 # Check that their components are short enough.
80 assert (A1 + A2)/2 < sqrt(N)
81 assert B1 < sqrt(N)
82 assert B2 < sqrt(N)
83 84 # Verify connection to Eisenstein integers Z[w] where w = (-1 + sqrt(-3))/2.
85 # The group order N factors as N = pi * conj(pi) in Z[w], where pi = A - B*w
86 # is an Eisenstein prime with norm A^2 + A*B + B^2. The GLV endomorphism
87 # eigenvalue LAMBDA equals B/A mod N, which is the image of w^2 under the
88 # isomorphism Z[w]/(pi) -> Z/NZ (since w -> A/B and (A/B)^2 = B/A in Z/NZ).
89 A_EIS, B_EIS = -B1, A1
90 assert A_EIS**2 + A_EIS*B_EIS + B_EIS**2 == N
91 assert Z(B_EIS / A_EIS) == LAMBDA
92 93 G1 = round((2**384)*B2/N)
94 G2 = round((2**384)*(-B1)/N)
95 96 def rnddiv2(v):
97 if v & 1:
98 v += 1
99 return v >> 1
100 101 def scalar_lambda_split(k):
102 """Equivalent to secp256k1_scalar_lambda_split()."""
103 c1 = rnddiv2((k * G1) >> 383)
104 c2 = rnddiv2((k * G2) >> 383)
105 c1 = (c1 * -B1) % N
106 c2 = (c2 * -B2) % N
107 r2 = (c1 + c2) % N
108 r1 = (k + r2 * -LAMBDA) % N
109 return (r1, r2)
110 111 # The result of scalar_lambda_split can depend on the representation of k (mod n).
112 SPECIAL = (2**383) // G2 + 1
113 assert scalar_lambda_split(SPECIAL) != scalar_lambda_split(SPECIAL + N)
114 115 print(' A1 =', hex(A1))
116 print(' -B1 =', hex(-B1))
117 print(' A2 =', hex(A2))
118 print(' -B2 =', hex(-B2))
119 print(' =', hex(Z(-B2)))
120 print(' -LAMBDA =', hex(-LAMBDA))
121 122 print(' G1 =', hex(G1))
123 print(' G2 =', hex(G2))
124