def descent(A, B):
    """
    Uses Lagrange's method to finda a solution to $w^2 = Ax^2 + By^2$.
    Output a tuple $(x_0, y_0, z_0)$ which a solution to the above equation.
    This solution is used as a base to calculate all the solutions.
    """
    if abs(A) > abs(B):
        y, x, w = descent(B, A)
        return (x, y, w)
    
    if A == 1:
        return (1, 0, 1)
    
    if B == 1:
        return (0, 1, 1)
    
    if quadratic_congruence(A, B, 0.5) != None:
        r = quadratic_congruence(A, B, 0.5)
        Q = (r**2 - A)//B
        div = divisors(Q)
        
        B_0 = 0

        for i in div:
            if isinstance(sqrt(Q // i), Integer):
                B_0, d = i, sqrt(Q // i)
                break

        if B_0 > 0:
            X, Y, W = descent(A, B_0)
            return (r*X + W, B*d*Y, A*X + r*W)

   
def quadratic_congruence(a, m, k=1.0):
    """
    Solves the quadratic congruence $x^2 \equiv a \ (mod \ m)$. Returns the
    first solution in the range $0 .. \lfloor k*m \rfloor$.
    Return None if solutions do not exist. Currently uses bruteforce.
    Good enough for ``m`` sufficiently small.

    TODO: An efficient algorithm should be implemented.
    """
    upper_bound = floor(abs(k*m)) + 1
    
    for i in range(upper_bound):
        if (i**2 - a)% abs(m) == 0:
            return i
    
    return None