Minimum eccentricity of ellipses around another ellipse

exact_code.sage

# Exact computation to solve http://math.stackexchange.com/q/688861/35416,
# heavily inspired by https://groups.google.com/d/msg/sci.math/LYtIdRhk2ac/mQtEBcgCjZkJ

def printPoly(name, poly):
    print(name + ":")
    print(poly.denominator()*poly)
    print("")

# a and b are the alpha and beta from http://math.stackexchange.com/a/698656/35416.
# We use approximate values to choose the right alternatives in some situations.
approx = dict(a = 0.147739677661122668796455680683973472705,
              b = 0.779540453991525714639344299875163938705)

PR.<a,b,x,y> = QQ[]                                                      # polynomial ring in 4 variables
PR2.<mu> = PR[]                                                          # another ring on top of that
E1 = diagonal_matrix(PR, [a, 1, -1])                                     # center ellipse
RotScale = matrix(PR, [[1, b, 0], [-b, 1, 0], [0, 0, 1]])                # rotation with scaling
# We undo the scaling using the diagonal marix in the following step:
E2 = diagonal_matrix(PR, [1, 1, 1+b^2])*RotScale.transpose()*E1*RotScale # rotated
Transl = matrix(PR, [[1, 0, -x], [0, 1, -y], [0, 0, 1]])                 # translation
E3 = Transl.transpose()*E2*Transl                                        # rotated and translated ellipse
print("E1: " + latex(E1))
print("E3: " + latex(E3))
print("")

# Start with conditions that E3 touches the axes of the coordinate system:
eqs = [(v.row()*E3.adjoint()*v.column())[0, 0] for v in [vector([1,0,0]), vector([0,1,0])]]
# Add condition that the ellipses touch one another:
eqs.append(det(mu*E1 + E3).discriminant())
print("Equations:")
for eq in eqs:
    print(eq)
print("")

# Eliminate x and y (i.e. translation), resulting in one polynomial in a and b:
abPoly = eqs[2].resultant(eqs[0], x).resultant(eqs[1], y)
print("abPoly has {} terms".format(len(list(abPoly))))
abFactors = abPoly.change_ring(QQ).factor()
abFactor = sorted(abFactors, key = lambda f: abs(f[0](**approx)))[0][0]  # the relevant factor is the last one
printPoly("abFactor", abFactor)
aPoly = abFactor.polynomial(b).discriminant().univariate_polynomial()    # discriminant wrt. b
print("aPoly has {} terms".format(len(list(aPoly))))
aFactors = aPoly.factor()
aFactor = sorted(aFactors, key = lambda f: abs(f[0](**approx)))[0][0]    # relevant factor is this one
printPoly("aFactor", aFactor)
aRoots = aFactor.roots(AA, False)
a0 = sorted(aRoots, key = lambda a: abs(float(a) - approx['a']))[0]      # find the one we seek
printPoly("alpha minpoly", a0.minpoly())
bRoots = abFactor(a=a0).univariate_polynomial().roots(AA, False)
b0 = max(bRoots)
printPoly("beta minpoly", b0.minpoly())
x0 = max(eqs[0](a=a0, b=b0).univariate_polynomial().roots(AA, False))
y0 = max(eqs[1](a=a0, b=b0).univariate_polynomial().roots(AA, False))

# Now change parameters from our specific ones to more common names:
a1 = AA(1)                                                               # semi-major axis
b1 = sqrt(a0)                                                            # semi-minor axis
c1 = sqrt(1 - a0)                                                        # numerical eccentricity
x1 = x0*b1                                                               # horizontal translation
y1 = y0*b1                                                               # vertical translation

# So far everything is exact, but angles will not be exact any more.
RF = RealField(512)                                                      # VERY precise approximations
phiRad = atan(RF(b0))                                                    # rotation angle in radians
phiDeg = RF(phiRad*180/pi)                                               # likewise in degrees

# Print results, with lots of digits:
for varname in ["b1", "c1", "phiRad", "phiDeg", "x1", "y1", "a0", "b0"]:
    print("{}: {}".format(varname, locals()[varname].n(digits=60)))

exact_output.txt

E1: \left(\begin{array}{rrr}
a & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & -1
\end{array}\right)
E3: \left(\begin{array}{rrr}
b^{2} + a & a b -  b & - b^{2} x -  a b y -  a x + b y \\
a b -  b & a b^{2} + 1 & - a b^{2} y -  a b x + b x -  y \\
- b^{2} x -  a b y -  a x + b y & - a b^{2} y -  a b x + b x -  y & a b^{2} y^{2} + b^{2} x^{2} + 2 a b x y + a x^{2} - 2 b x y -  b^{2} + y^{2} - 1
\end{array}\right)

Equations:
a*b^4*x^2 - a*b^4 + 2*a*b^2*x^2 - a*b^2 + a*x^2 - b^2 - 1
a*b^4*y^2 + 2*a*b^2*y^2 - b^4 - a*b^2 + a*y^2 - b^2 - a
a^8*b^12*x^4*y^4 + 2*a^7*b^12*x^6*y^2 + 4*a^8*b^11*x^5*y^3 + 2*a^7*b^12*x^2*y^6 + a^6*b^12*x^8 + 4*a^7*b^11*x^7*y - 2*a^8*b^12*x^4*y^2 + 6*a^8*b^10*x^6*y^2 - 4*a^7*b^11*x^5*y^3 - 2*a^8*b^12*x^2*y^4 + 4*a^8*b^10*x^4*y^4 + 4*a^6*b^12*x^4*y^4 + 8*a^7*b^11*x^3*y^5 + a^6*b^12*y^8 + 2*a^7*b^12*x^6 + 2*a^7*b^10*x^8 - 4*a^8*b^11*x^5*y + 4*a^8*b^9*x^7*y - 4*a^6*b^11*x^7*y - 6*a^7*b^12*x^4*y^2 + 2*a^5*b^12*x^6*y^2 - 8*a^8*b^11*x^3*y^3 + 16*a^8*b^9*x^5*y^3 + 8*a^6*b^11*x^5*y^3 - 6*a^7*b^12*x^2*y^4 + 14*a^7*b^10*x^4*y^4 - 8*a^6*b^11*x^3*y^5 + 2*a^7*b^12*y^6 + 8*a^7*b^10*x^2*y^6 + 2*a^5*b^12*x^2*y^6 + 4*a^6*b^11*x*y^7 + a^8*b^12*x^4 - 2*a^8*b^10*x^6 - 4*a^6*b^12*x^6 + a^8*b^8*x^8 + 4*a^6*b^10*x^8 - 8*a^7*b^11*x^5*y + 12*a^7*b^9*x^7*y + 4*a^8*b^12*x^2*y^2 - 20*a^8*b^10*x^4*y^2 + 2*a^6*b^12*x^4*y^2 + 24*a^8*b^8*x^6*y^2 + 10*a^6*b^10*x^6*y^2 - 8*a^7*b^11*x^3*y^3 - 4*a^7*b^9*x^5*y^3 - 8*a^5*b^11*x^5*y^3 + a^8*b^12*y^4 - 6*a^8*b^10*x^2*y^4 + 2*a^6*b^12*x^2*y^4 + 6*a^8*b^8*x^4*y^4 + a^4*b^12*x^4*y^4 + 16*a^7*b^11*x*y^5 + 32*a^7*b^9*x^3*y^5 + 4*a^5*b^11*x^3*y^5 - 4*a^6*b^12*y^6 + 10*a^6*b^10*x^2*y^6 - 4*a^5*b^11*x*y^7 + 4*a^6*b^10*y^8 - 6*a^7*b^12*x^4 + 10*a^7*b^10*x^6 - 4*a^5*b^12*x^6 + 8*a^7*b^8*x^8 + 8*a^8*b^11*x^3*y - 24*a^8*b^9*x^5*y + 16*a^8*b^7*x^7*y - 16*a^6*b^9*x^7*y + 12*a^7*b^12*x^2*y^2 - 28*a^7*b^10*x^4*y^2 + 2*a^5*b^12*x^4*y^2 - 16*a^7*b^8*x^6*y^2 + 8*a^5*b^10*x^6*y^2 + 4*a^8*b^11*x*y^3 - 24*a^8*b^9*x^3*y^3 + 40*a^6*b^11*x^3*y^3 + 24*a^8*b^7*x^5*y^3 + 20*a^6*b^9*x^5*y^3 - 6*a^7*b^12*y^4 + 16*a^7*b^10*x^2*y^4 + 2*a^5*b^12*x^2*y^4 + 56*a^7*b^8*x^4*y^4 + 14*a^5*b^10*x^4*y^4 - 28*a^6*b^11*x*y^5 - 20*a^6*b^9*x^3*y^5 - 4*a^4*b^11*x^3*y^5 + 6*a^7*b^10*y^6 - 4*a^5*b^12*y^6 + 12*a^7*b^8*x^2*y^6 + 16*a^6*b^9*x*y^7 + 2*a^5*b^10*y^8 - 2*a^8*b^12*x^2 + 8*a^8*b^10*x^4 - 10*a^8*b^8*x^6 - 26*a^6*b^10*x^6 + 2*a^4*b^12*x^6 + 4*a^8*b^6*x^8 + 6*a^6*b^8*x^8 + 12*a^7*b^11*x^3*y - 36*a^7*b^9*x^5*y + 28*a^5*b^11*x^5*y + 8*a^7*b^7*x^7*y - 2*a^8*b^12*y^2 + 18*a^8*b^10*x^2*y^2 - 6*a^6*b^12*x^2*y^2 - 48*a^8*b^8*x^4*y^2 + 78*a^6*b^10*x^4*y^2 - 6*a^4*b^12*x^4*y^2 + 36*a^8*b^6*x^6*y^2 + 40*a^6*b^8*x^6*y^2 - 36*a^7*b^11*x*y^3 + 8*a^7*b^9*x^3*y^3 - 40*a^5*b^11*x^3*y^3 + 24*a^7*b^7*x^5*y^3 - 32*a^5*b^9*x^5*y^3 + 2*a^8*b^10*y^4 - 6*a^8*b^8*x^2*y^4 - 112*a^6*b^10*x^2*y^4 - 6*a^4*b^12*x^2*y^4 + 4*a^8*b^6*x^4*y^4 - 34*a^6*b^8*x^4*y^4 + 4*a^4*b^10*x^4*y^4 + 48*a^7*b^9*x*y^5 + 48*a^7*b^7*x^3*y^5 + 4*a^5*b^9*x^3*y^5 - 12*a^6*b^10*y^6 + 2*a^4*b^12*y^6 + 40*a^6*b^8*x^2*y^6 + 6*a^4*b^10*x^2*y^6 - 12*a^5*b^9*x*y^7 + 6*a^6*b^8*y^8 + 4*a^7*b^12*x^2 - 24*a^7*b^10*x^4 + 10*a^5*b^12*x^4 + 20*a^7*b^8*x^6 - 12*a^5*b^10*x^6 + 12*a^7*b^6*x^8 - 4*a^8*b^11*x*y + 28*a^8*b^9*x^3*y - 48*a^8*b^7*x^5*y + 36*a^6*b^9*x^5*y - 16*a^4*b^11*x^5*y + 24*a^8*b^5*x^7*y - 24*a^6*b^7*x^7*y + 4*a^7*b^12*y^2 - 20*a^5*b^12*x^2*y^2 - 42*a^7*b^8*x^4*y^2 - 112*a^5*b^10*x^4*y^2 - 2*a^3*b^12*x^4*y^2 - 24*a^7*b^6*x^6*y^2 + 12*a^5*b^8*x^6*y^2 + 8*a^8*b^9*x*y^3 + 12*a^6*b^11*x*y^3 - 24*a^8*b^7*x^3*y^3 - 24*a^6*b^9*x^3*y^3 + 8*a^4*b^11*x^3*y^3 + 16*a^8*b^5*x^5*y^3 - 12*a^7*b^10*y^4 + 10*a^5*b^12*y^4 + 84*a^7*b^8*x^2*y^4 + 78*a^5*b^10*x^2*y^4 - 2*a^3*b^12*x^2*y^4 + 84*a^7*b^6*x^4*y^4 + 56*a^5*b^8*x^4*y^4 - 72*a^6*b^9*x*y^5 + 8*a^4*b^11*x*y^5 - 16*a^4*b^9*x^3*y^5 + 6*a^7*b^8*y^6 - 26*a^5*b^10*y^6 + 8*a^7*b^6*x^2*y^6 - 16*a^5*b^8*x^2*y^6 + 24*a^6*b^7*x*y^7 - 4*a^4*b^9*x*y^7 + 8*a^5*b^8*y^8 + a^8*b^12 - 8*a^8*b^10*x^2 + 19*a^8*b^8*x^4 - 18*a^6*b^10*x^4 - 18*a^8*b^6*x^6 - 64*a^6*b^8*x^6 + 6*a^4*b^10*x^6 + 6*a^8*b^4*x^8 + 4*a^6*b^6*x^8 + 16*a^7*b^11*x*y + 36*a^7*b^9*x^3*y - 40*a^5*b^11*x^3*y - 72*a^7*b^7*x^5*y + 72*a^5*b^9*x^5*y - 8*a^7*b^5*x^7*y - 4*a^8*b^10*y^2 + 24*a^8*b^8*x^2*y^2 + 66*a^6*b^10*x^2*y^2 - 6*a^4*b^12*x^2*y^2 - 44*a^8*b^6*x^4*y^2 + 180*a^6*b^8*x^4*y^2 + 16*a^4*b^10*x^4*y^2 + 24*a^8*b^4*x^6*y^2 + 60*a^6*b^6*x^6*y^2 - 84*a^7*b^9*x*y^3 + 40*a^5*b^11*x*y^3 + 72*a^7*b^7*x^3*y^3 + 24*a^5*b^9*x^3*y^3 + 8*a^3*b^11*x^3*y^3 + 56*a^7*b^5*x^5*y^3 - 48*a^5*b^7*x^5*y^3 + a^8*b^8*y^4 - 24*a^6*b^10*y^4 - 2*a^8*b^6*x^2*y^4 - 348*a^6*b^8*x^2*y^4 - 28*a^4*b^10*x^2*y^4 + a^8*b^4*x^4*y^4 - 56*a^6*b^6*x^4*y^4 + 6*a^4*b^8*x^4*y^4 + 48*a^7*b^7*x*y^5 - 36*a^5*b^9*x*y^5 + 4*a^3*b^11*x*y^5 + 32*a^7*b^5*x^3*y^5 - 24*a^5*b^7*x^3*y^5 - 12*a^6*b^8*y^6 + 10*a^4*b^10*y^6 + 60*a^6*b^6*x^2*y^6 + 24*a^4*b^8*x^2*y^6 - 8*a^5*b^7*x*y^7 + 4*a^6*b^6*y^8 + a^4*b^8*y^8 - 4*a^7*b^12 + 14*a^7*b^10*x^2 - 2*a^5*b^12*x^2 - 30*a^7*b^8*x^4 + 68*a^5*b^10*x^4 - 6*a^3*b^12*x^4 + 16*a^7*b^6*x^6 - 12*a^5*b^8*x^6 + 8*a^7*b^4*x^8 - 8*a^8*b^9*x*y - 24*a^6*b^11*x*y + 32*a^8*b^7*x^3*y + 96*a^6*b^9*x^3*y - 12*a^4*b^11*x^3*y - 40*a^8*b^5*x^5*y + 120*a^6*b^7*x^5*y - 48*a^4*b^9*x^5*y + 16*a^8*b^3*x^7*y - 16*a^6*b^5*x^7*y - 2*a^7*b^10*y^2 - 2*a^5*b^12*y^2 - 36*a^7*b^8*x^2*y^2 - 168*a^5*b^10*x^2*y^2 + 12*a^3*b^12*x^2*y^2 - 18*a^7*b^6*x^4*y^2 - 348*a^5*b^8*x^4*y^2 - 6*a^3*b^10*x^4*y^2 - 6*a^7*b^4*x^6*y^2 + 8*a^5*b^6*x^6*y^2 + 4*a^8*b^7*x*y^3 - 108*a^6*b^9*x*y^3 - 8*a^8*b^5*x^3*y^3 - 336*a^6*b^7*x^3*y^3 - 8*a^4*b^9*x^3*y^3 + 4*a^8*b^3*x^5*y^3 - 40*a^6*b^5*x^5*y^3 - 6*a^7*b^8*y^4 + 68*a^5*b^10*y^4 - 6*a^3*b^12*y^4 + 96*a^7*b^6*x^2*y^4 + 180*a^5*b^8*x^2*y^4 - 20*a^3*b^10*x^2*y^4 + 56*a^7*b^4*x^4*y^4 + 84*a^5*b^6*x^4*y^4 - 48*a^6*b^7*x*y^5 + 36*a^4*b^9*x*y^5 + 40*a^6*b^5*x^3*y^5 - 24*a^4*b^7*x^3*y^5 + 2*a^7*b^6*y^6 - 64*a^5*b^8*y^6 - 2*a^3*b^10*y^6 + 2*a^7*b^4*x^2*y^6 - 24*a^5*b^6*x^2*y^6 + 16*a^6*b^5*x*y^7 - 16*a^4*b^7*x*y^7 + 12*a^5*b^6*y^8 + 2*a^8*b^10 + 4*a^6*b^12 - 10*a^8*b^8*x^2 + 6*a^6*b^10*x^2 - 2*a^4*b^12*x^2 + 18*a^8*b^6*x^4 - 63*a^6*b^8*x^4 - 24*a^4*b^10*x^4 + a^2*b^12*x^4 - 14*a^8*b^4*x^6 - 76*a^6*b^6*x^6 + 6*a^4*b^8*x^6 + 4*a^8*b^2*x^8 + a^6*b^4*x^8 + 20*a^7*b^9*x*y + 20*a^5*b^11*x*y + 60*a^7*b^7*x^3*y - 344*a^5*b^9*x^3*y + 36*a^3*b^11*x^3*y - 80*a^7*b^5*x^5*y + 48*a^5*b^7*x^5*y - 12*a^7*b^3*x^7*y - 2*a^8*b^8*y^2 + 30*a^6*b^10*y^2 - 2*a^4*b^12*y^2 + 10*a^8*b^6*x^2*y^2 + 60*a^6*b^8*x^2*y^2 + 66*a^4*b^10*x^2*y^2 + 4*a^2*b^12*x^2*y^2 - 14*a^8*b^4*x^4*y^2 + 80*a^6*b^6*x^4*y^2 + 84*a^4*b^8*x^4*y^2 + 6*a^8*b^2*x^6*y^2 + 40*a^6*b^4*x^6*y^2 - 60*a^7*b^7*x*y^3 + 344*a^5*b^9*x*y^3 - 12*a^3*b^11*x*y^3 + 88*a^7*b^5*x^3*y^3 + 336*a^5*b^7*x^3*y^3 + 24*a^3*b^9*x^3*y^3 + 44*a^7*b^3*x^5*y^3 - 32*a^5*b^5*x^5*y^3 - 57*a^6*b^8*y^4 - 18*a^4*b^10*y^4 + a^2*b^12*y^4 - 352*a^6*b^6*x^2*y^4 - 42*a^4*b^8*x^2*y^4 - 24*a^6*b^4*x^4*y^4 + 4*a^4*b^6*x^4*y^4 + 16*a^7*b^5*x*y^5 - 120*a^5*b^7*x*y^5 + 24*a^3*b^9*x*y^5 + 8*a^7*b^3*x^3*y^5 - 56*a^5*b^5*x^3*y^5 - 4*a^6*b^6*y^6 + 20*a^4*b^8*y^6 + 40*a^6*b^4*x^2*y^6 + 36*a^4*b^6*x^2*y^6 + 8*a^5*b^5*x*y^7 + a^6*b^4*y^8 + 4*a^4*b^6*y^8 - 4*a^7*b^10 + 4*a^5*b^12 + 4*a^7*b^8*x^2 - 8*a^5*b^10*x^2 + 136*a^5*b^8*x^4 - 12*a^3*b^10*x^4 - 2*a^7*b^4*x^6 - 4*a^5*b^6*x^6 + 2*a^7*b^2*x^8 - 4*a^8*b^7*x*y + 12*a^6*b^9*x*y - 20*a^4*b^11*x*y + 12*a^8*b^5*x^3*y + 228*a^6*b^7*x^3*y + 108*a^4*b^9*x^3*y - 4*a^2*b^11*x^3*y - 12*a^8*b^3*x^5*y + 144*a^6*b^5*x^5*y - 48*a^4*b^7*x^5*y + 4*a^8*b*x^7*y - 4*a^6*b^3*x^7*y - 16*a^7*b^8*y^2 - 28*a^5*b^10*y^2 - 24*a^7*b^6*x^2*y^2 - 96*a^5*b^8*x^2*y^2 + 8*a^7*b^4*x^4*y^2 - 352*a^5*b^6*x^4*y^2 - 6*a^3*b^8*x^4*y^2 + 8*a^7*b^2*x^6*y^2 + 2*a^5*b^4*x^6*y^2 - 336*a^6*b^7*x*y^3 - 96*a^4*b^9*x*y^3 - 8*a^2*b^11*x*y^3 - 464*a^6*b^5*x^3*y^3 - 72*a^4*b^7*x^3*y^3 - 40*a^6*b^3*x^5*y^3 + 136*a^5*b^8*y^4 - 24*a^3*b^10*y^4 + 34*a^7*b^4*x^2*y^4 + 80*a^5*b^6*x^2*y^4 - 48*a^3*b^8*x^2*y^4 + 14*a^7*b^2*x^4*y^4 + 56*a^5*b^4*x^4*y^4 + 8*a^6*b^5*x*y^5 + 72*a^4*b^7*x*y^5 + 40*a^6*b^3*x^3*y^5 - 16*a^4*b^5*x^3*y^5 - 76*a^5*b^6*y^6 - 10*a^3*b^8*y^6 - 6*a^5*b^4*x^2*y^6 + 4*a^6*b^3*x*y^7 - 24*a^4*b^5*x*y^7 + 8*a^5*b^4*y^8 + a^8*b^8 - 16*a^6*b^10 - 10*a^4*b^12 - 4*a^8*b^6*x^2 + 60*a^6*b^8*x^2 - 28*a^4*b^10*x^2 + 4*a^2*b^12*x^2 + 6*a^8*b^4*x^4 - 96*a^6*b^6*x^4 - 57*a^4*b^8*x^4 + 2*a^2*b^10*x^4 - 4*a^8*b^2*x^6 - 44*a^6*b^4*x^6 + 2*a^4*b^6*x^6 + a^8*x^8 - 8*a^7*b^7*x*y - 80*a^5*b^9*x*y + 24*a^3*b^11*x*y + 60*a^7*b^5*x^3*y - 712*a^5*b^7*x^3*y + 84*a^3*b^9*x^3*y - 48*a^7*b^3*x^5*y - 8*a^5*b^5*x^5*y - 4*a^7*b*x^7*y + 48*a^6*b^8*y^2 - 8*a^4*b^10*y^2 + 4*a^2*b^12*y^2 - 186*a^6*b^6*x^2*y^2 + 60*a^4*b^8*x^2*y^2 + 18*a^2*b^10*x^2*y^2 - 90*a^6*b^4*x^4*y^2 + 96*a^4*b^6*x^4*y^2 + 10*a^6*b^2*x^6*y^2 - 12*a^7*b^5*x*y^3 + 712*a^5*b^7*x*y^3 - 36*a^3*b^9*x*y^3 + 32*a^7*b^3*x^3*y^3 + 464*a^5*b^5*x^3*y^3 + 24*a^3*b^7*x^3*y^3 + 12*a^7*b*x^5*y^3 - 8*a^5*b^3*x^5*y^3 - 42*a^6*b^6*y^4 - 63*a^4*b^8*y^4 + 8*a^2*b^10*y^4 - 118*a^6*b^4*x^2*y^4 - 18*a^4*b^6*x^2*y^4 + 8*a^6*b^2*x^4*y^4 + a^4*b^4*x^4*y^4 - 144*a^5*b^5*x*y^5 + 48*a^3*b^7*x*y^5 - 44*a^5*b^3*x^3*y^5 + 16*a^4*b^6*y^6 + 10*a^6*b^2*x^2*y^6 + 24*a^4*b^4*x^2*y^6 + 12*a^5*b^3*x*y^7 + 6*a^4*b^4*y^8 + 4*a^7*b^8 + 68*a^5*b^10 + 4*a^3*b^12 - 18*a^7*b^6*x^2 - 82*a^5*b^8*x^2 + 30*a^3*b^10*x^2 - 2*a*b^12*x^2 + 24*a^7*b^4*x^4 + 120*a^5*b^6*x^4 - 6*a^3*b^8*x^4 - 10*a^7*b^2*x^6 + 96*a^6*b^7*x*y + 80*a^4*b^9*x*y - 16*a^2*b^11*x*y + 168*a^6*b^5*x^3*y + 336*a^4*b^7*x^3*y - 8*a^2*b^9*x^3*y + 72*a^6*b^3*x^5*y - 16*a^4*b^5*x^5*y - 10*a^7*b^6*y^2 - 2*a^5*b^8*y^2 + 6*a^3*b^10*y^2 - 2*a*b^12*y^2 + 400*a^5*b^6*x^2*y^2 - 36*a^3*b^8*x^2*y^2 + 6*a^7*b^2*x^4*y^2 - 118*a^5*b^4*x^4*y^2 - 2*a^3*b^6*x^4*y^2 + 4*a^7*x^6*y^2 - 300*a^6*b^5*x*y^3 - 228*a^4*b^7*x*y^3 - 28*a^2*b^9*x*y^3 - 216*a^6*b^3*x^3*y^3 - 88*a^4*b^5*x^3*y^3 - 12*a^6*b*x^5*y^3 + 120*a^5*b^6*y^4 - 30*a^3*b^8*y^4 - 90*a^5*b^4*x^2*y^4 - 44*a^3*b^6*x^2*y^4 + 14*a^5*b^2*x^4*y^4 + 12*a^6*b^3*x*y^5 + 80*a^4*b^5*x*y^5 + 12*a^6*b*x^3*y^5 - 4*a^4*b^3*x^3*y^5 - 44*a^5*b^4*y^6 - 18*a^3*b^6*y^6 + 8*a^5*b^2*x^2*y^6 - 16*a^4*b^3*x*y^7 + 2*a^5*b^2*y^8 - 44*a^6*b^8 - 100*a^4*b^10 + 4*a^2*b^12 + 102*a^6*b^6*x^2 - 2*a^4*b^8*x^2 - 2*a^2*b^10*x^2 - 75*a^6*b^4*x^4 - 42*a^4*b^6*x^4 + a^2*b^8*x^4 - 10*a^6*b^2*x^6 - 12*a^7*b^5*x*y - 220*a^5*b^7*x*y - 12*a^3*b^9*x*y + 4*a*b^11*x*y + 24*a^7*b^3*x^3*y - 552*a^5*b^5*x^3*y + 60*a^3*b^7*x^3*y - 12*a^7*b*x^5*y - 12*a^5*b^3*x^5*y + 6*a^6*b^6*y^2 - 82*a^4*b^8*y^2 + 14*a^2*b^10*y^2 - 258*a^6*b^4*x^2*y^2 - 186*a^4*b^6*x^2*y^2 + 24*a^2*b^8*x^2*y^2 - 78*a^6*b^2*x^4*y^2 + 34*a^4*b^4*x^4*y^2 + 552*a^5*b^5*x*y^3 - 60*a^3*b^7*x*y^3 + 216*a^5*b^3*x^3*y^3 + 8*a^3*b^5*x^3*y^3 - 9*a^6*b^4*y^4 - 96*a^4*b^6*y^4 + 19*a^2*b^8*y^4 + 8*a^4*b^4*x^2*y^4 + 6*a^6*x^4*y^4 - 72*a^5*b^3*x*y^5 + 40*a^3*b^5*x*y^5 - 12*a^5*b*x^3*y^5 - 2*a^4*b^4*y^6 + 6*a^4*b^2*x^2*y^6 + 4*a^5*b*x*y^7 + 4*a^4*b^2*y^8 + 4*a^7*b^6 + 124*a^5*b^8 + 68*a^3*b^10 - 4*a*b^12 - 12*a^7*b^4*x^2 - 148*a^5*b^6*x^2 + 48*a^3*b^8*x^2 - 4*a*b^10*x^2 + 12*a^7*b^2*x^4 + 54*a^5*b^4*x^4 - 4*a^7*x^6 + 60*a^6*b^5*x*y + 220*a^4*b^7*x*y - 20*a^2*b^9*x*y + 36*a^6*b^3*x^3*y + 300*a^4*b^5*x^3*y - 4*a^2*b^7*x^3*y + 12*a^6*b*x^5*y + 72*a^5*b^6*y^2 + 60*a^3*b^8*y^2 - 8*a*b^10*y^2 + 516*a^5*b^4*x^2*y^2 - 24*a^3*b^6*x^2*y^2 - 84*a^6*b^3*x*y^3 - 168*a^4*b^5*x*y^3 - 32*a^2*b^7*x*y^3 - 24*a^6*b*x^3*y^3 - 32*a^4*b^3*x^3*y^3 + 54*a^5*b^4*y^4 - 78*a^5*b^2*x^2*y^4 - 14*a^3*b^4*x^2*y^4 + 48*a^4*b^3*x*y^5 - 10*a^5*b^2*y^6 - 14*a^3*b^4*y^6 + 4*a^5*x^2*y^6 - 4*a^4*b*x*y^7 - 24*a^6*b^6 - 170*a^4*b^8 - 16*a^2*b^10 + b^12 + 48*a^6*b^4*x^2 + 72*a^4*b^6*x^2 - 16*a^2*b^8*x^2 - 24*a^6*b^2*x^4 - 9*a^4*b^4*x^4 - 120*a^5*b^5*x*y - 96*a^3*b^7*x*y + 8*a*b^9*x*y - 144*a^5*b^3*x^3*y + 12*a^3*b^5*x^3*y - 12*a^6*b^4*y^2 - 148*a^4*b^6*y^2 + 4*a^2*b^8*y^2 - 84*a^6*b^2*x^2*y^2 - 258*a^4*b^4*x^2*y^2 + 10*a^2*b^6*x^2*y^2 - 12*a^6*x^4*y^2 + 144*a^5*b^3*x*y^3 - 60*a^3*b^5*x*y^3 + 24*a^5*b*x^3*y^3 - 75*a^4*b^4*y^4 + 18*a^2*b^6*y^4 + 6*a^4*b^2*x^2*y^4 - 12*a^5*b*x*y^5 + 12*a^3*b^3*x*y^5 - 10*a^4*b^2*y^6 + a^4*y^8 + 60*a^5*b^6 + 124*a^3*b^8 - 4*a*b^10 - 72*a^5*b^4*x^2 + 6*a^3*b^6*x^2 - 2*a*b^8*x^2 + 12*a^5*b^2*x^4 + 120*a^4*b^5*x*y + 8*a^2*b^7*x*y + 84*a^4*b^3*x^3*y + 48*a^5*b^4*y^2 + 102*a^3*b^6*y^2 - 10*a*b^8*y^2 + 168*a^5*b^2*x^2*y^2 - 36*a^4*b^3*x*y^3 - 12*a^2*b^5*x*y^3 + 12*a^5*b^2*y^4 + 24*a^3*b^4*y^4 - 12*a^5*x^2*y^4 + 12*a^4*b*x*y^5 - 4*a^3*b^2*y^6 - 80*a^4*b^6 - 44*a^2*b^8 + 2*b^10 + 48*a^4*b^4*x^2 - 10*a^2*b^6*x^2 - 60*a^3*b^5*x*y + 4*a*b^7*x*y - 72*a^4*b^4*y^2 - 18*a^2*b^6*y^2 - 84*a^4*b^2*x^2*y^2 - 24*a^3*b^3*x*y^3 - 24*a^4*b^2*y^4 + 6*a^2*b^4*y^4 - 4*a^4*y^6 + 60*a^3*b^6 + 4*a*b^8 - 12*a^3*b^4*x^2 + 12*a^2*b^5*x*y + 48*a^3*b^4*y^2 - 4*a*b^6*y^2 + 12*a^3*b^2*y^4 - 24*a^2*b^6 + b^8 - 12*a^2*b^4*y^2 + 4*a*b^6

abPoly has 2069 terms
abFactor:
16*a^12*b^18 - 104*a^11*b^18 - 144*a^12*b^16 + 233*a^10*b^18 + 1008*a^11*b^16 - 64*a^9*b^18 + 432*a^12*b^14 - 3179*a^10*b^16 - 340*a^8*b^18 - 3520*a^11*b^14 + 5312*a^9*b^16 - 24*a^7*b^18 - 432*a^12*b^12 + 13194*a^10*b^14 - 84*a^8*b^16 + 822*a^6*b^18 + 5136*a^11*b^12 - 27912*a^9*b^14 - 17008*a^7*b^16 - 24*a^5*b^18 - 24902*a^10*b^12 + 47720*a^8*b^14 + 30494*a^6*b^16 - 340*a^4*b^18 - 2520*a^11*b^10 + 81968*a^9*b^12 - 73848*a^7*b^14 - 17008*a^5*b^16 - 64*a^3*b^18 + 23709*a^10*b^10 - 175528*a^8*b^12 + 97084*a^6*b^14 - 84*a^4*b^16 + 233*a^2*b^18 - 84992*a^9*b^10 + 270592*a^7*b^12 - 73848*a^5*b^14 + 5312*a^3*b^16 - 104*a*b^18 - 2367*a^10*b^8 + 200812*a^8*b^10 - 292164*a^6*b^12 + 47720*a^4*b^14 - 3179*a^2*b^16 + 16*b^18 + 25360*a^9*b^8 - 314344*a^7*b^10 + 270592*a^5*b^12 - 27912*a^3*b^14 + 1008*a*b^16 - 74484*a^8*b^8 + 386926*a^6*b^10 - 175528*a^4*b^12 + 13194*a^2*b^14 - 144*b^16 + 328*a^9*b^6 + 151472*a^7*b^8 - 314344*a^5*b^10 + 81968*a^3*b^12 - 3520*a*b^14 + 10880*a^8*b^6 - 167706*a^6*b^8 + 200812*a^4*b^10 - 24902*a^2*b^12 + 432*b^14 - 18696*a^7*b^6 + 151472*a^5*b^8 - 84992*a^3*b^10 + 5136*a*b^12 + 2288*a^8*b^4 + 36480*a^6*b^6 - 74484*a^4*b^8 + 23709*a^2*b^10 - 432*b^12 + 448*a^7*b^4 - 18696*a^5*b^6 + 25360*a^3*b^8 - 2520*a*b^10 + 3744*a^6*b^4 + 10880*a^4*b^6 - 2367*a^2*b^8 + 1408*a^7*b^2 + 448*a^5*b^4 + 328*a^3*b^6 - 512*a^6*b^2 + 2288*a^4*b^4 + 1408*a^5*b^2 + 256*a^6

aPoly has 331 terms
aFactor:
151632*a^12 - 556632*a^11 - 4029183*a^10 + 5710568*a^9 + 2456300*a^8 - 8614032*a^7 - 40073338*a^6 - 8614032*a^5 + 2456300*a^4 + 5710568*a^3 - 4029183*a^2 - 556632*a + 151632

alpha minpoly:
151632*x^12 - 556632*x^11 - 4029183*x^10 + 5710568*x^9 + 2456300*x^8 - 8614032*x^7 - 40073338*x^6 - 8614032*x^5 + 2456300*x^4 + 5710568*x^3 - 4029183*x^2 - 556632*x + 151632

beta minpoly:
16*x^12 - 100*x^10 + 615*x^8 - 2468*x^6 + 3186*x^4 - 1728*x^2 + 351

b1: 0.384369194474690789828391313191545078531026485622293845876818
c1: 0.923179463776614417385720356966316449484089189973058324576705
phiRad: 0.662140513907384715377580828031180874719561299125611851446970
phiDeg: 37.9378568915165327424085683976386530101541980484316304472328
x1: 0.823320543612498211991170347260337013360462788036571557634014
y1: 0.685480094624740433428838597906083940535454930045490771832954
a0: 0.147739677661122668796455680683973472704788159227719722064253
b0: 0.779540453991525714639344299875163938705026558053650718026349

numeric

sage: # This is an attempt to obtain numeric results for
sage: # http://math.stackexchange.com/q/688861/35416
sage: a0 = RDF(0.14778)
sage: b0 = RDF(0.77656)
sage: a1 = AA(14778/100000)
sage: b1 = AA(77656/100000)
sage: def cpm(v):
...       """Cross product matrix.
...       Multiplying with the returned matrix gives the
...       cross product with the input vector."""
...       x, y, z = v
...       return matrix(v.base_ring(), [
...           [0, z, -y],
...           [-z, 0, x],
...           [y, -x, 0]])
sage: def rank2to1(polyRing, rank2, c, sqrtIndex):
...       indices = [0, 1, 2]
...       del indices[sqrtIndex]
...       m1 = rank2.matrix_from_rows_and_columns(indices, indices).change_ring(polyRing)
...       m2 = c.matrix_from_rows_and_columns(indices, indices).change_ring(polyRing)
...       p = list((m1 + polyRing.gens()[0]*m2).det())
...       rank1 = sqrt(p[2])*rank2 + sqrt(-p[0])*c
...       return rank1
sage: def conicLineIsect(polyRing, conic, line, row = None, col = None, sqrtIndex = 1):
...       cpml = cpm(line)
...       rank2 = cpml.transpose()*conic*cpml
...       rank1 = rank2to1(polyRing, rank2, cpml, sqrtIndex=sqrtIndex)
...       if row is None and col is None:
...           return rank1
...       if row is not None:
...           row = rank1.row(row)
...           if col is None:
...               return row
...       if col is not None:
...           col = rank1.column(col)
...           if row is None:
...               return col
...       return (row, col)
sage: def simplify_full(x):
...       try:
...           f = f.apply_map
...       except:
...           return x.simplify_full()
...       else:
...           return f(simplify_full)
sage: RDFx = RDF['x']
sage: AAx = AA['x']
sage: def dehom(v):
...       return v[:-1]/v[-1]
sage: def myadj(m):
...       idxs = [[1, 2], [0, 2], [0, 1]]
...       return matrix(3, 3, lambda i, j: (-1)^(i+j)*m.matrix_from_rows_and_columns(idxs[j], idxs[i]).det())
sage: def compute1(polyRing, a, b, allvars=False):
...       field = polyRing.base_ring()
...       farpoint1 = vector(field, [1, b, 0])
...       farpoint2 = vector(field, [-b, 1, 0])
...       ellipse1 = diagonal_matrix(field, [a, 1, -1])
...       polar1 = ellipse1*farpoint1
...       polar2 = ellipse1*farpoint2
...       p1 = conicLineIsect(polyRing, ellipse1, polar1, col=0)
...       p2 = conicLineIsect(polyRing, ellipse1, polar2, row=0)
...       p3 = farpoint1.cross_product(p1).cross_product(farpoint2.cross_product(p2))
...       tr1 = matrix(field, [[-p3[2], 0, p3[0]], [0, -p3[2], p3[1]], [0, 0, p3[2]]])
...       tr2 = matrix(field, [[1, -b, 0], [b, 1, 0], [0, 0, sqrt(1+b^2)]])
...       tr3 = tr1*tr2
...       ellipse2 = (tr3.transpose()*ellipse1*tr3)
...       poly1 = (ellipse2.change_ring(polyRing) + polyRing.gens()[0]*ellipse1.change_ring(polyRing)).det()
...       roots1 = poly1.roots(field, False)
...       lam1 = min(roots1)
...       rank2 = ellipse2 + lam1*ellipse1
...       rank2adj = rank2.adjoint()
...       rank1 = rank2to1(polyRing, rank2, cpm(rank2adj.row(0)), 1)
...       g = rank1.column(2)
...       ellipse1adj = ellipse1.adjoint()
...       v = (g.row()*ellipse1adj*g.column())[0,0]
...       if allvars:
...           return locals().copy()
...       else:
...           return v
...
sage: res1 = compute1(AAx, a1, b1, allvars=True)
sage: view(res1['v'])
sage: def ell(m):
...       def f(x, y):
...           v = vector(m.base_ring(), [x, y, 1])
...           return (v.row()*m*v.column())[0,0] < 0
...       return f
sage: def draw1(polyRing, a, b):
...       res = compute1(polyRing, a, b, allvars=True)
...       r2 = sqrt(res['a'])
...       ctr = r2*dehom(res['tr2'].adjoint()*res['p3'])
...       p1 = r2*dehom(res['p1'])
...       p2 = r2*dehom(res['p2'])
...       p3 = r2*dehom(res['p3'])
...       # g = region_plot(ell(res['ellipse2'].change_ring(RDF)),
...       #                 (-4, 4), (-4, 4), plot_points=100, zorder=0, incol='cyan')
...       g = ellipse((0, 0), 1, r2, color='green', zorder=2)
...       g += ellipse(ctr, 1, r2, -atan(res1['b']), color='green', zorder=2)
...       g += point(p1, color='red', pointsize=50, zorder=3)
...       g += point(p2, color='red', pointsize=50, zorder=3)
...       g += point(p3, color='red', pointsize=50, zorder=3)
...       g += line2d([p1, p3], color='blue', zorder=2)
...       g += line2d([p2, p3], color='blue', zorder=2)
...       if 'g' in res:
...           g += line2d([r2*dehom(res['g'].cross_product(vector((1,0,0)))),
...                        r2*dehom(res['g'].cross_product(vector((0,1,0))))], color='magenta', zorder=2)
...       return g
...
sage: draw1(RDFx, a0, b0)

sage: compute1(RDFx, a0, b0)
-2.38178017815e-13
sage: [a0*1024, b0*1024]
[151.32672, 795.19744]
sage: plot(lambda x: compute1(RDFx, a0, x/1024), (790, 800))

sage: plot(lambda x: compute1(RDFx, x/1024, b0), (150, 155))

sage: compute1(RDFx, RDF(150/1024), b0)['roots1']
[-0.0112868672003, 0.0350920110764, 0.0455234757971]
sage: draw1(RDFx, RDF(150/1024), b0)

sage: poly2 = compute1(RDFx, RDF(151/1024), b0, allvars=True)['poly1']
sage: plot(lambda x: poly2(x), (-0.02, 0.05))

sage: def run1(cnt, R=None, Rx=None, prnt=True):
...       if R is None:
...           R = RDF
...           Rx = RDFx
...       if Rx is None:
...           Rx = R['x']
...       a1 = R(151/1024)
...       a3 = R(152/1024)
...       a2 = (a1 + a3)/2
...       b1 = R(792/1024)
...       b5 = R(794/1024)
...       b3 = (b1 + b5)/2
...       while cnt > 0:
...           cnt -= 1
...           v1 = compute1(Rx, a2, b1)
...           v3 = compute1(Rx, a2, b3)
...           v5 = compute1(Rx, a2, b5)
...           if (v5 - v3)*(v3 - v1) < 0:
...               # we certainly have an extremum in the range b1 .. b5
...               b2 = (b1 + b3)/2
...               b4 = (b3 + b5)/2
...               v2 = compute1(Rx, a2, b2)
...               v4 = compute1(Rx, a2, b4)
...               if (v4 - v3)*(v3 - v2) < 0:
...                   # 3 is more extreme than 2 and 4, so the
...                   b1 = b2
...                   b5 = b4
...               elif (v2 - v1)*(v3 - v2) < 0:
...                   b5 = b3
...                   b3 = b2
...               elif (v4 - v3)*(v5 - v4) < 0:
...                   b1 = b3
...                   b3 = b4
...               else:
...                   return ('Bad b adjustment: ', (v1, v2, v3, v4, v5), (b1, b2, b3, b4, b5))
...           else:
...               if v5 - v1 > 0:
...                   b3 = b5
...                   b5 = 2*b3 - b1
...               else:
...                   b3 = b1
...                   b1 = 2*b3 - b5
...           assert b3.parent() is R
...           v2 = compute1(Rx, a2, b3)
...           if v2 < 0:
...               v1 = compute1(Rx, a1, b3)
...               a3 = a2
...               if v1 > 0:
...                   a2 = (a1 + a3)/2
...               else:
...                   a2 = a1
...                   a1 = 2*a2 - a3
...           else:
...               v3 = compute1(Rx, a3, b3)
...               a1 = a2
...               if v3 < 0:
...                   a2 = (a1 + a3)/2
...               else:
...                   a2 = a3
...                   a3 = 2*a2 - a1
...           if prnt:
...               print((a2, a3-a1, b3, b5-b1))
...       return ((a2, a1, a3), (b3, b1, b5))
...
sage: run1(5)
(0.147705078125, 0.00048828125, 0.7734375, 0.00390625)
(0.147827148438, 0.000244140625, 0.775390625, 0.0078125)
(0.147766113281, 0.0001220703125, 0.771484375, 0.015625)
(0.147735595703, 6.103515625e-05, 0.779296875, 0.03125)
(0.147750854492, 3.0517578125e-05, 0.779296875, 0.015625)
((0.147750854492, 0.147735595703, 0.147766113281), (0.779296875, 0.771484375, 0.787109375))
sage: def commonDigits(a, b):
...       res = '?'
...       d = 2
...       Na = N(a, digits=d)
...       while Na == N(b, digits=d):
...           res = Na
...           d += 1
...           Na = N(a, digits=d)
...       return res
sage: res2 = run1(150, RealField(8*1024), prnt=False)
sage: print(commonDigits(res2[0][1], res2[0][2]))
0.1477396776611226687964556806839734727047881592
sage: print(commonDigits(res2[1][1], res2[1][2]))
0.779540453991525714639344299875163938705027
sage: res2 = run1(150, RealField(32*1024), prnt=False)
sage: print(commonDigits(res2[0][1], res2[0][2]))
0.1477396776611226687964556806839734727047881592
sage: print(commonDigits(res2[1][1], res2[1][2]))
0.779540453991525714639344299875163938705027
sage: a2 = res2[0][0]
sage: b2 = res2[1][0]
sage: draw1(RDFx, a2, b2).show(figsize=6)

sage: commonDigits(sqrt(res2[0][1]), sqrt(res2[0][2]))
0.3843691944746907898283913131915450785310264856
sage: commonDigits(sqrt(1-res2[0][1]), sqrt(1-res2[0][2]))
0.92317946377661441738572035696631644948408918997
sage: commonDigits(atan(res2[1][1]), atan(res2[1][2]))
0.662140513907384715377580828031180874719561
sage: commonDigits(atan(res2[1][1])*180/pi, atan(res2[1][2])*180/pi)
37.937856891516532742408568397638653010154
sage: res3 = compute1(a2.parent()['x'], a2, b2, allvars=True)
sage: N(sqrt(a2)*dehom(res3['tr2'].adjoint()*res3['p3']), digits=42)
(0.823320543612498211991170347260337013360463, 0.685480094624740433428838597906083940535455)