  We introduce an algorithm to decompose orthogonal matrix representations of
the symmetric group over the reals into irreducible representations, which as a
by-product also computes the multiplicities of the irreducible representations.
The algorithm applied to a $d$-dimensional representation of $S_n$ is shown to
have a complexity of $O(n^2 d^3)$ operations for determining which irreducible
representations are present and their corresponding multiplicities and a
further $O(n d^4)$ operations to fully decompose representations with
non-trivial multiplicities. These complexity bounds are pessimistic and in a
practical implementation using floating point arithmetic and exploiting
sparsity we observe better complexity. We demonstrate this algorithm on the
problem of computing multiplicities of two tensor products of irreducible
representations (the Kronecker coefficients problem) as well as higher order
tensor products. For hook and hook-like irreducible representations the
algorithm has polynomial complexity as $n$ increases. We also demonstrate an
application to constructing a basis of multivariate orthogonal polynomials with
respect to a tensor product weight so that applying a permutation of variables
induces an irreducible representation.
