The authors compute the next-to-next-to-leading order (NNLO) contributions to the splitting functions governing the evolution of the unpolarized flavour-singlet parton densities in perturbative QCD. The exact expressions are presented in both Mellin-N and Bjorken-x space, along with accurate parametrizations for practical applications. Their results agree with all partial results available in the literature. Unlike the non-singlet case, the leading logarithmic predictions for small x do not provide good estimates of the complete splitting functions. The size of the corrections and the stability of the NNLO evolution under variation of the renormalization scale are investigated. The perturbative expansion converges rapidly at x ≳ 10⁻³, but relatively large third-order corrections are found at smaller x. Parton distributions are essential for analyzing hard-scattering processes involving initial-state hadrons. Their scale dependence (evolution) is derived from first principles in terms of an expansion in powers of the strong coupling constant αs. The n-th order coefficients governing the evolution are the n-loop anomalous dimensions or splitting functions. Parton distributions evolved up to order αsⁿ⁺¹ in this expansion, together with the corresponding results for the partonic cross sections, constitute the NⁿLO approximation of perturbative QCD. The NNLO corrections are needed for quantitatively reliable predictions for hard processes at present and future high-energy colliders. These corrections are known for structure functions in deep-inelastic scattering (DIS) and for Drell-Yan lepton-pair and gauge-boson production in proton-(anti-)proton collisions. Work on NNLO cross sections for jet production is under way. For the three-loop splitting functions, only partial results had been computed until recently, especially the lowest six/seven Mellin moments and the leading (ln x)/x small-x terms of three of the four singlet splitting functions. These results have been used to improve the analysis of DIS data and hadron-collider predictions. However, this information is not sufficient for quantitative predictions at small x. The authors recently published the non-singlet part of the unpolarized three-loop splitting functions. In this article, they compute the corresponding singlet quantities. The article is organized as follows: In section 2, notations and the method are discussed. In section 3, the Mellin-N space results are presented. In section 4, the exact results and parametrizations for the x-space splitting functions are given, and their behavior at small x is studied. In section 5, the numerical implications of the results for the scale dependence of the singlet-quark and gluon distributions are illustrated. In section 6, the findings are summarized.The authors compute the next-to-next-to-leading order (NNLO) contributions to the splitting functions governing the evolution of the unpolarized flavour-singlet parton densities in perturbative QCD. The exact expressions are presented in both Mellin-N and Bjorken-x space, along with accurate parametrizations for practical applications. Their results agree with all partial results available in the literature. Unlike the non-singlet case, the leading logarithmic predictions for small x do not provide good estimates of the complete splitting functions. The size of the corrections and the stability of the NNLO evolution under variation of the renormalization scale are investigated. The perturbative expansion converges rapidly at x ≳ 10⁻³, but relatively large third-order corrections are found at smaller x. Parton distributions are essential for analyzing hard-scattering processes involving initial-state hadrons. Their scale dependence (evolution) is derived from first principles in terms of an expansion in powers of the strong coupling constant αs. The n-th order coefficients governing the evolution are the n-loop anomalous dimensions or splitting functions. Parton distributions evolved up to order αsⁿ⁺¹ in this expansion, together with the corresponding results for the partonic cross sections, constitute the NⁿLO approximation of perturbative QCD. The NNLO corrections are needed for quantitatively reliable predictions for hard processes at present and future high-energy colliders. These corrections are known for structure functions in deep-inelastic scattering (DIS) and for Drell-Yan lepton-pair and gauge-boson production in proton-(anti-)proton collisions. Work on NNLO cross sections for jet production is under way. For the three-loop splitting functions, only partial results had been computed until recently, especially the lowest six/seven Mellin moments and the leading (ln x)/x small-x terms of three of the four singlet splitting functions. These results have been used to improve the analysis of DIS data and hadron-collider predictions. However, this information is not sufficient for quantitative predictions at small x. The authors recently published the non-singlet part of the unpolarized three-loop splitting functions. In this article, they compute the corresponding singlet quantities. The article is organized as follows: In section 2, notations and the method are discussed. In section 3, the Mellin-N space results are presented. In section 4, the exact results and parametrizations for the x-space splitting functions are given, and their behavior at small x is studied. In section 5, the numerical implications of the results for the scale dependence of the singlet-quark and gluon distributions are illustrated. In section 6, the findings are summarized.