I'm writing a Python script for testing C implementation of Edwards elliptic curves point addition and doubling. I'm following this paper and implementing the formula set (5) for "unified addition" (which means any 2 points can be added), and formula set (7) for "dedicated doubling" (which is more efficient than (5) if 2 points are the same). But my code doesn't seem to compute the doubling correctly. Maybe it's also addition, which I have no way of telling, since I don't have another reference to compare to. #!/usr/bin/env python3 import sys, secrets from functools import reduce def os2int(v): ret = 0 for b in v: ret = (ret << 8) | b return ret # curve parameters from https://datatracker.ietf.org/doc/html/rfc8032#section-5.1 p = (1 << 255) - 19 a = -1 d_over = -121665 d_under = 121666 d = d_over * pow(d_under, p-2, p) % p # product of d_over with the modular inverse of d_under mod p def point_add_ref(x1, y1, t1, z1, x2, y2, t2, z2): x1y2 = x1 * y2 % p x2y1 = x2 * y1 % p x1x2 = x1 * x2 % p y1y2 = y1 * y2 % p z1z2 = z1 * z2 % p t1t2 = t1 * t2 % p x3 = (x1y2 + x2y1) * (z1z2 - d * t1t2) % p y3 = (y1y2 - a * x1x2) * (z1z2 + d * t1t2) % p t3 = (y1y2 - a * x1x2) * (x1y2 + x2y1) % p z3 = (z1z2 - d * t1t2) * (z1z2 + d * t1t2) % p return (x3, y3, t3, z3) def point_dbl_ref(x1, y1, t1, z1): xx = x1 * x1 % p yy = y1 * y1 % p zz = z1 * z1 % p xy = x1 * y1 % p t = 2 * xy % p u = (yy + a * xx) % p v = (yy - a * xx) % p w = (2 * zz - yy - a * xx) % p x3 = t * w % p y3 = u * v % p t3 = t * v % p z3 = u * w % p return (x3, y3, t3, z3) def xytz_cmp(P,Q): vec = (P * Q[3] % p != P[3] * Q % p for i in range(3)) return reduce(lambda a, b: a or b, vec) if __name__ == "__main__": fails = 0 slen = 12 for i in range(100): P = (os2int(secrets.token_bytes(slen)), os2int(secrets.token_bytes(slen)), os2int(secrets.token_bytes(slen)), os2int(secrets.token_bytes(slen))) R3 = point_add_ref(*P, *P) R4 = point_dbl_ref(*P) if xytz_cmp(R3, R4): fails += 1 print("{} test(s) failed.".format(fails)) Continue reading...