#! /usr/bin/env python3 # This file is dual-licensed. Choose whichever licence you want from # the two licences listed below. # # The first licence is a regular 2-clause BSD licence. The second licence # is the CC-0 from Creative Commons. It is intended to release Monocypher # to the public domain. The BSD licence serves as a fallback option. # # SPDX-License-Identifier: BSD-2-Clause OR CC0-1.0 # # ------------------------------------------------------------------------ # # Copyright (c) 2020, Loup Vaillant and Andrew Moon # All rights reserved. # # # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions are # met: # # 1. Redistributions of source code must retain the above copyright # notice, this list of conditions and the following disclaimer. # # 2. Redistributions in binary form must reproduce the above copyright # notice, this list of conditions and the following disclaimer in the # documentation and/or other materials provided with the # distribution. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS # "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR # A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT # HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, # SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT # LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, # DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY # THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. # # ------------------------------------------------------------------------ # # Written in 2020 by Loup Vaillant and Andrew Moon # # To the extent possible under law, the author(s) have dedicated all copyright # and related neighboring rights to this software to the public domain # worldwide. This software is distributed without any warranty. # # You should have received a copy of the CC0 Public Domain Dedication along # with this software. If not, see # #################### # Field arithmetic # #################### class fe: """Prime field over 2^255 - 19""" p = 2**255 - 19 def __init__(self, x): self.val = x % self.p # Basic arithmetic operations def __neg__ (self ): return fe(-self.val ) def __add__ (self, o): return fe( self.val + o.val ) def __sub__ (self, o): return fe( self.val - o.val ) def __mul__ (self, o): return fe((self.val * o.val ) % self.p) def __truediv__ (self, o): return fe((self.val * o.invert().val) % self.p) def __floordiv__(self, o): return fe( self.val // o ) def __pow__ (self, s): return fe(pow(self.val, s , self.p)) def invert (self ): return fe(pow(self.val, self.p-2, self.p)) def __eq__(self, other): return self.val % self.p == other.val % self.p def __ne__(self, other): return self.val % self.p != other.val % self.p def is_positive(self) : return self.val % self.p <= (p-1) // 2 def is_negative(self) : return self.val % self.p > (p-1) // 2 def abs(self): if self.is_positive(): return self else : return -self def print(self): """prints a field element in little endian""" m = self.val % self.p for _ in range(32): print(format(m % 256, '02x'), end='') m //= 256 if m != 0: raise ValueError('number is too big!!') print(':') def print_raw(raw): """prints a raw element in little endian""" for _ in range(32): print(format(raw % 256, '02x'), end='') raw //= 256 if raw != 0: raise ValueError('number is too big!!') print(':') ######################## # Curve25519 constants # ######################## p = fe.p A = fe(486662) # B = 1 ############### # Square root # ############### # Legendre symbol: # - 0 if n is zero # - 1 if n is a non-zero square # - -1 if n is not a square # We take for granted that n^((p-1)/2) does what we want def chi (n): return n**((p-1)//2) def is_square(n): return n == fe(0) or chi(n) == fe(1) # square root of -1 sqrtm1 = (fe(2)**((p-1) // 4)).abs() if sqrtm1 * sqrtm1 != fe(-1): raise ValueError('Wrong sqrtm1') # The square root of n, if n is a square. # # Note that p is congruent to 5 modulo 8, so (p+3)/8 is an integer. # If n is zero, then n^((p+3)/8) is zero (zero is its own square root). # Otherwise: # (n^((p+3)/8))^4 = n^((p+3)/2) # (n^((p+3)/8))^4 = n^((p-1)/2) * n^2 # (n^((p+3)/8))^4 = chi(n) * n^2 -- chi(n) == 1 # (n^((p+3)/8))^4 = n^2 -- because n is a non-zero square # (n^((p+3)/8))^2 = n or -n # case n: # (n^((p+3)/8))^2 = n # n^((p+3)/8) = sqrt(n) or -sqrt(n) # case -n: # (n^((p+3)/8))^2 = -n # -1 * (n^((p+3)/8))^2 = n # sqrt(-1) * n^((p+3)/8) = sqrt(n) or -sqrt(n) # # We then choose the positive square root, between 0 and (p-1)/2 def sqrt(n): if not is_square(n) : raise ValueError('Not a square!') root = n**((p+3) // 8) if root * root != n: root = (root * sqrtm1) if root * root != n: raise ValueError('Should be a square!!') return root.abs() ########################### # Elligator 2 (reference) # ########################### # Elligator: Elliptic-curve points indistinguishable from uniform random strings # by Daniel J. Bernstein, Mike Hamburg, Anna Krasnova, and Tanja Lange # 2013 # https://elligator.cr.yp.to/ # Arbitrary non square, typically chosen to minimise computation. # 2 and sqrt(-1) both work fairly well, but 2 seems to be more popular. # We stick to 2 for compatibility. non_square = fe(2) # Representative to curve point, straight from the paper. # Unlike the paper, curve coordinates are called (u, v) to follow # established conventions. Thus, "v" in the paper is called "w" here. def hash_to_curve(r): w = -A / (fe(1) + non_square * r**2) e = chi(w**3 + A*w**2 + w) u = e*w - (fe(1)-e)*(A//2) v = -e * sqrt(u**3 + A*u**2 + u) return (u, v) # Test whether a point has a representative, straight from the paper. def can_curve_to_hash(u): return u != -A and is_square(-non_square * u * (u+A)) # Computes the representative of a point, straight from the paper. def curve_to_hash(u, v_is_negative): if not can_curve_to_hash(u): raise ValueError('cannot curve to hash') sq1 = sqrt(-u / (non_square * (u+A))) sq2 = sqrt(-(u+A) / (non_square * u )) if v_is_negative: return sq2 else : return sq1 ##################### # Elligator2 (fast) # ##################### # Inverse square root. # Returns (0 , True ) if x is zero. # Returns (sqrt(1/x) , True ) if x is non-zero square. # Returns (sqrt(sqrt(-1)/x), False) if x is not a square. # We do not guarantee the sign of the square root. # # Notes: # Let quartic = x^((p-1)/4) # # x^((p-1)/2) = chi(x) # quartic^2 = chi(x) # quartic = sqrt(chi(x)) # quartic = 1 or -1 or sqrt(-1) or -sqrt(-1) # # Note that x is a square if quartic is 1 or -1 # There are 4 cases to consider: # # if quartic = 1 (x is a square) # then x^((p-1)/4) = 1 # x^((p-5)/4) * x = 1 # x^((p-5)/4) = 1/x # x^((p-5)/8) = sqrt(1/x) or -sqrt(1/x) # # if quartic = -1 (x is a square) # then x^((p-1)/4) = -1 # x^((p-5)/4) * x = -1 # x^((p-5)/4) = -1/x # x^((p-5)/8) = sqrt(-1) / sqrt(x) # x^((p-5)/8) * sqrt(-1) = sqrt(-1)^2 / sqrt(x) # x^((p-5)/8) * sqrt(-1) = -1/sqrt(x) # x^((p-5)/8) * sqrt(-1) = -sqrt(1/x) or sqrt(1/x) # # if quartic = sqrt(-1) (x is not a square) # then x^((p-1)/4) = sqrt(-1) # x^((p-5)/4) * x = sqrt(-1) # x^((p-5)/4) = sqrt(-1)/x # x^((p-5)/8) = sqrt(sqrt(-1)/x) or -sqrt(sqrt(-1)/x) # # Note that the product of two non-squares is always a square: # For any non-squares a and b, chi(a) = -1 and chi(b) = -1. # Since chi(x) = x^((p-1)/2), chi(a)*chi(b) = chi(a*b) = 1. # Therefore a*b is a square. # # Since sqrt(-1) and x are both non-squares, their product is a # square, and we can compute their square root. # # if quartic = -sqrt(-1) (x is not a square) # then x^((p-1)/4) = -sqrt(-1) # x^((p-5)/4) * x = -sqrt(-1) # x^((p-5)/4) = -sqrt(-1)/x # x^((p-5)/8) = sqrt(-sqrt(-1)/x) # x^((p-5)/8) = sqrt( sqrt(-1)/x) * sqrt(-1) # x^((p-5)/8) * sqrt(-1) = sqrt( sqrt(-1)/x) * sqrt(-1)^2 # x^((p-5)/8) * sqrt(-1) = sqrt( sqrt(-1)/x) * -1 # x^((p-5)/8) * sqrt(-1) = -sqrt(sqrt(-1)/x) or sqrt(sqrt(-1)/x) def invsqrt(x): isr = x**((p - 5) // 8) quartic = x * isr**2 if quartic == fe(-1) or quartic == -sqrtm1: isr = isr * sqrtm1 is_square = quartic == fe(1) or quartic == fe(-1) or x == fe(0) return isr, is_square # From the paper: # w = -A / (fe(1) + non_square * r^2) # e = chi(w^3 + A*w^2 + w) # u = e*w - (fe(1)-e)*(A//2) # v = -e * sqrt(u^3 + A*u^2 + u) # # Note that e is either 0, 1 or -1 # if e = 0 # (u, v) = (0, 0) # if e = 1 # u = w # v = -sqrt(u^3 + A*u^2 + u) # if e = -1 # u = -w - A = w * non_square * r^2 # v = sqrt(u^3 + A*u^2 + u) # # Let r1 = non_square * r^2 # Let r2 = 1 + r1 # Note that r2 cannot be zero, -1/non_square is not a square. # We can (tediously) verify that: # w^3 + A*w^2 + w = (A^2*r1 - r2^2) * A / r2^3 # Therefore: # chi(w^3 + A*w^2 + w) = chi((A^2*r1 - r2^2) * (A / r2^3)) # chi(w^3 + A*w^2 + w) = chi((A^2*r1 - r2^2) * (A / r2^3)) * 1 # chi(w^3 + A*w^2 + w) = chi((A^2*r1 - r2^2) * (A / r2^3)) * chi(r2^6) # chi(w^3 + A*w^2 + w) = chi((A^2*r1 - r2^2) * (A / r2^3) * r2^6) # chi(w^3 + A*w^2 + w) = chi((A^2*r1 - r2^2) * A * r2^3) # Corollary: # e = 1 if (A^2*r1 - r2^2) * A * r2^3) is a non-zero square # e = -1 if (A^2*r1 - r2^2) * A * r2^3) is not a square # Note that w^3 + A*w^2 + w (and therefore e) can never be zero: # w^3 + A*w^2 + w = w * (w^2 + A*w + 1) # w^3 + A*w^2 + w = w * (w^2 + A*w + A^2/4 - A^2/4 + 1) # w^3 + A*w^2 + w = w * (w + A/2)^2 - A^2/4 + 1) # which is zero only if: # w = 0 (impossible) # (w + A/2)^2 = A^2/4 - 1 (impossible, because A^2/4-1 is not a square) # # Let isr = invsqrt((A^2*r1 - r2^2) * A * r2^3) # isr = sqrt(1 / ((A^2*r1 - r2^2) * A * r2^3)) if e = 1 # isr = strt(sqrt(-1) / ((A^2*r1 - r2^2) * A * r2^3)) if e = -1 # # if e = 1 # let u1 = -A * (A^2*r1 - r2^2) * A * r2^2 * isr^2 # u1 = w # u1 = u # let v1 = -(A^2*r1 - r2^2) * A * isr # v1 = -sqrt((A^2*r1 - r2^2) * A / r2^3) # v1 = -sqrt(w^3 + A*w^2 + w) # v1 = -sqrt(u^3 + A*u^2 + u) (because u = w) # v1 = v # # if e = -1 # let ufactor = -non_square * sqrt(-1) * r^2 # let vfactor = sqrt(ufactor) # let u2 = -A * (A^2*r1 - r2^2) * A * r2^2 * isr^2 * ufactor # u2 = w * -1 * -non_square * r^2 # u2 = w * non_square * r^2 # u2 = u # let v2 = (A^2*r1 - r2^2) * A * isr * vfactor # v2 = sqrt(non_square * r^2 * (A^2*r1 - r2^2) * A / r2^3) # v2 = sqrt(non_square * r^2 * (w^3 + A*w^2 + w)) # v2 = sqrt(non_square * r^2 * w * (w^2 + A*w + 1)) # v2 = sqrt(u (w^2 + A*w + 1)) # v2 = sqrt(u ((-u-A)^2 + A*(-u-A) + 1)) # v2 = sqrt(u (u^2 + A^2 + 2*A*u - A*u -A^2) + 1)) # v2 = sqrt(u (u^2 + A*u + 1)) # v2 = sqrt(u^3 + A*u^2 + u) # v2 = v ufactor = -non_square * sqrtm1 vfactor = sqrt(ufactor) def fast_hash_to_curve(r): t1 = r**2 * non_square # r1 u = t1 + fe(1) # r2 t2 = u**2 t3 = (A**2 * t1 - t2) * A # numerator t1 = t2 * u # denominator t1, is_square = invsqrt(t3 * t1) u = r**2 * ufactor v = r * vfactor if is_square: u = fe(1) if is_square: v = fe(1) v = v * t3 * t1 t1 = t1**2 u = u * -A * t3 * t2 * t1 if is_square != v.is_negative(): # XOR v = -v return (u, v) # From the paper: # Let sq = -non_square * u * (u+A) # if sq is not a square, or u = -A, there is no mapping # Assuming there is a mapping: # if v is positive: r = sqrt(-u / (non_square * (u+A))) # if v is negative: r = sqrt(-(u+A) / (non_square * u )) # # We compute isr = invsqrt(-non_square * u * (u+A)) # if it wasn't a square, abort. # else, isr = sqrt(-1 / (non_square * u * (u+A)) # # If v is positive, we return isr * u: # isr * u = sqrt(-1 / (non_square * u * (u+A)) * u # isr * u = sqrt(-u / (non_square * (u+A)) # # If v is negative, we return isr * (u+A): # isr * (u+A) = sqrt(-1 / (non_square * u * (u+A)) * (u+A) # isr * (u+A) = sqrt(-(u+A) / (non_square * u) def fast_curve_to_hash(u, v_is_negative): t = u + A r = -non_square * u * t isr, is_square = invsqrt(r) if not is_square: return None if v_is_negative: u = t r = u * isr r = r.abs() return r