1 files changed, 96 insertions, 0 deletions
diff --git a/python_originals/bloch_messiah.py b/python_originals/bloch_messiah.py
new file mode 100644
index 0000000..a4cc310
--- /dev/null
+++ b/python_originals/bloch_messiah.py
@@ -0,0 +1,96 @@
+def bloch_messiah(S, tol=1e-10, rounding=9):
+    r"""Bloch-Messiah decomposition of a symplectic matrix.
+
+    See :ref:`bloch_messiah`.
+
+    Decomposes a symplectic matrix into two symplectic unitaries and squeezing transformation.
+    It automatically sorts the squeezers so that they respect the canonical symplectic form.
+
+    Note that it is assumed that the symplectic form is
+
+    .. math:: \Omega = \begin{bmatrix}0&I\\-I&0\end{bmatrix}
+
+    where :math:`I` is the identity matrix and :math:`0` is the zero matrix.
+
+    As in the Takagi decomposition, the singular values of N are considered
+    equal if they are equal after np.round(values, rounding).
+
+    If S is a passive transformation, then return the S as the first passive
+    transformation, and set the the squeezing and second unitary matrices to
+    identity. This choice is not unique.
+
+    For more info see:
+    https://math.stackexchange.com/questions/1886038/finding-euler-decomposition-of-a-symplectic-matrix
+
+    Args:
+        S (array[float]): symplectic matrix
+        tol (float): the tolerance used when checking if the matrix is symplectic:
+            :math:`|S^T\Omega S-\Omega| \leq tol`
+        rounding (int): the number of decimal places to use when rounding the singular values
+
+    Returns:
+        tuple[array]: Returns the tuple ``(ut1, st1, vt1)``. ``ut1`` and ``vt1`` are symplectic orthogonal,
+            and ``st1`` is diagonal and of the form :math:`= \text{diag}(s1,\dots,s_n, 1/s_1,\dots,1/s_n)`
+            such that :math:`S = ut1  st1  v1`
+    """
+    (n, m) = S.shape
+
+    if n != m:
+        raise ValueError("The input matrix is not square")
+    if n % 2 != 0:
+        raise ValueError("The input matrix must have an even number of rows/columns")
+
+    n = n // 2
+    omega = sympmat(n)
+    if np.linalg.norm(np.transpose(S) @ omega @ S - omega) >= tol:
+        raise ValueError("The input matrix is not symplectic")
+
+    if np.linalg.norm(np.transpose(S) @ S - np.eye(2 * n)) >= tol:
+
+        u, sigma = polar(S, side="left")
+        ss, uss = takagi(sigma, tol=tol, rounding=rounding)
+
+        # Apply a permutation matrix so that the squeezers appear in the order
+        # s_1,...,s_n, 1/s_1,...1/s_n
+        perm = np.array(list(range(0, n)) + list(reversed(range(n, 2 * n))))
+
+        pmat = np.identity(2 * n)[perm, :]
+        ut = uss @ pmat
+
+        # Apply a second permutation matrix to permute s
+        # (and their corresonding inverses) to get the canonical symplectic form
+        qomega = np.transpose(ut) @ (omega) @ ut
+        st = pmat @ np.diag(ss) @ pmat
+
+        # Identifying degenerate subspaces
+        result = []
+        for _k, g in groupby(np.round(np.diag(st), rounding)[:n]):
+            result.append(list(g))
+
+        stop_is = list(np.cumsum([len(res) for res in result]))
+        start_is = [0] + stop_is[:-1]
+
+        # Rotation matrices (not permutations) based on svd.
+        # See Appendix B2 of Serafini's book for more details.
+        u_list, v_list = [], []
+
+        for start_i, stop_i in zip(start_is, stop_is):
+            x = qomega[start_i:stop_i, n + start_i : n + stop_i].real
+            u_svd, _s_svd, v_svd = np.linalg.svd(x)
+            u_list = u_list + [u_svd]
+            v_list = v_list + [v_svd.T]
+
+        pmat1 = block_diag(*(u_list + v_list))
+
+        st1 = pmat1.T @ pmat @ np.diag(ss) @ pmat @ pmat1
+        ut1 = uss @ pmat @ pmat1
+        v1 = np.transpose(ut1) @ u
+
+    else:
+        ut1 = S
+        st1 = np.eye(2 * n)
+        v1 = np.eye(2 * n)
+
+    return ut1.real, st1.real, v1.real
+
+