Calculating the point of intersection of two parabolas

I recently wrote a function to calculate the intersection point(s) of two parabolas. Working it out on paper was a nightmare, but I eventually did it (I think it's correct), and ended up with:

\$ x = \frac{\sqrt{c_2 - c_1 + \frac{b_1^2}{4a_1} - \frac{b_2^2}{4a_2} + E}}{\sqrt{a_2 - a_1}} - \frac{b1 - b2}{2(a1 - a2)} \$

with

\$ E = \frac{b_1^2a_2^3(a_2 - 2a_1^2) - a_1^3b_2^2(a_1 - 2a_2^2) - 2a_1^2a_2^2b_2(b_2 - b_1)}{(4a_1^2a_2^2)(a_1 - a_2)^2} \$

For any two parabolas \$ a_1x^2 + b_1x + c_1 = 0\$ and \$a_2x^2 + b_2x + c_2 = 0\$.

That can't be the fastest way to do it, I must be missing something. Does anyone know of any better ways I could implement parabola to parabola intersection point calculations?

The function I wrote is below:

def parabola_to_parabola_poi(a1, b1, c1, a2, b2, c2):
    """
    Calculate the intersection point(s) of two parabolas.
    """
    a1_sqrd, a2_sqrd = a1*a1, a2*a2
    b1_sqrd, b2_sqrd = b1*b1, b2*b2
    c1_sqrd, c2_sqrd = c1*c1, c2*c2

    E1 = b1_sqrd*a2_sqrd*a2*(a2 - 2*a1_sqrd)
    E2 = a1_sqrd*a1*b2_sqrd*(a1 - 2*a2_sqrd)
    E3 = 2*a1_sqrd*a2_sqrd*b2*(b2 - b1)
    E4 = (4*a1_sqrd*a2_sqrd)*(a1 - a2)*(a1 - a2)

    E = (E1 - E2 - E3)/E4

    F = (c2 - c1 + b1_sqrd/(4*a1) - b2_sqrd/(4*a2) + E)/(a2 - a1)
    G = (b1 - b2)/(2*(a1 - a2))

    if not F:
        px = math.sqrt(F) - G
        py = a1*px*px + b1*px + c1
        return [
            (px, py)
            ]
    elif F < 0:
        return []
    sqrt_F = math.sqrt(F)

    px1 = sqrt_F - G
    px2 = -sqrt_F - G

    py1 = a1*px1*px1 + b1*px1 + c1
    py2 = a1*px2*px2 + b1*px2 + c1

    return [
        (px1, py1),
        (px2, py2)
        ]

Well it seems as though I SIGNIFICANTLY overcomplicated the equation for some reason. I just did the equation it a second time and I got something much nicer:

For any two parabolas \$ ax^2 + bx + c = 0\$ and \$dx^2 + ex + f = 0\$, we get:

\$
ax^2 - dx^2 + bx - ex + c - f = 0 \\
x^2(a - d) + x(b - e) = f - c
\$

Now it's just a matter of completing the square.

\$
x^2(a - d) + x(b - e) + \frac{(b - e)^2}{4(a - d)} = f - c + \frac{(b - e)^2}{4(a - d)} \\
(x\sqrt{a - d} + \frac{b - e}{2\sqrt{a - d}})^2 = f - c + \frac{(b - e)^2}{4(a - d)} \\
(a - d)(x + \frac{b - e}{2(a - d)})^2 = f - c + \frac{(b - e)^2}{4(a - d)} \\
x + \frac{b - e}{2(a - d)} = \sqrt{\frac{f - c + \frac{(b - e)^2}{a - d}}{a - d}} \\
x = \sqrt{\frac{f - c + \frac{(b - e)^2}{a - d}}{a - d}} - \frac{b - e}{2(a - d)} \\
\$

I don't know how I ended up with the monster I originally posted, but this seems much more reasonable.

StackExchange Code Review Q#51011, answer score: 4