If I have understood the code correctly (ignoring the "$-1$" in the computation of $m$), its input is an $n$-vector $\mu = (\mu_1, \ldots, \mu_n)$ and a symmetric $n$ by $n$ matrix $\Sigma = (\sigma_{ij})$. The output is an $n$-vector $m$ with
$$m_i = \exp(\mu_i + \sigma_{ii}/2)$$
and an $n$ by $n$ matrix $S$ with
$$S_{ij} = \exp(\mu_i + \mu_j + (\sigma_{ii}+\sigma_{jj})/2)(\exp(\sigma_{ij})-1) = m_i(\exp(\sigma_{ij})-1)m_j.$$
If this is correct, then we can solve readily for $\mu$ and $\Sigma$ in terms of $m$ and $S$ essentially by reversing these operations. Begin by forming the diagonal matrix $M$ whose diagonal entries are $1/m_i$: that is, $M_{ii}=1/m_i$ and $M_{ij}=0$ for $i\ne j$. From the right hand side of the preceding formula it follows immediately that
$$M S M + 1_n = \exp(\sigma_{ij})$$
and we easily recover $\Sigma$ by taking the logarithms term-by-term. With these values in hand,
$$\mu_i = \log(m_i) - \sigma_{ii}/2.$$
Edit
The code in the question uses "linear returns" rather than means. There's no problem with that: starting with the "returns" $m_i$ computed as $\exp(\mu_i + \sigma_{ii}/2)-1$, first add back the $1$ and proceed as above.
Since I haven't found a duplicate and nobody else has pointed one out, I have expanded my comment into an answer and added an example to show you how it works.
Any invertible matrix $A$ for which $AA^T=S$ will work by taking $C=A^{−1}$ .
So for example, Cholesky decomposition is one common way to do it - it's a simple algorithm for finding a lower triangular $A$ (or equivalently, in many implementations, an upper triangular $A^T$). That is, it finds $S=LL^T$, whence $L^{−1}S(L^{−1})^T=I$.
Example:
Consider the matrix
$$S=\left(\begin{matrix} 4 & 1 \\ 1 & 2.5 \end{matrix}\right) $$
We write $S = LL^T$, which implies that
$$s_{11} = l_{11}^2$$
That is, $l_{11}=2$ (we'll take positive square roots each time to make our diagonals on $L$ positive).
It also implies that $s_{21} = l_{21}l_{11}$, so $l_{21} = s_{21}/l_{11}$, i.e. $l_{21} = 1/2$. Then we see that $s_{22} = l_{21}^2+l_{22}^2$, so $l_{22}^2 = s_{22} - l_{21}^2 = 2.5 - 0.5^2 = 2.25$, in which case $l_{22} = 1.5$, so
$$L=\left(\begin{matrix} 2 & 0 \\ 0.5 & 1.5 \end{matrix}\right) $$
You can confirm that
$$C = L^{-1}=\left(\begin{matrix} \frac{1}{2} & 0 \\ -\frac{1}{6} & \frac{2}{3} \end{matrix}\right) $$
and that $CSC^T = I$.
In practice we don't invert $L$, however, but solve a linear system, it's more stable.
Another possibility is to use Singular Value Decomposition to find a matrix to do the diagonalizing, and then scale by the singular values.
But a number of other choices are available. (As whuber points out, a very large number.)
In practice, frequently we know $S$ is of the form $X^TX$ for some $X$. In that case, it's usually better to work on $X$ than $S$, via the QR decomposition.
Best Answer
There is no single number that encompasses all of the covariance information - there are 6 pieces of information, so you'd always need 6 numbers.
However there are a number of things you could consider doing.
Firstly, the error (variance) in any particular direction $i$, is given by
$\sigma_i^2 = \mathbf{e}_i ^ \top \Sigma \mathbf{e}_i$
Where $\mathbf{e}_i$ is the unit vector in the direction of interest.
Now if you look at this for your three basic coordinates $(x,y,z)$ then you can see that:
$\sigma_x^2 = \left[\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}\right]^\top \left[\begin{matrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{xz} & \sigma_{yz} & \sigma_{zz} \end{matrix}\right] \left[\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}\right] = \sigma_{xx}$
$\sigma_y^2 = \sigma_{yy}$
$\sigma_z^2 = \sigma_{zz}$
So the error in each of the directions considered separately is given by the diagonal of the covariance matrix. This makes sense intuitively - if I am only considering one direction, then changing just the correlation should make no difference.
You are correct in noting that simply stating:
$x = \mu_x \pm \sigma_x$
$y = \mu_x \pm \sigma_y$
$z = \mu_z \pm \sigma_z$
Does not imply any correlation between those three statement - each statement on its own is perfectly correct, but taken together some information (correlation) has been dropped.
If you will be taking many measurements each with the same error correlation (supposing that this comes from the measurement equipment) then one elegant possibility is to rotate your coordinates so as to diagonalise your covariance matrix. Then you can present errors in each of those directions separately since they will now be uncorrelated.
As to taking the "vector error" by adding in quadrature I'm not sure I understand what you are saying. These three errors are errors in different quantities - they don't cancel each other out and so I don't see how you can add them together. Do you mean error in the distance?