We compute analytically, for large N, the probability distribution of the number of positive eigenvalues (the index N+) of a random N x N matrix belonging to Gaussian orthogonal (beta=1), unitary (beta=2) or symplectic (beta=4) ensembles. The distribution of the fraction of positive eigenvalues c=N+/N scales, for large N, as P(c,N) approximately = exp[-betaN(2)Phi(c)] where the rate function Phi(c), symmetric around c=1/2 and universal (independent of beta), is calculated exactly. The distribution has non-Gaussian tails, but even near its peak at c=1/2 it is not strictly Gaussian due to an unusual logarithmic singularity in the rate function.