The resolvent method for determining Galois groups works over any field – not just $\mathbb Q$ – and is enough to solve this question. To illustrate how consider coordinate values $a$ and $b$ in the following unit-distance embedding of $F_{56}A$ with $D_{14}$ symmetry:

$a$ and $b$ satisfy
$$a^{4} \left(16 t^{4} + 16 t^{2}\right) - 64 a^{3} t^{3} + a^{2} \left(- 8 t^{4} + 24 t^{2}\right) + 16 a t + t^{4} + 5 t^{2} - 4=0\\
b^{4} \left(16 t^{4} + 16 t^{2}\right) + 64 b^{3} t^{3} + b^{2} \left(- 8 t^{4} + 24 t^{2}\right) - 16 b t + t^{4} + 5 t^{2} - 4=0$$
where $t=\tan\frac\pi{14}$. To compute the resolvent polynomial in $\mathbb Q[t][y]$ exactly I use Stauduhar's method of computing its coefficients numerically to high precision and "rounding" them to the nearest element in $\mathbb Q[t]$ – where rounding means applying an integer relation algorithm. Concretely, using mpmath and SymPy for the resolvent of the $a$-polynomial with respect to $x_1-x_2$:
#!/usr/bin/env python3
from mpmath import *
from itertools import combinations, permutations
mp.dps = 100
from sympy import symbols, primitive
tt, y = symbols('t y')
t = tan(pi/14)
roots1 = polyroots([16t4 + 16t2, -64*t3, -8t4 + 24t2, 16*t, t4 + 5t*2 - 4]) # the a-polynomial
[16t4 + 16t2, 64*t3, -8t4 + 24t2, -16*t, t4 + 5t*2 - 4] is the b-polynomial
resultant_roots = [a-b for (a,b) in permutations(roots1, 2)]
ydeg = len(resultant_roots)
P = yydeg
for r in range(1, ydeg+1):
Z = chop(fsum(fprod(x) for x in combinations(resultant_roots, r)), tol=1e-70)
if Z == 0:
continue
basis = [1, t, t2, t3, t4, t5, Z]
coords = pslq(basis, maxcoeff=1012, maxsteps=1000)
print(r, coords, nstr(fdot(basis, coords)))
Z_expr = sum(ctti / -coords[-1] for (i, c) in enumerate(coords[:-1]))
P += (-1)r Z_expr * y**(ydeg-r)
P = primitive(P)[1]
print("ring r=(0,t),y,dp; minpoly=7t6-35t4+21t2-1;")
print(f"poly p={P};")
print("factorize(p);")
$P$ at the end is the exact resolvent polynomial, and the last lines of the program's output are Singular commands to factor the resolvent in $\mathbb Q[t]$.
> ring r=(0,t),y,dp; minpoly=7t6-35t4+21t2-1;
> poly p=91*t**4*y**10 - 70*t**4*y**8 + 3570*t**4*y**6 - 581*t**4*y**4 - 113862*t**4*y**2 - 230846*t**4 - 462*t**2*y**10 + 392*t**2*y**8 - 17752*t**2*y**6 + 2982*t**2*y**4 + 563220*t**2*y**2 + 1142232*t**2 + 16*y**12 + 167*y**10 - 2*y**8 + 9750*y**6 - 1909*y**4 - 311542*y**2 - 633178;
> factorize(p);
[1]:
_[1]=16
_[2]=y4+(35/16t4-91/8t2+119/16)*y2+(-399/16t4+987/8t2-1099/16)
_[3]=y8+(7/2t4-35/2t2+3)*y6+(175/16t4-441/8t2+831/16)*y4+(1225/8t4-1505/2t2+3083/8)*y2+(3521/16t4-8757/8t2+10117/16)
[2]:
1,1,1
The orbit-length partition corresponding to $x_1-x_2$ is $(4,8)$. Now I can consult Soicher and McKay (1985), p. 277 for the possible Galois groups of the $a$-polynomial; the only possibility for degree $4$ and partition $(4,8)$ is $D_4$, the dihedral group of order $8$. Checking LMFDB then shows that the $a$-polynomial splits over a quadratic extension of $\mathbb Q[t]$, which in this case is $\mathbb Q[t,u=\sqrt{-21t^4+98t^2+71}]$. Indeed the $b$-polynomial also splits over $\mathbb Q[t,u]$ and
$$z_1=2\sqrt{14(31t^5u-20t^4-154t^3u+104t^2+87tu-68)}\qquad z_2=7t(t^2-3)^2-4u$$
$$a=\frac{z_1+z_2}{64}\qquad b=\frac{z_1-z_2}{64}$$
Returning to the original question, I compute the resolvent of $(1)$ with respect to $x_1+x_2$:
#!/usr/bin/env python3
from mpmath import *
from itertools import combinations
mp.dps = 100
from sympy import symbols, primitive
tt, y = symbols('t y')
t = tan(pi/7)
roots1 = polyroots([256, 16t5 - 320t3 + 304t, -16t2 - 144, -32t5 + 656t3 - 912t, -4t4 + 112t2 - 188, 7t5 - 136*t3 + 57t, -t4 + 13t**2 + 14])
print(roots1[0])
resultant_roots = [a+b for (a,b) in combinations(roots1, 2)]
ydeg = len(resultant_roots)
P = yydeg
for r in range(1, ydeg+1):
Z = chop(fsum(fprod(x) for x in combinations(resultant_roots, r)), tol=1e-70)
if Z == 0:
continue
basis = [1, t, t2, t3, t4, t5, Z]
coords = pslq(basis, maxcoeff=1012, maxsteps=1000)
print(r, coords, nstr(fdot(basis, coords)))
Z_expr = sum(ctti / -coords[-1] for (i, c) in enumerate(coords[:-1]))
P += (-1)r Z_expr * y**(ydeg-r)
P = primitive(P)[1]
print("ring r=(0,t),y,dp; minpoly=t6-21t4+35t2-7;")
print(f"poly p={P};")
print("factorize(p);")
then factor:
> ring r=(0,t),y,dp; minpoly=t6-21t4+35t2-7;
> poly p=160*t**5*y**14 - 328*t**5*y**12 + 38*t**5*y**10 - 425*t**5*y**8 + 123*t**5*y**6 + 16*t**5*y**4 - 65*t**5*y**2 - 13*t**5 - 80*t**4*y**13 - 14*t**4*y**11 + 516*t**4*y**9 - 118*t**4*y**7 + 91*t**4*y**5 + 27*t**4*y**3 - 18*t**4*y - 3200*t**3*y**14 + 6576*t**3*y**12 - 864*t**3*y**10 + 8922*t**3*y**8 - 2902*t**3*y**6 - 520*t**3*y**4 + 1394*t**3*y**2 + 258*t**3 + 1632*t**2*y**13 - 612*t**2*y**11 - 8040*t**2*y**9 + 1824*t**2*y**7 - 3218*t**2*y**5 - 314*t**2*y**3 + 328*t**2*y + 3040*t*y**14 - 7944*t*y**12 + 3802*t*y**10 - 6877*t*y**8 + 1759*t*y**6 + 3832*t*y**4 - 2653*t*y**2 - 161*t + 512*y**15 - 592*y**13 - 2294*y**11 + 3668*y**9 - 3418*y**7 + 3443*y**5 - 469*y**3 - 406*y;
> factorize(p);
[1]:
_[1]=512
_[2]=y15+(5/16t5-25/4t3+95/16t)*y14+(-5/32t4+51/16t2-37/32)*y13+(-41/64t5+411/32t3-993/64t)*y12+(-7/256t4-153/128t2-1147/256)*y11+(19/256t5-27/16t3+1901/256t)*y10+(129/128t4-1005/64t2+917/128)*y9+(-425/512t5+4461/256t3-6877/512t)*y8+(-59/256t4+57/16t2-1709/256)*y7+(123/512t5-1451/256t3+1759/512t)*y6+(91/512t4-1609/256t2+3443/512)*y5+(1/32t5-65/64t3+479/64t)*y4+(27/512t4-157/256t2-469/512)*y3+(-65/512t5+697/256t3-2653/512t)*y2+(-9/256t4+41/64t2-203/256)*y+(-13/512t5+129/256t3-161/512t)
[2]:
1,1
The orbit-length partition is $(15)$; Soicher and McKay indicate that the Galois group is one of $PSL(2,5),PGL(2,5),A_6,S_6$. However, LMFDB shows that none of these four groups allow subfields, so the conclusion is that $(1)$ has no factorisation over an extension of $\mathbb Q[t]$. Indeed, the results of further resolvent tests with respect to $p_\Delta=\prod_{i<j}(x_i-x_j)$ and $x_1+x_2+x_3+p_\Delta$ (cf. this "tutorial" by Curtis Bright) shows that the Galois group of $(1)$ over $\mathbb Q[t]$ is exactly $S_6$.
IsSolvable
on the degree-$36$ polynomial $p$ takes a while but returns False. If $(1)$ split, the extension would be of degree $2$ or $3$, and since $t$ is the root of a solvable polynomial, $p$ would be solvable, which contradicts the first result. I am not sure if this logic is sound, though. – Parcly Taxel Apr 11 '21 at 03:44