permissible.py

# Jacobi symbol calculation from https://eprint.iacr.org/2021/1271.pdf

def batch_matrix (x, y, delta, batch):
    """Compute matrix using (batch + 2) bits of (x,y)"""
    (ai, bi, ci, di, u) = (1, 0, 0, 1, 0)
    for i in range (batch):
        yi = y
        if delta >= 0 and x & 1:
            (delta, x, y) = (-delta, (x - y) >> 1, x)
            (ai, bi, ci, di) = (ai - ci, bi - di, 2 * ai, 2 * bi)
        elif x & 1:
            (delta, x, y) = (1 + delta, (x + y) >> 1, y)
            (ai, bi, ci, di) = (ai + ci, bi + di, 2 * ci, 2 * di)
        else:
            (delta, x, y) = (1 + delta, x >> 1, y)
            (ai, bi, ci, di) = (ai, bi, 2 * ci, 2 * di)
        u += ((yi & y) ^ (y >> 1)) & 2 # quadratic reciprocity
        u += (u & 1) ^ int(ci < 0) # count sign changes
        u %= 4
    return (u, delta, (ai, bi, ci, di))

def jacobi (x, y, batch = 32):
    """Return jacobi symbol (x on y)"""
    assert (y > 0) and (y & 1) # y-input must be positive and odd
    x, delta, t = x % y, 0, 0
    # Compute number of iterations per:
    # https://eprint.iacr.org/2021/549.pdf page 14
    # https://github.com/sipa/safegcd-bounds
    nbits = y.bit_length()
    niters = (45907 * nbits + 26313) // 19929
    mask = (1 << (batch + 2)) - 1 # low bits
    for i in range (0, niters, batch):
        u, delta, (ai, bi, ci, di) = \
            batch_matrix(x & mask, y & mask, delta, batch)
        (x, y) = ((ai * x + bi * y) >> batch, (ci * x + di * y) >> batch)
        t = (t + u) % 4
        t = (t + ((t & 1) ^ int(y < 0))) % 4
    t = (t + (t & 1)) % 4 # snap to [0 ,2]
    if y in [-1, 1]: jacobi = 1 - t
    else : jacobi = 0 # gcd != 1
    return jacobi

p = 2**255-19

def expmod(b, e, m):
    if e == 0: return 1
    t = expmod(b, e//2, m)**2 % m
    if e & 1: t = (t * b) % m
    return t

# Finite field mod p
class FieldElement:
    def __init__(self, val):
        self.val = val
    def __add__(self, other):
        return FieldElement((self.val + other.val) % p)
    def __sub__(self, other):
        return FieldElement((p + self.val - other.val) % p)
    def __mul__(self, other):
        return FieldElement((self.val * other.val) % p)
    def __neg__(self):
        return FieldElement((p - self.val) % p)
    def inv(self):
        return FieldElement(expmod(self.val, p - 2, p))
    def legendre(self):
        return jacobi(self.val, p)
    def __truediv__(self, other):
        return self * other.inv()
    def __eq__(self, other):
        return self.val == other.val
    def __str__(self):
        return str(self.val)

c = FieldElement(51042569399160536130206135233146329284152202253034631822681833788666877215207)
d = -FieldElement(121665) / FieldElement(121666)

# d is non-square
assert d.legendre() == -1

# Twisted Edwards in projective coordinates
class GroupElement:
    def __init__(self, x, y, z = FieldElement(1)):
        self.x = x
        self.y = y
        self.z = z
    def __add__(self, other):
        X1, Y1, Z1 = self.x, self.y, self.z
        X2, Y2, Z2 = other.x, other.y, other.z
        # Formula from https://hyperelliptic.org/EFD/g1p/auto-twisted-projective.html#addition-add-2008-bbjlp
        A = Z1 * Z2
        B = A * A
        C = X1 * X2
        D = Y1 * Y2
        E = d * C * D
        F = B - E
        G = B + E
        X3 = A * F * ((X1 + Y1) * (X2 + Y2) - C - D)
        Y3 = A * G * (D + C)
        Z3 = F * G
        return GroupElement(X3, Y3, Z3)
    def __neg__(self):
        return GroupElement(-self.x, self.y, self.z)
    def __sub__(self, other):
        return self + (-other)
    def normalize(self):
        invz = self.z.inv()
        X = self.x * invz
        Y = self.y * invz
        return GroupElement(X, Y)
    def __eq__(self, other):
        return self.x * other.z == other.x * self.z and self.y * other.z == other.y * self.z
    def __str__(self):
        return "[%s:%s:%s]" % (self.x, self.y, self.z)

def is_permissible(P):
    # Convert to projective Montgomery form (U, V, W) where u = U/W and v = V/W
    # Note: v-coordinate in Montgomery form (Curve25519) is equal to the y-coordinate in short Weierstrass form (Wei25519).
    # U is not needed to evaluate permissibility
    # U = P.x * (P.z + P.y)
    V = c * P.z * (P.z + P.y)
    W = P.x * (P.z - P.y)
    lW = W.legendre()
    # check that (1 + v) is square
    if (W + V).legendre() != lW:
        return False
    # check that (1 - v) is non-square
    if (W - V).legendre() != -lW:
        return False
    return True

def make_permissible(P, Q):
    """Add Q to P until the result is permissible"""
    i = 0
    while not is_permissible(P):
        P = P + Q
        i = i + 1
    P = P.normalize()
    return (P, i)

Identity = GroupElement(FieldElement(0), FieldElement(1))
G = GroupElement(FieldElement(15112221349535400772501151409588531511454012693041857206046113283949847762202), \
    FieldElement(46316835694926478169428394003475163141307993866256225615783033603165251855960))
T = GroupElement(FieldElement(31919686468363732305077775836587378998231646534201478966360489369744205619), \
    FieldElement(37881138124249956567716014034273041038361541292759991175595923373144531920593))

assert (G - G) == Identity

Gp, it = make_permissible(G, T)
print("Permissible point: G + %i T" % it)
print(Gp)
assert is_permissible(Gp)
assert not is_permissible(-Gp)