Given any real monic polynomial $p(x)$ with only real zeros, there is a finite algebraic (no root finding) algorithm to construct a real symmetric and tri-diagonal matrix $A$ so that $p(x)$ is the characteristic polynomial of $A$. This proves $p(x)$ has only real roots and $A$ gives a certificate. Interestingly this also solves the finite version of Hilbert-Polya problem, given real rooted $p(x)$ find a Hermitian operator whose eigenvalues gives the roots of $p(x)$.

If $p,q$ are two real polynomials with $deg(p)=deg(q)+1$ and only real and simple roots which interlace , then applying Euclidean algorithm (divide out and flip, and this is the same as the usual Strum computation), we get a recursion.
If we write $p/q=x+a-b/(q/q_1),$ where $a,b$ are real and $q_1$ is a monic polynomial with $deg(q_1)=deg(q)-1$ which we can, we will have $b>0$ and the roots of $q_1$ are all real and simple and interlaces those of $q$. By induction, we get a continued fraction expansion applying this to $(p,p'/n)$,

$\frac{p(x)}{p'(x)/n}=x+a_1-\cfrac{b_1}{x+a_2-\cfrac{b_2}{...-\cfrac{...}{x+a_{n-1}-\cfrac{b_{n-1}}{x+a_n}}}},$

with $a_j$ real and $b_j>0.$ But a generic continued fraction is a Cauchy interlacing ratio of a symmetric tri-diagonal matrix which can be made symmetric, so we have

$\frac{p(x)}{p'(x)/n}=\frac{\det(xI-A)}{\det(xI-A_1)},$ where $A$ is the tridiagonal matrix

$A=\begin{pmatrix} -a_1 & \sqrt{b_1} & 0 & ... & 0\cr \sqrt{b_1} & -a_2 & \sqrt{b_2} & 0 &
0\cr 0 & \sqrt{b_2} & -a_3 & \sqrt{b_3} & 0 \cr ... & ... &.... & .... &...
\cr 0 &...& \sqrt{b_{n-2}} & -a_{n-1} & \sqrt{b_{n-1}} \cr 0 & ... & 0 & \sqrt{b_{n-1}} & -a_n
\end{pmatrix},$

and $A_1$ is $A$ with the first row and column deleted. This in fact gives $2^{n-1}$ solutions as we can replace each $\sqrt{b_j}$ by $-\sqrt{b_j}$ independently. We can also replace $p'(x)/n$ by any polynomial $q(x)$ with roots interlacing those of $p$. If we start with a $p$ with non real roots, either the algorithm fail or we have some $b_j \le 0$. If $p$ has multiple root, it will be cancelled out by $p'(x)$. To certify that a given $g(z)$ has all roots on the unit circle, we apply the above to $p(x)=(x+i)^{deg(g)}g(\frac{x-i}{x+i})$.

1more comment