Dealing with negative roots in code

Question

I have an equation in the following form that needs to be solved for $x$, with $x > 0$ and constants $a,b,c,d \geq 0$.

$\left(\frac{x - a}{b}\right)^2 + \left(\frac{x - c}{d}\right)^2 = 1$

From the following paper, which i am trying to implement: https://grail.cs.washington.edu/projects/crowd-flows/78-treuille.pdf

Wolfram Alpha gives the following solution:

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

http://www.wolframalpha.com/input/?i=((x+-+a)+%2F+b)+%5E+2+%2B+((x+-+c)+%2F+d)+%5E+2++%3D+1+solve+for+x

However this can result in the root of a negative number, and resulting in a complex number, which is not what the paper seems to want. I am unsure on how to solve this (in code) to get a real value for x.

So my question is, how do i solve this equation as desired by the paper?

---

PS. Limiting to real x still has the same problem

http://www.wolframalpha.com/input/?i=((x+-+a)+%2F+b)+%5E+2+%2B+((x+-+c)+%2F+d)+%5E+2++%3D+1+solve+for+real+x

Answer 1

Consider the equation $\left(\frac{x-a}{b}\right)^2 + \left(\frac{y-c}{d}\right)^2 = 1$. It is an ellipsis centered at $(a,c)$ with axes parallel to origin of lengths $b$ and $d$ respectively. The solution of your task are the points where this ellipsis intersects with the line $y=x$. There may be no intersection - and that is exactly the case of complex roots.

Answer 2

It's $$d^2(x^2-2ax+a^2)+b^2(x^2-2cx+c^2)=b^2d^2$$ or $$(b^2+d^2)x^2-2(d^2a+b^2c)x+a^2d^2+b^2c^2-b^2d^2=0,$$ which gives $$x=\frac{d^2a+b^2c\pm\sqrt{(d^2a+b^2c)^2-(b^2+d^2)(a^2d^2+b^2c^2-b^2d^2)}}{b^2+d^2}$$ or $$x=\frac{d^2a+b^2c\pm\sqrt{d^4a^2+2b^2d^2ac+b^4c^2-b^2a^2d^2-b^4c^2+b^4d^2-d^4a^2-b^2c^2d^2+d^4b^2}}{b^2+d^2}$$ or $$x=\frac{d^2a+b^2c\pm bd\sqrt{2ac-a^2-c^2+b^2+d^2}}{b^2+d^2}.$$ We need $b^2+d^2\neq0$ and $2ac-a^2-c^2+b^2+d^2\geq0$.

Thus, for the biggest root we obtain: $$x=\frac{d^2a+b^2c+bd\sqrt{2ac-a^2-c^2+b^2+d^2}}{b^2+d^2}.$$