We consider N × N Gaussian random matrices, whose average density of eigenvalues has the Wigner semicircle form over [-√2],√2]. For such matrices, using a Coulomb gas technique, we compute the large N behavior of the probability P(N,L)(N(L)) that N(L) eigenvalues lie within the box [-L,L]. This probability scales as P(N,L)(N(L) = κ(L)N) ≈ exp(-βN(2)ψ(L)(κ(L))), where β is the Dyson index of the ensemble and ψ(L)(κ(L)) is a β-independent rate function that we compute exactly. We identify three regimes as L is varied: (i) N(-1)≪L < √2 (bulk), (ii) L∼√2 on a scale of O(N(-2/3)) (edge), and (iii) L > sqrt[2] (tail). We find a dramatic nonmonotonic behavior of the number variance V(N)(L) as a function of L: after a logarithmic growth ∝ln(NL) in the bulk (when L∼O(1/N)), V(N)(L) decreases abruptly as L approaches the edge of the semicircle before it decays as a stretched exponential for L > sqrt[2]. This "dropoff" of V(N)(L) at the edge is described by a scaling function V(β) that smoothly interpolates between the bulk (i) and the tail (iii). For β = 2 we compute V(2) explicitly in terms of the Airy kernel. These analytical results, verified by numerical simulations, directly provide for β = 2 the full statistics of particle-number fluctuations at zero temperature of 1D spinless fermions in a harmonic trap.