diff --git a/CHANGELOG.md b/CHANGELOG.md
index 24e8911..58c9d3c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # News
 
+## v1.2.7 - 2025-01-25
+
+- Use `blochmessiah` decomposition from SymplecticFactorizations.jl.
+
 ## v1.2.6 - 2025-01-14
 
 - Add `changebasis` for changing symplectic bases of Gaussian types.
diff --git a/Project.toml b/Project.toml
index 21af7b0..0a10806 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Gabs"
 uuid = "0eb812ee-a11f-4f5e-b8d4-bb8a44f06f50"
 authors = ["Andrew Kille"]
-version = "1.2.6"
+version = "1.2.7"
 
 [deps]
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
@@ -23,7 +23,7 @@ LinearAlgebra = "1.9"
 Makie = "0.21, 0.22"
 QuantumInterface = "0.3.4"
 StaticArrays = "1.9.7"
-SymplecticFactorizations = "0.1.1"
+SymplecticFactorizations = "0.1.4"
 julia = "1.6.7, 1.10.0"
 
 [extras]
diff --git a/docs/src/manual.md b/docs/src/manual.md
index e05b79f..df1d81c 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -452,11 +452,14 @@ symplectic: 4×4 Matrix{Float64}:
   0.0       0.0       -0.479426  0.877583
 ```
 Various symplectic decompositions are supported in Gabs through the symplectic linear algebra package [SymplecticFactorizations.jl](https://github.com/apkille/SymplecticFactorizations.jl). Particularly important
-ones are the Williamson decomposition ([`williamson`](@ref)) and symplectic polar decomposition ([`polar`](@ref)):
+ones are the Williamson decomposition ([`williamson`](@ref)), Bloch-Messiah/Euler decomposition ([`blochmessiah`](@ref)), and the symplectic polar decomposition ([`polar`](@ref)):
 ```@docs; canonical = false
 williamson
 ```
 ```@docs; canonical = false
+blochmessiah
+```
+```@docs; canonical = false
 polar
 ```
 Let's see an example with the Williamson decomposition:
diff --git a/src/Gabs.jl b/src/Gabs.jl
index 6a7205e..dce9e90 100644
--- a/src/Gabs.jl
+++ b/src/Gabs.jl
@@ -7,8 +7,8 @@ using LinearAlgebra: I, det, mul!, diagm, diag, qr, eigvals
 
 import QuantumInterface: StateVector, AbstractOperator, apply!, tensor, ⊗, directsum, ⊕
 
-import SymplecticFactorizations: williamson, Williamson, polar, Polar, randsymplectic, symplecticform, issymplectic
-using SymplecticFactorizations: williamson, Williamson, polar, Polar, BlockForm, PairForm
+import SymplecticFactorizations: williamson, Williamson, polar, Polar, blochmessiah, BlochMessiah, randsymplectic, symplecticform, issymplectic
+using SymplecticFactorizations: williamson, Williamson, polar, Polar, blochmessiah, BlochMessiah, BlockForm, PairForm
 
 export 
     # types
@@ -29,7 +29,7 @@ export
     # symplectic form and checks
     symplecticform, issymplectic, isgaussian, sympspectrum,
     # factorizations
-    williamson, Williamson, polar, Polar,
+    williamson, Williamson, polar, Polar, blochmessiah, BlochMessiah,
     # metrics
     purity
 
diff --git a/src/factorizations.jl b/src/factorizations.jl
index 3851612..4322f9f 100644
--- a/src/factorizations.jl
+++ b/src/factorizations.jl
@@ -2,7 +2,7 @@
     williamson(state::GaussianState) -> Williamson
 
 Compute the williamson decomposition of the `covar` field
-of a `Gaussian state` object and return a `Williamson` object.
+of a `GaussianState` object and return a `Williamson` object.
 
 A symplectic matrix `S` and symplectic spectrum `spectrum` can be obtained
 via `F.S` and `F.spectrum`.
@@ -21,10 +21,10 @@ function williamson(x::GaussianState{<:QuadPairBasis,M,V}) where {M,V}
 end
 
 """
-    polar(state::GaussianState) -> Polar
+    polar(state::GaussianUnitary) -> Polar
 
 Compute the Polar decomposition of the `symplectic` field
-of a `Gaussian unitary` object and return a `Polar` object.
+of a `GaussianUnitary` object and return a `Polar` object.
 
 `O` and `P` can be obtained from the factorization `F` via `F.O` and `F.P`, such that `S = O * P`.
 For the symplectic polar decomposition case, `O` is an orthogonal symplectic matrix and `P` is a positive-definite
@@ -32,11 +32,24 @@ symmetric symplectic matrix.
 
 Iterating the decomposition produces the components `O` and `P`.
 """
-function polar(x::GaussianUnitary{<:QuadBlockBasis,D,S}) where {D,S}
+polar(x::GaussianUnitary{B,D,S}) where {B<:SymplecticBasis,D,S} = polar(x.symplectic)
+
+"""
+    blochmessiah(state::GaussianUnitary) -> BlochMessiah
+
+Compute the Bloch-Messiah/Euler decomposition of the `symplectic` field
+of a `GaussianUnitary` and return a `BlockMessiah` object.
+
+The orthogonal symplectic matrices `O` and `Q` as well as the singular values `values` can be obtained
+via `F.O`, `F.Q`, and `F.values`, respectively.
+
+Iterating the decomposition produces the components `O`, `values`, and `Q`, in that order.
+"""
+function blochmessiah(x::GaussianUnitary{<:QuadBlockBasis,D,S}) where {D,S}
     basis = x.basis
-    return polar(x.symplectic)
+    return blochmessiah(BlockForm(basis.nmodes), x.symplectic)
 end
-function polar(x::GaussianUnitary{<:QuadPairBasis,D,S}) where {D,S}
+function blochmessiah(x::GaussianUnitary{<:QuadPairBasis,D,S}) where {D,S}
     basis = x.basis
-    return polar(x.symplectic)
+    return blochmessiah(PairForm(basis.nmodes), x.symplectic)
 end
\ No newline at end of file
diff --git a/test/test_factorizations.jl b/test/test_factorizations.jl
index 0fb47ad..8dd6bb9 100644
--- a/test/test_factorizations.jl
+++ b/test/test_factorizations.jl
@@ -50,4 +50,33 @@
         @test isapprox(P_block, transpose(P_block), atol = 1e-5) && all(i > 0 for i in eigvals(P_block))
         @test isapprox(P_pair, transpose(P_pair), atol = 1e-5) && all(i > 0 for i in eigvals(P_pair))
     end
+
+    @testset "Bloch-Messiah/Euler" begin
+        nmodes = rand(1:5)
+        qpairbasis = QuadPairBasis(nmodes)
+        qblockbasis = QuadBlockBasis(nmodes)
+        
+        uni_pair = randunitary(qpairbasis)
+        S_pair = uni_pair.symplectic
+        uni_block = randunitary(qblockbasis)
+        S_block = uni_block.symplectic
+
+        F_pair = blochmessiah(uni_pair)
+        O_pair, vals_pair, Q_pair = F_pair
+        F_block = blochmessiah(uni_block)
+        O_block, vals_block, Q_block = F_block
+
+        @test F_block.O == O_block && F_block.values == vals_block && F_block.Q == Q_block
+        @test F_pair.O == O_pair && F_pair.values == vals_pair && F_pair.Q == Q_pair
+        @test issymplectic(qpairbasis, O_pair, atol = 1e-5) && issymplectic(qblockbasis, O_block, atol = 1e-5)
+        @test issymplectic(qpairbasis, Q_pair, atol = 1e-5) && issymplectic(qblockbasis, Q_block, atol = 1e-5)
+        @test isapprox(inv(O_pair), transpose(O_pair), atol = 1e-5) && isapprox(inv(O_block), transpose(O_block), atol = 1e-5)
+        @test isapprox(inv(Q_pair), transpose(Q_pair), atol = 1e-5) && isapprox(inv(Q_block), transpose(Q_block), atol = 1e-5)
+
+        @test all(x -> x > 0.0, vals_block) && all(x -> x > 0.0, vals_pair)
+        D_block = Diagonal(vcat(vals_block, vals_block .^ (-1)))
+        D_pair = Diagonal(collect(Iterators.flatten(zip(vals_pair, vals_pair .^ (-1)))))
+        @test isapprox(O_block * D_block * Q_block, S_block, atol = 1e-5)
+        @test isapprox(O_pair * D_pair * Q_pair, S_pair, atol = 1e-5)
+    end
 end
\ No newline at end of file