import numpy as np
def a( d1, d2, x1, y1, x2, y2 ):
    numerator=d1**2-d2**2 - ((x1**2+y1**2)-(x2**2+y2**2))
    denominator=2*(y2-y1)
    return numerator/denominator

def b( x1, y1, x2, y2 ):
    return -(x2-x1)/(y2-y1)

def solve_xy( x1, y1, x2, y2, d1, d2 ):
    bb=b(x1,y1,x2,y2)
    aa=a(d1,d2,x1,y1,x2,y2)
    rad=4*(bb*(aa-y1)-x1)**2 - 4*(1+bb**2)*(x1**2-d1**2+(y1-aa)**2)
    pre=2*(x1-bb*(aa-y1))
    den=2*(1+bb**2)
    xp=(pre+np.sqrt(rad))/den
    xm=(pre-np.sqrt(rad))/den
    yp=aa+xp*bb
    ym=aa+xm*bb
    return xm,ym,xp,yp