Python function for Edwards curves point doubling and addition

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[i] * Q[3] % p != P[3] * Q[i] % 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))

It's neither the addition nor the doubling that's incorrect. It's that you're generating points with random coordinates which are not valid curve points