McCallum-Relyea exchange on x25519 elliptic curve

Background

I have a degree in electronics engineering, I'm not bad at programming with many languages, but I'm not that good at cryptography, I'm just able to use that from libraries, and I understand just a few basic principles. I'm evaluating the McCallum-Relyea key exchange protocol (the one used in the tang protocol) in order to use it in a project, so I became aware of a few things about elliptic curve cryptography, further details can be found below. I apologize in advance for the atrocities I'm surely going to say regarting this subject.

The problem

I would like to use the protocol to derive an encryption key used for storing other sensitive keys on a device, allowing decryption only if it can be trusted (for now suppose it's a manual process, the server waits for confirmation before responding to the device's request). The point is protecting those sensitive keys in case of theft. So far I implemented a simple test using Python and tinyec as below, however I need to do the same on top of the X25519 protocol, but I can't find a library with point addition and subtraction methods.

Example Python code with tinyec, working (as far as I can say in my limited knowledge):

# McCallum-Relyea exchange simple test
from tinyec import registry
from tinyec.ec import *

curve = registry.get_curve('brainpoolP256r1')
server = make_keypair(curve)

# provisioning step
client = make_keypair(curve)  # public must remain secret
derived = ECDH(client).get_secret(server)
client.priv = None  # discarded
client.can_sign = False
print(derived.x, derived.y)

# recovery step
ephemeral = make_keypair(curve)  # both must be kept secret
temp = client.pub + ephemeral.pub
# send temp key to server
temp = ECDH(server).get_secret(Keypair(curve, priv=None, pub=temp))  # server side operation
# get back updated temp key from server
recovered = temp - ECDH(ephemeral).get_secret(server)
print(recovered.x, recovered.y)

print("Keys%s match." % ("" if recovered == derived else " do not"))

As far as I've been able to understand, in the X25519 protocol only the $x$ coordinate is used, and I found formulas for adding points using only their $x$ coordinates but they also require to know the difference between those points (I think because they are only used in Montgomery ladders), thus I haven't been able to use those. So I tried to calculate the $y$ coordinate using the actual curve equation, then selecting the sign for the calculated $y$ based on the sign of the input $x$. Then I tried to add or subtract two points (subtraction should be the same as addition with flipped $y$ coordinate) using the affine coordinates formulas.

Edit: here is the updated code based on Daniel's suggestion, and there was an error in addition/subtraction too, now the exchange works half the time because there is still the problem of selecting the correct root for $y$ when decoding the point.

So the problem now is: how to make it work the other half times? (partially solved)

And is there a way of computing $\sqrt{y}$ using exponentiation instead?

Or even better is there a way of adding/subtracting without having to recover $y$? It seems that there is no indication of the sign to pick for $y$ in $x$ (the MSB seems always 0).

Edit 2: it seems to be working in any case actually, in fact, if both addition and subtraction are computed at the final step then either one of them has to be the correct result. But still, the problem remains: how to determine in advance which one to compute?

from secrets import randbelow
from curve25519 import *
from libnum import sqrtmod_prime_power

def decodePoint(x):
    x = unpack_number(x)
    y = [y for y in sqrtmod_prime_power(x ** 3 + A * (x ** 2) + x, P, 1)]
    y = y[-1]  # just picking a random one since I don't know which one to pick
    return x % P, y % P


def groupAddSub(a, b, sub):
    xa, ya = decodePoint(a)
    xb, yb = decodePoint(b)
    yb = -yb if sub else yb
    x = pow(yb - ya, 2, P) * pow(xb - xa, -2, P) - A - xa - xb
    return pack_number(x % P)


server = X25519PrivateKey(fix_secret(randbelow(P)))

# provisioning
client = X25519PrivateKey(fix_secret(randbelow(P)))  # public must remain secret
derived = client.exchange(server.public_key())
client = client.public_key()
print(derived.hex())

# recovery
ephemeral = X25519PrivateKey(fix_secret(randbelow(P)))  # both must be kept secret
temp = groupAddSub(client.public_bytes(), ephemeral.public_key_bytes(), False)
# send temp key to server
temp = server.exchange(temp)  # server side
# get back updated temp key from server
recoveredAdd = groupAddSub(temp, ephemeral.exchange(server.public_key()), False)
recoveredSub = groupAddSub(temp, ephemeral.exchange(server.public_key()), True)

if derived == recoveredAdd:
    print(recoveredAdd.hex())
elif derived == recoveredSub:
    print(recoveredSub.hex())
else:
    print("no match")

Please note that as of now I still haven't considered security issues like constant time execution, I'm still trying to understand the process, let's consider one issue at a time. And here is the curve25519.py file I've been using, I just removed the trailing _ from functions containing it to use them in my script.

Context details

This is not relevant to the question, but here is the full story anyway (suggestions are always welcome). I'm trying to implement an access control system, since it has to work with NFC tags supporting only 3DES authentication, I have to store their keys on the reader, and being symmetric keys I have to encrypt them "at rest". The reader device has to work even offline, except for the initial connection where it decrypts its main key, keeping it in volatile memory so if it gets stolen it will still be safe. The one I'm discussing is the only applicable mechanism I've been able to find, but again I may be missing a lot of better ways to do it since I have a limited knowledge on the subject.

I think that your main issue is in the use of the math.isqrt (there may or may not be others). The function returns the floor of the positive, real square root of a positive number i.e. you are computing $$y=\lfloor\sqrt{x^3+Ax^2+x}\rfloor.$$ Such a $y$ is unlikely to satisfy the elliptic curve equation modulo $p$. Instead you should be computing a modular square root such as sqrtmod_prime_power in the libnum library. Calling libnum.sqrtmod_prime_power(x**3+A*x**2+x,P,1) will return a $y$ such that the curve equation is satisfied modulo $p$, but you may need to adjust the circumstances under which you negate $y$ to account for the high bit of $y$.