The zero-temperature, classical XY model on an L×L square lattice is studied by exploring the distribution Φ_{L}(y) of its centered and normalized magnetization y in the large-L limit. An integral representation of the cumulant generating function, known from earlier works, is used for the numerical evaluation of Φ_{L}(y), and the limit distribution Φ_{L→∞}(y)=Φ_{0}(y) is obtained with high precision. The two leading finite-size corrections Φ_{L}(y)-Φ_{0}(y)≈a_{1}(L)Φ_{1}(y)+a_{2}(L)Φ_{2}(y) are also extracted both from numerics and from analytic calculations. We find that the amplitude a_{1}(L) scales as ln(L/L_{0})/L^{2} and the shape correction function Φ_{1}(y) can be expressed through the low-order derivatives of the limit distribution, Φ_{1}(y)=[yΦ_{0}(y)+Φ_{0}^{'}(y)]^{'}. Thus, Φ_{1}(y) carries the same universal features as the limit distribution and can be used for consistency checks of universality claims based on finite-size systems. The second finite-size correction has an amplitude a_{2}(L)∝1/L^{2} and one finds that a_{2}Φ_{2}(y)≪a_{1}Φ_{1}(y) already for small system size (L>10). We illustrate the feasibility of observing the calculated finite-size corrections by performing simulations of the XY model at low temperatures, including T=0.