http://math.stackexchange.com/a/861455/35416

mx861389a.sage

def rotmat(a, t=(0,0,1), l=False):
    # Rotate by a degrees, then shift origin to point t.
    # If l=True then transform lines not points.
    x, y, z = t
    q = QQ(a/360)
    w = QQbar.zeta(q.denominator())^(q.numerator())
    c = AA(z*w.real())
    s = AA(z*w.imag())
    m = matrix(AA, [[c, -s, x], [s, c, y], [0, 0, z]])
    if l:
        m = m.adjoint().transpose()
    return m
B = vector((-1, 0, 1))
C = vector((1, 0, 1))
a = vector((0, 1, 0))
b = rotmat(-80, C, True)*a
c = rotmat(80, B, True)*a
A = b.cross_product(c)
BD = rotmat(30, B, True)*a
CE = rotmat(-20, C, True)*a
D = b.cross_product(BD)
E = c.cross_product(CE)
DE = D.cross_product(E)
cosX = DE[:2].normalized().dot_product(BD[:2].normalized())
print cosX.minpoly()*1216 # 1216*x^6 - 3264*x^4 + 2916*x^2 - 867

Q80.<sin80, cos80> = QQ[]
I80 = ideal(sin80^2 + cos80^2 - 1)
rmd = dict()
rmd[80] = (cos80, sin80)
rmd[90] = (0, 1)
def addrmd(a, b):
    global rmd
    ca, sa = rmd[abs(a)]
    cb, sb = rmd[abs(b)]
    if a < 0: sa = -sa
    if b < 0: sb = -sb
    rmd[a + b] = (ca*cb - sa*sb, sa*cb + ca*sb)
addrmd(90, -80)
addrmd(10, 10)
addrmd(10, 20)
def rotmat(a, t=(0,0,1), l=False):
    x, y, z = t
    c, s = rmd[abs(a)]
    if a < 0: s = -s
    m = matrix(Q80, [[c, -s, x], [s, c, y], [0, 0, z]])
    if l:
        m = m.adjoint().transpose()
    return m
B = vector((-1, 0, 1))
C = vector((1, 0, 1))
a = vector((0, 1, 0))
b = rotmat(-80, C, True)*a
c = rotmat(80, B, True)*a
A = b.cross_product(c)
BD = rotmat(30, B, True)*a
CE = rotmat(-20, C, True)*a
D = b.cross_product(BD)
E = c.cross_product(CE)
DE = D.cross_product(E)
cosX2 = (DE[:2].dot_product(BD[:2]))^2/(DE[:2]*DE[:2])/(BD[:2]*BD[:2])
sinX2 = 1 - cosX2
tanX2 = sinX2/cosX2
print (tanX2.numerator().reduce(I80)/16).factor()
print (tanX2.denominator().reduce(I80)/16).factor()
# (4) * cos80^2 * (2*cos80^2 - 1)^2 * (8*cos80^4 - 8*cos80^2 + 1)^2
# (-1) * (cos80 - 1) * (cos80 + 1) * (32*cos80^6 - 32*cos80^4 + 8*cos80^2 - 1)^2
tanX = (2*cos80*(2*cos80^2-1)*(8*cos80^4-8*cos80^2+1))/(
       sin80*(32*cos80^6-32*cos80^4+8*cos80^2-1))
assert (tanX^2 - tanX2).numerator().reduce(I80).is_zero()