From 96fd2c794e31381339fbf1d9372f9f11cd945691 Mon Sep 17 00:00:00 2001 From: Daira Hopwood Date: Tue, 29 Dec 2020 12:53:13 +0000 Subject: [PATCH] [WIP] Add a prototype implementation of hash-to-curve. This intends to implement the Internet Draft but has not been checked. Signed-off-by: Daira Hopwood --- hashtocurve.sage | 232 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100755 hashtocurve.sage diff --git a/hashtocurve.sage b/hashtocurve.sage new file mode 100755 index 0000000..1113853 --- /dev/null +++ b/hashtocurve.sage @@ -0,0 +1,232 @@ +#!/usr/bin/env sage + +# Simplified SWU for a = 0 as described in [WB19] and +# . + +import sys +from math import ceil, log +from struct import pack + +import hashlib +if sys.version_info < (3, 6): + try: + import sha3 + except ImportError: + print('Please run:\n`sage -c "import sys; print(sys.executable)"` -m pip install pysha3\n') + raise +from hashlib import shake_128 + +if sys.version_info[0] == 2: + range = xrange + as_byte = ord +else: + as_byte = lambda x: x + + +class Cost: + def __init__(self, sqrs=0, muls=0, invs=0): + self.sqrs = sqrs + self.muls = muls + self.invs = invs + + def sqr(self, x): + self.sqrs += 1 + return x^2 + + def mul(self, x, y): + self.muls += 1 + return x * y + + def div(self, x, y): + self.invs += 1 + self.muls += 1 + return x / y + + def inv0(self, x): + self.invs += 1 + return 0 if x == 0 else x^-1 + + def sqrt(self, x): + self.sqrs += 247 + self.muls += 35 + return x.sqrt() if x.is_square() else 0 + + def __add__(self, other): + return Cost(self.sqrs + other.sqrs, self.muls + other.muls, self.invs + other.invs) + + def __repr__(self): + return "%dS + %dM + %dI" % (self.sqrs, self.muls, self.invs) + +# E: a short Weierstrass elliptic curve +def find_z_sswu(E): + (0, 0, 0, A, B) = E.a_invariants() + F = E.base_field() + + R. = F[] # Polynomial ring over F + g = x^3 + F(A) * x + F(B) # y^2 = g(x) = x^3 + A * x + B + ctr = F.gen() + while True: + for Z_cand in (F(ctr), F(-ctr)): + if Z_cand.is_square(): + # Criterion 1: Z is non-square in F. + continue + if Z_cand == F(-1): + # Criterion 2: Z != -1 in F. + continue + if not (g - Z_cand).is_irreducible(): + # Criterion 3: g(x) - Z is irreducible over F. + continue + if g(B / (Z_cand * A)).is_square(): + # Criterion 4: g(B / (Z * A)) is square in F. + return Z_cand + ctr += 1 + + +p = 0x40000000000000000000000000000000224698fc094cf91b992d30ed00000001 +E_isop = EllipticCurve(GF(p), [10949663248450308183708987909873589833737836120165333298109615750520499732811, 1265]) +E_p = EllipticCurve(GF(p), [0, 5]) +Z_isop = find_z_sswu(E_isop) +assert Z_isop == Mod(-13, p) + +k = 128 +L = (len(format(p, 'b')) + k + 7) // 8 +assert L == 48 + +CONSTANT_TIME = True + +def select_z_nz(s, ifz, ifnz): + # This should be constant-time in a real implementation. + return ifz if (s == 0) else ifnz + +def map_to_curve_simple_swu(E, Z, u, c): + # would be precomputed + (0, 0, 0, A, B) = E.a_invariants() + mBdivA = -B / A + BdivZA = B / (Z * A) + #print("A = ", A) + #print("B = ", B) + #print("Z = ", Z) + #print("-B/A = ", mBdivA) + #print("B/ZA = ", BdivZA) + + # 1. tv1 = inv0(Z^2 * u^4 + Z * u^2) + Z2 = c.sqr(Z) + u2 = c.sqr(u) + u4 = c.sqr(u2) + ta = c.mul(Z2, u4) + c.mul(Z, u2) + tv1 = c.inv0(ta) + + # 2. x1 = (-B / A) * (1 + tv1) + # 3. If tv1 == 0, set x1 = B / (Z * A) + x1 = select_z_nz(tv1, BdivZA, mBdivA * (1 + tv1)) + + # 4. gx1 = x1^3 + A * x1 + B + # = x1*(x1^2 + A) + B + x1_2 = c.sqr(x1) + gx1 = c.mul(x1, x1_2 + A) + B + + # 5. x2 = Z * u^2 * x1 + tb = c.mul(Z, u2) + x2 = c.mul(tb, x1) + + # 6. gx2 = x2^3 + A * x2 + B + # = x2*(x2^2 + A) + B + x2_2 = c.sqr(x2) + gx2 = c.mul(x2, x2_2 + A) + B + + # 7. If is_square(gx1), set x = x1 and y = sqrt(gx1) + # 8. Else set x = x2 and y = sqrt(gx2) + y1 = c.sqrt(gx1) + y1_2 = c.sqr(y1) + if CONSTANT_TIME or y1_2 != gx1: + y2 = c.sqrt(gx2) + x = select_z_nz(y1_2 - gx1, x1, x2) + y = select_z_nz(y1_2 - gx1, y1, y2) + else: + (x, y) = (x1, y1) + + # 9. If sgn0(u) != sgn0(y), set y = -y + y = select_z_nz((int(u) % 2) - (int(y) % 2), y, -y) + + return E((x, y)) + + +# iso_Ep = Isogeny of degree 3 from Elliptic Curve defined by y^2 = x^3 + 10949663248450308183708987909873589833737836120165333298109615750520499732811*x + 1265 over Fp +def iso_map(x, y, c): + c.muls += 2+1+1 + 2+1+1+2 + # batch inversion + c.muls += 3 + c.invs += 1 + return (((( 6432893846517566412420610278260439325191790329320346825767705947633326140075 *x + + 23989696149150192365340222745168215001509815558210986772351135915822265203574)*x + + 10492611921771203378452795982353351666191589197598957448093274638589204800759)*x + + 12865787693035132824841220556520878650383580658640693651535411895266652280192) / + (( x + + 13271109177048389296812780941310096270046944650307955939477485891950613419807)*x + + 22768321103861051515190775253992702316905399997697804654926324362758820947460), + (((11793638718615538422771118843477472096184948937087302513907460903994431256804 *x + + 11994848074575096182670111372584107500754907779105493386175567957911132601787)*x + + 28823569610051396102362669851238297121581474897215657071023781420043761726004)*x + + 1072148974419594402070101713043406554198631721553391137627950991272221023311) * y / + ((( x + + 5432652610908059517272798285879155923388888734491153551238890455750936314542)*x + + 10408918692925056833786833257634153023990087029210292532869619559576527581706)*x + + 28948022309329048855892746252171976963363056481941560715954676764349967629797)) + + +def expand_message_xof(msg, DST, len_in_bytes): + assert len(DST) < 256 + len_in_bytes = int(len_in_bytes) + + # This is horrible but matches the reference code. + xof = shake_128() + xof.update(msg) + xof.update(pack(">H", len_in_bytes)) + xof.update(pack("B", len(DST))) + xof.update(DST) + return xof.digest(len_in_bytes) + +def hash_to_field(msg, DST, count): + uniform_bytes = expand_message_xof(msg, DST, L*count) + return [Mod(OS2IP(uniform_bytes[L*i : L*(i+1)]), p) for i in range(count)] + +def OS2IP(bs): + acc = 0 + for b in bs: + acc = (acc<<8) + as_byte(b) + return acc + +def hash_to_curve(msg, DST, uniform=True): + c = Cost() + u = hash_to_field(msg, DST, 2) + #print("u = ", u) + R = map_to_curve_simple_swu(E_isop, Z_isop, u[0], c) + + if uniform: + Q1 = map_to_curve_simple_swu(E_isop, Z_isop, u[1], c) + + # We could batch the two inv0 inversions (not done above for simplicity). + c.invs -= 1 + c.muls += 3 + + # Complete addition using affine coordinates: I + 2M + 2S + # (S for x1^2; compute numerator and denominator of the division for the correct case; + # I + M to divide; S + M to compute x and y of the result.) + R = R + Q1 + #print("R = ", R) + c.invs += 1 + c.sqrs += 2 + c.muls += 2 + + # no cofactor clearing needed since Pallas and Vesta are prime-order + (x, y) = R.xy() + P = E_p(iso_map(x, y, c)) + return (P, c) + + +#print(hash_to_curve("hello", "blah")) + +iters = 100 +for i in range(iters): + (res, cost) = hash_to_curve(pack(">I", i), "blah", uniform=False) + print(res, cost)