We study the probability density function for the fluctuations of the magnetic order parameter in the low-temperature phase of the XY model of finite size. In two dimensions, this system is critical over the whole of the low-temperature phase. It is shown analytically and without recourse to the scaling hypothesis that, in this case, the distribution is non-Gaussian and of universal form, independent of both system size and critical exponent eta. An exact expression for the generating function of the distribution is obtained, which is transformed and compared with numerical data from high-resolution molecular dynamics and Monte Carlo simulations. The asymptotes of the distribution are calculated and found to be of exponential and double exponential form. The calculated distribution is fitted to three standard functions: a generalization of Gumbel's first asymptote distribution from the theory of extremal statistics, a generalized log-normal distribution, and a chi(2) distribution. The calculation is extended to general dimension and an exponential tail is found in all dimensions less than 4, despite the fact that critical fluctuations are limited to D=2. These results are discussed in the light of similar behavior observed in models of interface growth and for dissipative systems driven into a nonequilibrium steady state.