An orthogonal companion matrix Let $P\in{\mathbb R}[X]$ be a monic polynomial with roots on the unit circle. For the problem below, we may assume wlog that the roots are simple and distinct from $\pm1$. It can be shown that there exists a matrix $M\in{\bf SO}_n({\mathbb R})$, whose characteristic polynomial is $P$ (an orthogonal companion matrix of $P$, in short OCM). See for instance Exercise 99 on my list http://www.umpa.ens-lyon.fr/~serre/DPF/exobis.pdf .
Regretfully, this exercise uses the square root of Hermitian positive definite matrices, which cannot be computed in finitely many operations.
Does there exist a construction of an OCM that uses only finitely many elementary operations (including the square root of complex numbers) ?
Thanks to the reduction to Hessenberg form, which can be done in finite time and which preserves the orthogonal group, we may restrict our attention to a Hessenberg orthogonal matrix $M$. It writes
$$\left( \begin{array}{ccccc}
 c_1 & s_1c_2 & s_1s_2c_3 & s_1s_2s_3c_4 &  \ldots   \\\\
-s_1 & c_1c_2 & c_1s_2c_3 & c_1s_2s_3c_4 &  \ldots   \\\\ 
0 & -s_2 & c_2c_3 & s_2s_3c_4 &  \ldots     \\\\  
0 & 0 & -s_3 & c_3c_4 & \ldots     \\\\ 
0 & 0 & 0 & -s_4 & \ldots    
 \end{array} \right)$$
where $(c_j,s_j)$ are cosine/sine pairs.
 A: I'd do this in three steps:


*

*Find any $2n \times 2n$ matrix $A$ whose eigenvalues are $e^{\pm i \theta}$.

*Find a positive definite quadratic form preserved by $A$. In equations, we want $A P A^T = P$. 

*Find an orthonormal basis for $P$, using the Gram-Schimdt algorithm. In equations, we want $S P S^T = \mathrm{Id}$. 
Then $S A S^{-1}$ is orthogonal and has the required eigenvalues.

I can think of two ways to do step 2. The first is more purely algebraic, the second I think would be much easier to implement.
Algebra: Let $f(x) = \prod_{j=1}^{n} (x-e^{i \theta_j}) (x - e^{-i \theta_j}) = \prod (x^2 - 2 \cos \theta_j + 1)$ be your characteristic polynomial. Let $V$ be the ring $\mathbb{R}[x]/f(x)$. Note $1$, $x$, ..., $x^{2n-1}$ is a basis for this ring, in which multiplication can be written down algebraically in terms of the coefficients of $f$. Also, multiplication by $x$ has the desired eigenvalues, so that accomplishes part 1.
For $y \in T$, let $T(y)$ be the trace of multiplication by $y$. Also, let $y \mapsto \overline{y}$ be the automorphism of $V$ induced by $x \mapsto x^{-1}$. Again, both of these can be written down, in the monomial basis, algebraically in terms of the coefficients of $f$.
Then $\langle y,z \rangle = T(y*\overline{z})$ is the desired positive definite quadratic form. Namely, observe that the ring $V$ is isomorphic to $\mathbb{C}^{\oplus n}$. In terms of this isomorphism, $\langle (z_1, \ldots, z_n), (z_1, \ldots, z_n) \rangle = 2 \sum |z_i|^2$.
In practice: The condition that $APA^T = P$ is a linear condition on $P$. Let $W$ be the subspace of the vector space of symmetric matrices where this condition is satisfied; finding $W$ is just algebra. Now, our goal is to find a positive definite element of $W$. For large $N$, $(1/N) \left( \mathrm{Id} + AA^T + A^2 (A^{T})^2 + \cdots + A^{N-1} (A^T)^{N-1} \right)$ is positive definite and is near $W$. I would guess that the orthogonal projection of this matrix onto $W$ would probably be positive definite for large $N$.
A: You can use the Cayley transform and reduce this problem to generating a skew-symmetric matrix with a prescribed characteristic polynomial. For example, this works in the $3\times 3$ case (although when $n$ is odd, $1$ is always an eigenvalue). I have not thought about it thoroughly, but presumably, methods used in Inverse Symmetric Eigenvalue Problem should apply in the skew-symmetric case. 
