Here is how i would proceed. There will be some human calculation routine too, but the "technique" is a combined human+machine activity. The structure will come in a second, but so far, let us experiment something.
First of all, i asked sage about the Galois group of the number field generated by a root of the polynomial
$$
P = x^8 + x^7 - 7x^6 - 6x^5 + 15x^4 + 10x^3 - 10x^2 - 4x + 1
$$
to be factorized. It is a cyclic group of order eight. I also asked for the subfields, and sage mentions a subfield of degree four over $\Bbb Q$, related to the polynomial
$$
Q = x^4 + x^3 - 6 x^2 - x + 1
\ .
$$
(Code and results are postponed to not disturb here the flow. A structural way to see it follows below.)
Let us compute its roots. This is a "twisted reciprocal" polynomial, one can substitute $X=ix$ first to get a reciprocal polynomial. The substitution
$$
t = x-\frac 1x
$$
allows to relate $Q/x^2$ with
$$
T = t^2 + t -4\ .
$$
(This polynomial is also present in the sage list of subfields...)
The roots of $T$ are $t_{\pm}=(-1\pm \sqrt{17})/2$. The polynomials $x^2-t_\pm x-1$ have roots
$$
\begin{aligned}
x_{\boxed\pm}
&=\frac 12\left(t_\pm\boxed\pm\sqrt{t_\pm^2+4}\right)
\\
&=\frac 12\left(t_\pm\boxed\pm\sqrt{8-t_\pm}\right)
\\
&=
\frac 14\left(-1\pm\sqrt{17}\right)
\boxed\pm
\frac 12\sqrt{\frac 12\left(17\mp\sqrt{17}\right)}
\ .
\\
4x_{\boxed\pm}
&=
\left(-1\pm\sqrt{17}\right)
\boxed\pm
\sqrt{2\left(17\mp\sqrt{17}\right)}
\ .
\end{aligned}
$$
Time for the structure. I learned it from an older book of Hans Rademacher,
typed it for myself some years ago. The idea is to use Gaussian periods.
Let $K=\Bbb Q(w)$ be the cyclotomic field generated by the primitive root of order $17$, here fixed and denoted by $w$ instead of the longer $\zeta_{17}$, which is hard to type. The Galois group of $K$ over $\Bbb Q$ has elements sending $w$ to some other primitive root $w^k$, so $k\in 1,2,\dots 16$.
My choice is $3$ as a generator of the cyclic multiplicatie group $(\Bbb Z/17)^\times$, so let us use the Frobenius morphism
$$ Fw = w^3\ .$$
So $F^2w = FFw=Fw^3=(w^3)^3=w^9$, and $F^3w = FF^2w=Fw^9=(w^3)^9 = w^{27}$, and in general $F^jw = w^{3^j}$.
The idea of the Gaussian periods is simple. We associate the following sums:
$$
\begin{aligned}
S &= F^0w + Fw + F^2w + F^3w+\dots+F^{14}w+F^{15}w\\
&=w+w^2+w^3+\dots+w^{16}\text{ (in other order)}\\
&=-1\ ,\\[2mm]
&\qquad\text{ and we split it in two subsums,}\\[2mm]
S_0 &= F^0w + F^2w + F^4w+\dots+F^{14}w\\
S_1 &= Fw + F^3w + F^5w+\dots+F^{15}w\\[2mm]
&\qquad\text{ and we split each in two subsums,}\\[2mm]
S_{00} &= F^0w + F^4w + F^8w+F^{12}w\\
S_{10} &= F^2w + F^6w + F^{10}w+F^{14}w\\[2mm]
S_{01} &= Fw + F^5w + F^9w+\dots+F^{13}w\\
S_{11} &= F^3w + F^7w + F^{11}w+\dots+F^{15}w\\[2mm]
&\qquad\text{ and we split each in two subsums,}\\[2mm]
S_{000}&= F^0w + F^8w\\
S_{100}&= F^4w + F^{12}w\\
&\qquad\text{and so on...}
\end{aligned}
$$
The indices are my choice, they suggest binary written numbers. So $S_{01}$ is a sum of terms $F^jw$ with $j$ congruent to (the binary phone number) $01_{[2]}=1$ modulo $4$, $S_{101}=F^5w+F^{13}w$ and so on.
There are obvious splitting relations, for instance $S_0+S_1=S=-1$, $S_{00}+S_{10}=S_0$, sum relations. But there are also product relations.
Instead of writing them here explicitly, a hard jax job, i will use sage to "type" them. Sage is a mathematically structurally thinking language, so i hope the dry exposition gives the information in a human way.
sage: K.<w> = CyclotomicField(17)
sage: K
Cyclotomic Field of order 17 and degree 16
sage: w^17
1
sage: S = sum( [ w^(3^j) for j in [0..15] ] )
sage: S0 = sum( [ w^(3^j) for j in [0..15] if j%2 == 0 ] )
sage: S1 = sum( [ w^(3^j) for j in [0..15] if j%2 == 1 ] )
sage: S00 = sum( [ w^(3^j) for j in [0..15] if j%4 == 0 ] )
sage: S01 = sum( [ w^(3^j) for j in [0..15] if j%4 == 1 ] )
sage: S10 = sum( [ w^(3^j) for j in [0..15] if j%4 == 2 ] )
sage: S11 = sum( [ w^(3^j) for j in [0..15] if j%4 == 3 ] )
sage: S000 = sum( [ w^(3^j) for j in [0..15] if j%8 == 0 ] )
sage: S001 = sum( [ w^(3^j) for j in [0..15] if j%8 == 1 ] )
sage: S010 = sum( [ w^(3^j) for j in [0..15] if j%8 == 2 ] )
sage: S011 = sum( [ w^(3^j) for j in [0..15] if j%8 == 3 ] )
sage: S100 = sum( [ w^(3^j) for j in [0..15] if j%8 == 4 ] )
sage: S101 = sum( [ w^(3^j) for j in [0..15] if j%8 == 5 ] )
sage: S110 = sum( [ w^(3^j) for j in [0..15] if j%8 == 6 ] )
sage: S111 = sum( [ w^(3^j) for j in [0..15] if j%8 == 7 ] )
sage: S0+S1, S0*S1
(-1, -4)
sage: S0.minpoly()
x^2 + x - 4
sage: S00*S10, S01*S11
(-1, -1)
sage: S00.minpoly()
x^4 + x^3 - 6*x^2 - x + 1
sage: S10.minpoly()
x^4 + x^3 - 6*x^2 - x + 1
sage: S01.minpoly()
x^4 + x^3 - 6*x^2 - x + 1
sage: S11.minpoly()
x^4 + x^3 - 6*x^2 - x + 1
sage: S000*S100
w^14 + w^12 + w^5 + w^3
sage: S00, S10, S01, S11
(-w^15 - w^14 - w^12 - w^11 - w^10 - w^9 - w^8 - w^7 - w^6 - w^5 - w^3 - w^2 - 1,
w^15 + w^9 + w^8 + w^2,
w^14 + w^12 + w^5 + w^3,
w^11 + w^10 + w^7 + w^6)
sage: S000*S100 == S01
True
(The first long expression staring with -w^15 - w^14 ...
is only the computer version to hide the "missing terms" of the minimal polynomial of $w$, the four that not appear in the list, since it thinks it is a good idea to not use $w^{16}$...)
Note that in the "naive approach" from the beginning was also dealing with polynomials that appear as minimal polynomials of the one or the other $S_?$.
The OP wants to understand the splitting of the polynomial $P$, which is the minimal polynomial (over rationals) of
$$S_{000}=w+w^{16}=w+\frac 1w\ .$$
And now we can give four factors of $P$. With the computer aid of sage:
sage: R.<x> = PolynomialRing(K)
sage: R
Univariate Polynomial Ring in x over Cyclotomic Field of order 17 and degree 16
sage: (x-S000)*(x-S100)
x^2 + (w^15 + w^14 + w^12 + w^11 + w^10 + w^9 + w^8 + w^7 + w^6 + w^5 + w^3 + w^2 + 1)*x + w^14 + w^12 + w^5 + w^3
sage: (x-S000)*(x-S100) == x^2 - S00*x + S01
True
sage: (x-S010)*(x-S110)
x^2 + (-w^15 - w^9 - w^8 - w^2)*x + w^11 + w^10 + w^7 + w^6
sage: (x-S010)*(x-S110) == x^2 - S10*x + S11
True
sage: (x-S001)*(x-S101)
x^2 + (-w^14 - w^12 - w^5 - w^3)*x + w^15 + w^9 + w^8 + w^2
sage: (x-S001)*(x-S101) == x^2 - S01*x + S10
True
sage: (x-S011)*(x-S111)
x^2 + (-w^11 - w^10 - w^7 - w^6)*x - w^15 - w^14 - w^12 - w^11 - w^10 - w^9 - w^8 - w^7 - w^6 - w^5 - w^3 - w^2 - 1
sage: (x-S011)*(x-S111) == x^2 - S11*x + S00
True
sage: (x-S00)*(x-S10)*(x-S01)*(x-S11)
x^4 + x^3 - 6*x^2 - x + 1
And here explicitly, because it is structural (the "traces", coefficients in $x$ on the RHSs are easy, we splitted them as the idea to proceed, for the computation of the free coefficients there is also a Galois pattern...):
$$
\begin{aligned}
(x-S_{000})(x-S_{100}) &= x^2 - S_{00}x+ S_{01}\ ,\\
(x-S_{010})(x-S_{110}) &= x^2 - S_{10}x+ S_{11}\ ,\\
(x-S_{001})(x-S_{101}) &= x^2 - S_{01}x+ S_{10}\ ,\\
(x-S_{011})(x-S_{111}) &= x^2 - S_{11}x+ S_{00}\ ,\\[2mm]
&\text{where}\\
(x-S_{00})
(x-S_{10})
(x-S_{01})
(x-S_{11})
&=x^4 + x^3 - 6x^2 - x + 1\ ,
\end{aligned}
$$
and its roots were computed in advance.
I have to submit, but there remains a lot to say, maybe just mention that this structure made the 19 years old Gauss to choose mathematics, not physics, one main reason for still missing a unification of matter in this field...
Postponed:
The sage code used to obtain the polynomial $Q$, the Galois theoretical information, and the decomposition of $P$.
sage: R.<x> = PolynomialRing( QQ )
sage: P = x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
sage: K.<v> = NumberField( P )
sage: v.minpoly()
x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
sage: G = K.galois_group()
sage: G
Galois group of Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
sage: G.is_cyclic()
True
sage: K.subfields()
[
(Number Field in v0 with defining polynomial x + 1, Ring morphism:
From: Number Field in v0 with defining polynomial x + 1
To: Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
Defn: -1 |--> -1, None),
(Number Field in v1 with defining polynomial x^2 + x - 4, Ring morphism:
From: Number Field in v1 with defining polynomial x^2 + x - 4
To: Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
Defn: v1 |--> -v^7 - v^6 + 6*v^5 + 6*v^4 - 10*v^3 - 9*v^2 + 5*v + 1, None),
(Number Field in v2 with defining polynomial x^4 + x^3 - 6*x^2 - x + 1, Ring morphism:
From: Number Field in v2 with defining polynomial x^4 + x^3 - 6*x^2 - x + 1
To: Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
Defn: v2 |--> v^4 - 4*v^2 + v + 2, None),
(Number Field in v3 with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1, Ring morphism:
From: Number Field in v3 with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
To: Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
Defn: v3 |--> v, Ring morphism:
From: Number Field in v with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
To: Number Field in v3 with defining polynomial x^8 + x^7 - 7*x^6 - 6*x^5 + 15*x^4 + 10*x^3 - 10*x^2 - 4*x + 1
Defn: v |--> v3)
]
sage: K2, morphism, _ = K.subfields()[2]
sage: K2
Number Field in v2 with defining polynomial x^4 + x^3 - 6*x^2 - x + 1
sage: K2.gen()
v2
sage: for f, mul in P.base_extend(K2).factor():
....: print f
....:
x^2 - v2*x - 1/2*v2^3 + 3*v2 - 3/2
x^2 + (-v2^3 - v2^2 + 6*v2 + 1)*x - 1/2*v2^3 - v2^2 + 2*v2 + 3/2
x^2 + (1/2*v2^3 - 3*v2 + 3/2)*x + v2^3 + v2^2 - 6*v2 - 1
x^2 + (1/2*v2^3 + v2^2 - 2*v2 - 3/2)*x + v2
Best Answer
The Fundamental Theorem of Algebra states precisely that:
This result does not hold in general. In some cases, you cannot factor them (you cannot in $\mathbb{R}$, for instance, where $x^2+1$ cannot be written as a product of linear terms; over $\mathbb{Q}$, there are polynomials of arbitrarily high degree that cannot be factored at all, let alone into a product of linear terms). In some cases, the factorization is not unique: in the quaternions, the polnomial $x^2+1$ can be factored in infinitely many essentially distinct ways as a product of linear terms, e.g., $x^2+1 = (x+i)(x-i) = (x+j)(x-j) = (x+k)(x-j) = \cdots$.
A field $F$ is said to be algebraically closed if and only if every nonconstant polynomial $p(x)$ with coefficients in $F$ has at least one root. It is then easy to verify that $F$ is algebraically closed if and only if every nonzero polynomial can be factored into a product of linear terms, as above. Subject to some technical assumptions (The Axiom of Choice), every field is contained in an algebraically closed field, and for every field there is a "smallest" algebraically closed field that contains it, so there is a "smallest" larger field $K$ where you can guarantee that every polynomial factors. This does not hold for arbitrary sets of coefficients (e.g., the quaternions, because they are non-commutative; or rings with zero divisors; and others).
Once you go beyond one variable, it is no longer true that a polynomial can always be factored into terms of degree one, even in $\mathbb{C}$. For example, the polynomial $xy-1$ in $\mathbb{C}[x,y]$ cannot be written as a product of polynomials of degree $1$, $(ax+by+c)(rx+sy+d)$. If you could, then $ar=bs=0$, and we cannot have $a=b=0$ or $r=s=0$, so without loss of generality we would have $b=0$ and $r=0$, so we would have $(ax+c)(sy+d) = xy-1$. Then $ad=0$, so $d=0$ (since $a\neq 0)$, but then we cannot have $cd=-1$. So no such factorization is possible.