Unitary equivalence to permutation matrix I would like to find a check whether a given matrix is unitarily equivalent to a permutation matrix. I split this question into two parts below. I expect the solution will use the second part, but if not that's fine.
${\bf(1)}$ Let $M$ be a unitary matrix. How can I check if there exists a unitary $U$ such that $U M U^\dagger = P$ is a permutation matrix?
I know that a permutation matrix has eigenvalues that form complete sets of roots of unity. The question then boils down to.
${\bf(2)}$ Given a set of complex numbers on the unit circle $S = \{z_1,z_2,z_3,\ldots\}$, is there a simple test to check if this contains complete sets of roots of unity. In other words, check that $S= S_1 \cup S_2 \cup \cdots$, where $S_i = \{e^{2\pi i \frac{n}{q}} | n \in \{0,\ldots,q-1\} \}$, for $q\in \mathbb{N}$.
Examples
$\{1,1,1,-1\} = \{1\} \cup \{1\} \cup \{1,-1\}$ and so this does have complete sets of roots of unity.
$\{1,i,-i,1\}$ does not contain complete sets of roots of unity.
 A: Main ideas
i.)  Unitary matrices (including permutation matrices as a special case) are normal, so they are unitarily similar to a diagonal matrix $D$, which is itself unitary.  Thus to test whether we can have $UMU^* = P$, with given $M$ and some $\text{n x n}$ Permutation matrix $P$ it suffices to check that there is some (unitary) diagonal matrix $\Lambda$ that they could both be similar to.  Since diagonalizability is guaranteed, it suffices to check some algebraic invariant like the characteristic polynomial.  
ii.)  With $M = QDQ^*$ and $P= V\Lambda V^*$ The algebraic invariant of interest is in fact the trace.  I.e. examine $\text{trace}\big(M^k\big)= \text{trace}\big(D^k\big)$ for $k\in\{1,2,...,n\}$ and checking whether it could be equal to  $\text{trace}\big(P^k\big) = \text{trace}\big(\Lambda^k\big)$ where $\Lambda$ is diagonal and has the constraint that it only has complete sets of $r_i$th roots of unity, with each $r_i \in \{1,2,...,n\}$. The trace is of interest since it counts closed walks on graphs and any permutation matrix $P$ is an adjacency matrix for a graph.  Of course $P$, if it exists, may be written in a particularly simple form-- counting the number of eigenvalues of 1 tells us how many communicating classes it has (up to graph isomorphism this should be thought of as disjoint blocked cyclic permutation submatrices);  this result follows e.g. from Perron-Frobenius Theory.  We can even do better and construct a simple $P$ here, if it exists, by paying attention to the specific rth roots of unity and associated 'multiplicities'.  The algorithm below actually gives a way to construct a 'nice' Permutation matrix $P$ that is similar to $M$ (and fails if some similarity is not possible).  
Part (2) of the OP is a good question in that it maps question (1) down to a simpler problem, but implies the result for (1).  I give an algorithmic solution to (2), below.  Once that is understood, one can re-run the below algorithm verbatim, except compute $\text{trace}\big(M^k \big)$ and $\text{trace}\big(M^k - \Lambda^k\big)$ instead of $\text{trace}\big(D^k\big)$ and $\text{trace}\big(D^k - \Lambda^k\big)$.   
For avoidance of doubt for (2), instead of a "set" I interpret this as a "multi-set" $S$ and $\big \vert S\big \vert = n$.  With part 1 in mind this means we are interested solely in finding complete sets of $r_i$th roots of unity for $r_i \in \{1,2,...,n\}$.  
The below algorithm assumes we are working with exact precision (or symbolically).  Issues surrounding round-off for arbitrarily large matrices... are way outside the scope. 

algorithm for (2) which implies (1)
Suppose that your multiset of eigenvalues / numbers on the unit circle correspond to a Permutation matrix.  Place them in diagonal matrix $D$ and focus on taking the trace for powers of $D$.  Also construct n x n diagonal matrix $\Lambda$ where $\lambda_i = 0$ for $i\in\{1,2,..,n\}$ at initialization (i.e. it is initialized to be the zero matrix).  
Now compute $\text{trace}\big(D\big) = \alpha$.  If $\alpha =0$ proceed to next step.  If $\alpha$ is not a natural number you can reject the hypothesis that $D$ is similar to a permutation matrix, and the algorithm terminates. If $\alpha\gt 0$ (i.e. $\alpha$ is a positive integer), then set the first available $\alpha$ components on the diagonal of $\Lambda$ equal to one, i.e. $\lambda_i :=1$ for $i\in\{1,2,...,\alpha\}$, and now proceed to the next step.  
Compute
$\text{trace}\big(D^k - \Lambda^k\big) = \text{trace}\big(D^k\big) - \text{trace}\big(\Lambda^k\big)$
for $k=2,3,...,n$.
It should be zero or a positive integer throughout -- if you see anything else, the algorithm terminates and we  reject the eigenvalues as belonging to any $\text{n x n}$ permutation matrix.  (This is essentially an assert line in a for loop.)   
As soon as you see some positive integer, you want to 'nullify it'.  (The idea to 'nullify' comes from the idea of reducibility of a graph in Perron Frobenius theory.)  Put differently, we have $\text{trace}\big(D^1 - \Lambda^1\big)=0$ and we want to replicate this behavior for higher powers.
Suppose e.g. $\text{trace}\big(D^2-\Lambda^2\big)=5$.  This implies $5$ 2nd roots of unity that did not 'show up' in $\text{trace}\big(D^1\big)$ -- this is a problem... because if they are complete sets we need $\text{5%2} =0$.  So if  $\text{trace}\big(D^2- \Lambda^2\big)=\alpha$ and $\alpha \text{% 2} \neq 0$ the algorithm terminates and we reject the eigenvalues as being associated with any $\text{n x n}$ permutation matrix $P$.  Otherwise if $\alpha \text{% 2} = 0$ and $\alpha$ is a positive integer, proceed to 'nullify' these roots of unity.  Say e.g. $\alpha = 6$.  Then place 3 complete sets of 2nd roots of unity, in 3 contiguous 'diagonal blocks' in the the next available slots of $\Lambda$ (i.e. first remaining zeros on the diagonal) .  
Now proceed to compute $\text{trace}\big(D^3-\Lambda^3\big)$ and repeat the above process, then $\text{trace}\big(D^4-\Lambda^4\big)$ and repeat above process, then $\text{trace}\big(D^5-\Lambda^5\big)$, and so on.  
After this has been done all the way up to and including $k=n$ without raising an error, then we know
$\text{trace}\big(D^k-\Lambda^k\big) = 0 \longrightarrow \text{trace}\big(D^k\big) = \text{trace}\big(\Lambda^k\big)$
for $k \in\{1,2,...,n\}$
By Newton's Identities, $D$ has the same characteristic polynomial (and hence eigenvalues) as $\Lambda$.  I.e. up to graph isomorphism $D = \Lambda$. 
The underlying idea here is that sums of sets of complete $r_i$th roots of unity ($r_i \geq 2)$ telescope for powers less than $r_i$ so each time we append new $r_i$th roots of unity to $\Lambda$ we don't impact the $\text{trace}\big(\Lambda^k\big)$ for $k\in \{1,2,...,r_i-1\}$, and this allows us to exploit a kind of optimal substructure.  
closing ideas:
Newton's Identities implies $\big \vert\det\big(\Lambda \big )\big \vert  = \big \vert\det\big(D\big )\big \vert = 1$ -- so if the algorithm made it through $k=n$ without raising an error, then $\Lambda$ has no zeros on its diagonal.  Furthermore, we constructed all the non-zero diagonal components of $\Lambda$ to come in diagonal 'blocks' $r_i$ at a time corresponding to complete sets of $r_i$th roots of unity.  Thus $\Lambda$ can be separated into $m$ distinct complete sets of $r_i$ roots of unity for $r_i \in \{1,2...,n\}$ for $i\in\{1,2,....,m\}$ where $m$ is the number of 1's in $\Lambda$ -- and e.g. Perron Frobenius Theory tells us the associated permutation matrix $P$ has $m$ distinct communicating classes, each of size $r_i$.  
This motivates a different approach to the algorithm -- rather than constructing $\Lambda$ as we go, we could construct a particularly nice $P$ as the algorithm proceeds.  e.g. if at step one we had $\alpha = 5$, then the next step had $\alpha = 4$, and then $\alpha = 3$, then after 3 iterations we'd have a candidate $P= \begin{bmatrix}
 I_{5} &  \mathbf 0 &  \mathbf 0& \mathbf 0 & \mathbf 0 &  \cdots &  \mathbf 0 \\ 
\mathbf 0&  C^{(2)} &  \mathbf 0& \mathbf 0 & \mathbf 0 &  \cdots &   \mathbf 0 \\ 
 \mathbf 0&   \mathbf 0&    C^{(2)}&\mathbf 0 & \mathbf 0 &  \cdots  &  \mathbf 0 \\ 
 \mathbf 0&   \mathbf 0& \mathbf 0&   C^{(3)} &  \mathbf 0 &\cdots  &  \mathbf 0 \\ 
 \mathbf 0&   \mathbf 0&   \mathbf 0& \mathbf 0& \mathbf 0 &  \cdots &  \mathbf 0 \\ 
 \vdots&   \vdots&   \vdots& \vdots& \vdots &  \ddots &  \vdots \\ 
 \mathbf 0&   \mathbf 0 &  \mathbf 0& \mathbf 0& \mathbf 0 &   \cdots&  \mathbf 0
\end{bmatrix}$ 
where $C^{(j)}$ is the $\text{j x j}$ 'cyclic permutation matrix' which is also the Companion matrix associated with the polynomial $x^j -1$.  Of course, re-visiting this entire process, it becomes clear that we never needed the diagonal matrix $D$ that our unitary matrix $M$ is similar to -- we could have simply worked with traces of powers of $M$.     
