A = fe(486662)
sqrtm1 = (fe(2)**((p - 1) // 4)).abs()
u = fe(2)
ufactor = -u * sqrtm1
ufactor_sqrt = (ufactor)**((p + 3) // 8)

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)
    return isr, is_square

def fast_hash_to_curve(q):
    r = u * q**2
    r1 = (r + fe(1))
    num = A * (A**2 * r - r1**2)
    den = r1**3
    # x = -A / (r + 1)
    # y = x^3 + A*x^2 + x
    # y = A^3/(r + 1)^2 - A^3/(r + 1)^3 - A/(r + 1)
    # y = (A^3*r - A*(r + 1)^2) / (r + 1)^3
    isr, is_square = invsqrt(num * den)
    # if is_square: isr = sqrt(1 / (num * den))
    # if not is_square: isr = sqrt(sqrtm1 / (num * den))
    x = -A * (num * r1**2 * isr**2)
    # x = -A * num * (r + 1)^2 * sqrt(1 / (num * den))^2
    # x = -A * num * (r + 1)^2 * 1 / (num * den)
    # x = -A * (r + 1)^2 * 1 / den
    # x = -A / (r + 1)
    y = num * isr
    # y = num * sqrt(1 / (num * den))
    # y = sqrt(num^2 / (num * den))
    # y = sqrt(num / den)
    qx = q**2 * ufactor
    qy = q * ufactor_sqrt
    if is_square: qx = fe(1)
    if is_square: qy = fe(1)
    x = qx * x
    # x = q^2 * -u * sqrt(-1) * -A * sqrt(-1) / (r + 1)
    # x = -A * u * q^2 / (r + 1)
    # x = -A * r / (r + 1)
    y = qy * y
    # y = q * sqrt(-u * sqrtm1) * sqrt(sqrt(-1) * num / den)
    # y = sqrt(q^2 * u * num / den)
    # y = sqrt(r * num / den)
    y = y.abs()
    if not is_square: y = -y
    return (x, y)

def fast_curve_to_hash(point):
    x, y = point
    t0 = A + x
    t1 = x
    # if is_positive(y): r = u*q^2 = -(A + x)/x
    # if is_negative(y): r = u*q^2 = -x/(A + x)
    isr, is_square = invsqrt(-t0 * t1 * u)
    # isr = sqrt(-1 / ((A + x) * x * u))
    if not is_square:
        return None
    num = t0
    if y.is_positive(): num = t1
    q = num * isr
    # if is_positive(y): q = (A + x) * sqrt(1 / (-x * (A + x) * u)) = sqrt(-(A + x) / (x * u))
    # if is_positive(y): q = sqrt(-(A + x) / (x * u))
    # if is_negative(y): q = x * sqrt(1 / (-x * (A + x) * u)) = sqrt(-x / ((A + x) * u))
    # if is_negative(y): q = sqrt(-x / ((A + x) * u))
    q = q.abs()
    return q