How to solve this algebraic equation by SageMath (rather than by hand)

edit

sage: u,r,t,a = var('u','r','t','a') 
sage: g = (sqrt(a^2 + (r+t)^2) - sqrt(a^2 + (r-t)^2))/(2*r)
sage: E = u == g ; E
u == 1/2*(sqrt(a^2 + (r + t)^2) - sqrt(a^2 + (r - t)^2))/r

This can't be solved "automagically" by Sage's default solver (i. e. Maxima's).

However, there are some worarounds. First, one can "pilot" Sage's default solver to the solution. First, we try to eliminate the radicals by squaring :

sage: E2 = (E^2).expand(); E2
u^2 == 1/2*a^2/r^2 + 1/2*t^2/r^2 - 1/2*sqrt(a^2 + r^2 + 2*r*t + t^2)*sqrt(a^2 + r^2 - 2*r*t + t^2)/r^2 + 1/2

isolating the radicals :

sage: RP=E2.rhs().operands()[2] ; RP
-1/2*sqrt(a^2 + r^2 + 2*r*t + t^2)*sqrt(a^2 + r^2 - 2*r*t + t^2)/r^2
sage: PP=E2.rhs()-RP
sage: E2-PP
u^2 - 1/2*a^2/r^2 - 1/2*t^2/r^2 - 1/2 == -1/2*sqrt(a^2 + r^2 + 2*r*t + t^2)*sqrt(a^2 + r^2 - 2*r*t + t^2)/r^2

and squaring again :

sage: E4=((E2-PP)^2).expand(); E4
u^4 - a^2*u^2/r^2 - t^2*u^2/r^2 - u^2 + 1/4*a^4/r^4 + 1/2*a^2/r^2 + 1/2*a^2*t^2/r^4 + 1/2*t^2/r^2 + 1/4*t^4/r^4 + 1/4 == 1/4*a^4/r^4 + 1/2*a^2/r^2 + 1/2*a^2*t^2/r^4 - 1/2*t^2/r^2 + 1/4*t^4/r^4 + 1/4

Any solution of E is solution of E4, which is now a polynomial equation in t ultimately of the form A*t^2+B==0. In turns out that Maxima can solve it:

sage: E4.solve(t)
[t == -sqrt(r^2*u^2/(u^2 - 1) - a^2/(u^2 - 1) - r^2/(u^2 - 1))*u, t == sqrt(r^2*u^2/(u^2 - 1) - a^2/(u^2 - 1) - r^2/(u^2 - 1))*u]

$$\left[t = -\sqrt{\frac{r^{2} u^{2}}{u^{2} - 1} - \frac{a^{2}}{u^{2} - 1} - \frac{r^{2}}{u^{2} - 1}} u, t = \sqrt{\frac{r^{2} u^{2}}{u^{2} - 1} - \frac{a^{2}}{u^{2} - 1} - \frac{r^{2}}{u^{2} - 1}} u\right]$$

Since we have squared our equations (twice..), we may have introduced "spurious solutions. The candidate solutions above must be checked. This is not trivial, and is therefore cowardly left to the reader "as an exercise".

It turns out that sympy can solve the original equation E, but that Sage is currently unable to translate it back to Sage: sage: SS=E.solve(t, algorithm="sympy"); SS

ConditionSet(t, Eq(-2*r*u - sqrt(a**2 + (r - t)**2) + sqrt(a**2 + (r + t)**2), 0), {-u*sqrt((-a**2 + r**2*u**2 - r**2)/((u - 1)*(u + 1))), u*sqrt((-a**2 + r**2*u**2 - r**2)/((u - 1)*(u + 1)))})
sage: SS._sage_()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-180-8b468472da8b> in <module>()
----> 1 SS._sage_()

AttributeError: 'ConditionSet' object has no attribute '_sage_'

One can manually extract the roots proposed by sympy:

sage: [u._sage_() for u in SS.args[2].args]
[u*sqrt((r^2*u^2 - a^2 - r^2)/((u + 1)*(u - 1))),
 -u*sqrt((r^2*u^2 - a^2 - r^2)/((u + 1)*(u - 1)))]

which turn out to be the same as those proposed by Maxima:

$$\left[u \sqrt{\frac{r^{2} u^{2} - a^{2} - r^{2}}{{\left(u + 1\right)} {\left(u - 1\right)}}}, -u \sqrt{\frac{r^{2} u^{2} - a^{2} - r^{2}}{{\left(u + 1\right)} {\left(u - 1\right)}}}\right]$$

EDIT : One can also directly call sympy and sk for a dictionary solution simpler to translate back in Sage:

sage: import sympy
sage: [{t._sage_():u[t]._sage_() for t in u.keys()} for u in sympy.solve(sympy.sympify(E),
....:  t, dict=True)]
[{t: u*sqrt((r^2*u^2 - a^2 - r^2)/(u^2 - 1))},
 {t: -u*sqrt((r^2*u^2 - a^2 - r^2)/(u^2 - 1))}]

EDIT 2: Along the same lines, one can use FriCAS to solve the same equation (with contorsions to to convert the result, again...):

sage: SG=fricas.solve(*[fricas(v) for v in [g-u, t]])
sage: [SG[w].lhs().sage()==SG[w].rhs().sage() for w in range(len(SG))]
[t == u*sqrt((r^2*u^2 - a^2 - r^2)/(u^2 - 1)),
 t == -u*sqrt((r^2*u^2 - a^2 - r^2)/(u^2 - 1))]

EDIT 3: Ditto for giac, minus the contorsions:

sage: giac.solve(giac(u-g), giac(t)).sage()
[sqrt(r^2*u^4 - a^2*u^2 - 2*r^2*u^2 + a^2 + r^2)*u/(u^2 - 1),
 -sqrt(r^2*u^4 - a^2*u^2 - 2*r^2*u^2 + a^2 + r^2)*u/(u^2 - 1)]

Mathematica gives the same solutions (again not automagically translatable to Sage):

sage: mathematica.Solve(E,t)
{{t -> -(Sqrt[-(a^2*u^2) - r^2*u^2 + r^2*u^4]/Sqrt[-1 + u^2])}, 
 {t -> Sqrt[-(a^2*u^2) - r^2*u^2 + r^2*u^4]/Sqrt[-1 + u^2]}}

HTH,