I have the following matrix $M$:
\begin{bmatrix}
P & p & 0 & 0 \\
pe\_1 & Pe\_1 & \sqrt{P^2 - p^2}q\_1 & \sqrt{P^2 - p^2}r\_1 \\
pe\_2 & Pe\_2 & \sqrt{P^2 - p^2}q\_2 & \sqrt{P^2 - p^2}r\_2 \\
pe\_3 & Pe\_3 & \sqrt{P^2 - p^2}q\_3 & \sqrt{P^2 - p^2}r\_3
\end{bmatrix}
with $P$, $p$ two real numbers such that $\lvert P\rvert > \lvert p\rvert$ and $(\mathbf{e}, \mathbf{q}, \mathbf{r})$ that forms an orthonormal triad in Euclidean space.

I want to invert the matrix $M$.

I tried the following, but got nothing:

```
Assuming[
  RealAbs[P] > RealAbs[p] && 
  (Indexed[e, 1])^2 + (Indexed[e, 2])^2 + (Indexed[e, 3])^2 == 1 && 
  (Indexed[q, 1])^2 + (Indexed[q, 2])^2 + (Indexed[q, 3])^2 == 1 && 
  (Indexed[r, 1])^2 + (Indexed[r, 2])^2 + (Indexed[r, 3])^2 == 1 && 
  Indexed[e, 1]*Indexed[q, 1] + Indexed[e, 2]*Indexed[q, 2] + 
    Indexed[e, 3]*Indexed[q, 3] == 0 && 
  Indexed[e, 1]*Indexed[r, 1] + Indexed[e, 2]*Indexed[r, 2] + 
    Indexed[e, 3]*Indexed[r, 3] == 0 && 
  Indexed[q, 1]*Indexed[r, 1] + Indexed[q, 2]*Indexed[r, 2] + 
    Indexed[q, 3]*Indexed[r, 3] == 0, 
  inverse[M]]

```

I also tried

```
PseudoInverse[M]

```

but I got a message saying that a very large output was generated and I'm unable to see the output.

Could you please explain how to invert $M$ with Mathematica?