Finding Euler decomposition of a symplectic matrix

A symplectic matrix is a $2n\times2n$ matrix $S$ with real entries that satisfies the condition

$$ S^T \Omega S = \Omega $$ where $\Omega$ is the symplectic form, typically chosen to be $\Omega=\left(\begin{smallmatrix}0 & I_N \\ -I_N & 0\end{smallmatrix}\right)$. Sympletic matrices form the symplectic group $Sp(2n,\mathbb{R})$. Any symplectic matrix S can be decomposed as a product of three matrices as

\begin{equation} S = O\begin{pmatrix}D & 0 \\ 0 & D^{-1}\end{pmatrix}O' \quad \quad \forall S \in Sp(2n,\mathbb{R}), \end{equation} where $O, O'$ are orthogonal and symplectic - $\operatorname{Sp}(2n,\mathbb{R})\cap \operatorname{O}(2n)$; $D$ is positive definite and diagonal. The form of a matrix that is both symplectic and orthogonal can be given in block form as $O=\left(\begin{smallmatrix}X & Y \\ -Y & X\end{smallmatrix}\right)$, where $XX^T+YY^T=I_N$ and $XY^T-YX^T=0$. The decomposition above is known as Euler decomposition or alternatively as Bloch-Messiah decomposition.

How can I find the matrices in the decomposition for a given symplectic matrix?

Apparently, the decomposition is closely related to the singular value decomposition and I think the elements of the matrices $D$ and $D^{-1}$ coincide with the singular values of $S$. Also, I have the impression that the case where it can be assumed that $S$ is also symmetric is easier. Any help, tips or pointers would be much appreciated!

Let $\varDelta$ stand for your diagonal matrix to write the factorization as $S=O\varDelta O'$. Now, rewrite a bit more: $$ S=(O\varDelta O^T)(OO')=\varSigma U, $$ with $\varSigma$ symplectic positive definite and $U$ symplectic and orthogonal. Thus we are looking for a polar decomposition. In fact such a decomposition is unique, and given by $$ \varSigma=(SS^T)^{1/2},\quad U=(SS^T)^{-1/2}S. $$ Here we use that a positive semidefinite (symplectic) matrix has a unique positive definite (symplectic) square root. Now, once we have $\varSigma$, it can be diagonalized with a symplectic orthogonal linear change $O$ (again quite constructive), and we get the data we sought. For the "basic facts" one can look at Gosson's book, for instance.