We evaluate, in the large-N limit, the complete probability distribution P(A,m) of the values A of the sum ∑i=1N|λi|m , where λ i ( i=1,2,…,N ) are the eigenvalues of a Gaussian random matrix, and m is a positive real number. Combining the Coulomb gas method with numerical simulations using a matrix variant of the Wang–Landau algorithm, we found that, in the limit of N→∞ , the rate function of P(A,m) exhibits phase transitions of different characters. The phase diagram of the system on the (A, m) plane is surprisingly rich, as it includes three regions: (i) a region with a single-interval support of the optimal spectrum of eigenvalues, (ii) a region emerging for m < 2 where the optimal spectrum splits into two separate intervals, and (iii) a region emerging for m > 2 where the maximum or minimum eigenvalue ‘evaporates’ from the rest of eigenvalues and dominates the statistics of A. The phase transition between regions (i) and (iii) is of second order. Analytical arguments and numerical simulations strongly suggest that the phase transition between regions (i) and (ii) is of (in general) fractional order p=1+1/|m−1| , where 0<m<2 . The transition becomes of infinite order in the special case of m = 1, where we provide a more complete analytical and numerical description. Remarkably, the transition between regions (i) and (ii) for m⩽1 and the transition between regions (i) and (iii) for m > 2 occur at the ground state of the Coulomb gas which corresponds to the Wigner’s semicircular distribution.