  Parton distribution functions (PDFs) are an essential ingredient for
theoretical predictions at colliders. Since their exact form is unknown, their
handling and delivery for practical applications relies on approximate
numerical methods. We discuss the implementation of PDFs based on a global
interpolation in terms of Chebyshev polynomials. We demonstrate that this
allows for significantly higher numerical accuracy at lower computational cost
compared with local interpolation methods such as splines. Whilst the numerical
inaccuracy of currently used local methods can become a nontrivial limitation
in high-precision applications, in our approach it is negligible for practical
purposes. This holds in particular for differentiation and for Mellin
convolution with kernels that have end point singularities. We illustrate our
approach for these and other important numerical operations, including DGLAP
evolution, and find that they are performed accurately and fast. Our results
are implemented in the C++ library ChiliPDF.
