def sign(v):
    if v < 0 :
        return -1
    if v > 0 :
        return 1
    return 0

def sqrt_base_phi(k, l, m, n, p = 100):  
    if sign(5*l**2*sign(l)+sign(2*k+l)*(2*k+l)**2) <= 0:
        raise Exception("numerator k + l * phi not positive")
    if sign(5*n**2*sign(n)+sign(2*m+n)*(2*m+n)**2) <= 0:
        raise Exception("denominator m + n * phi not positive")
    d=[]
    e, a, b, x, y = -1, -1, 1, 0, 0
    i, j = 2*(x+a)+(y+b), (y+b)
    u, v = 2*i*j*(2*m+n)+n*(i*i+5*j*j)-4*l, 8*k-10*n*i*j+4*l-(2*m+n)*(i*i+5*j*j)
    while 5*sign(u)*u**2 <= sign(v)*v**2: #find highest place value
        e, a, b, p = e+1, b, b+a, p+1
        i, j = 2*(x+a)+(y+b), (y+b)
        u, v = 2*i*j*(2*m+n)+n*(i*i+5*j*j)-4*l, 8*k-10*n*i*j+4*l-(2*m+n)*(i*i+5*j*j)
    while p > 0:
        i, j = 2*(x+a)+(y+b), (y+b)
        u, v = 2*i*j*(2*m+n)+n*(i*i+5*j*j)-4*l, 8*k-10*n*i*j+4*l-(2*m+n)*(i*i+5*j*j)
        if 5*sign(u)*u**2 <= sign(v)*v**2: #((x+a)+(y+b)*phi)** 2 <= (k+l*phi)/(m+n*phi)
            d.append([e, 1]) 
            x, y = x+a,y+b
        elif d:
            d.append([e, 0])
        e, a, b, p = e-1, b-a, a, p-1
    return d # [[exponent, digit]]

print(sqrt_base_phi(0, 1, 1, 0, 100))     # A352594 - Expansion of sqrt(phi) in base phi.  
#print(sqrt_base_phi(16, 0, 0, 1, 100))   # A202108 - Expansion of 4/sqrt(phi) in base phi. 
#print(sqrt_base_phi(9, 0, 4, 0, 100))    # A173857 - Expansion of 3/2 in base phi.
#print(sqrt_base_phi(100, 0, 81, 0, 100)) # A173856 - Expansion of 10/9 in base phi.
#print(sqrt_base_phi(16, 0, 9, 0, 100))   # A173858 - Expansion of 4/3 in base phi.
#print(sqrt_base_phi(k, l, m, n, p))      # *USAGE* - Expansion of (k+l*phi)/(m+n*phi) in base phi (upto p fractional places)