We study the statistics of the condition number $κ=λ_{\mathrm{max}}/λ_{\mathrm{min}}$ (the ratio between largest and smallest squared singular values) of $N\times M$ Gaussian random matrices. Using a Coulomb fluid technique, we derive analytically and for large $N$ the cumulative $\mathcal{P}[κ<x]$ and tail-cumulative $\mathcal{P}[κ>x]$ distributions of $κ$. We find that these distributions decay as $\mathcal{P}[κ<x]\approx\exp\left(-βN^2 Φ_{-}(x)\right)$ and $\mathcal{P}[κ>x]\approx\exp\left(-βN Φ_{+}(x)\right)$, where $β$ is the Dyson index of the ensemble. The left and right rate functions $Φ_{\pm}(x)$ are independent of $β$ and calculated exactly for any choice of the rectangularity parameter $α=M/N-1>0$. Interestingly, they show a weak non-analytic behavior at their minimum $\langleκ\rangle$ (corresponding to the average condition number), a direct consequence of a phase transition in the associated Coulomb fluid problem. Matching the behavior of the rate functions around $\langleκ\rangle$, we determine exactly the scale of typical fluctuations $\sim\mathcal{O}(N^{-2/3})$ and the tails of the limiting distribution of $κ$. The analytical results are in excellent agreement with numerical simulations.