The complex dielectric function of hexagonal gallium nitride (α-GaN) is obtained from a numerical solution of the excitonic Schrödinger equation taking into account the full 6 × 6 valence-band structure. The valence-band parametrization includes anisotropy, nonparabolicity, and the coupling of angular-momentum eigenstates. Spectra of excitonic eigenfunctions are obtained from a time-dependent initial-value problem, which is solved via an exponential split-operator method. In particular, we calculate the dielectric function and the reflectivity of a-plane GaN with polarization vectors parallel and perpendicular to the c-axis of the crystal. The simulated reflection spectra are in excellent agreement with recent experimental data and allow the unambiguous identification of the experimentally observed excitonic resonances. The binding energies of the FXA, FXB, and FXC excitons found in our calculation differ by up to 27%, depending on the chosen parameter set. An important consequence of this observation is that the experimentally observed splittings of the excitons cannot be used for the parametrization of the valence band near the Γ-point, but need to be corrected by the differences of the binding energies. This is of general relevance for all spectroscopic measurements in semiconductors with a wide bandgap.