Random systems of hard spheres are important models or starting points for models of many structures in physics and materials science. Today, they can be explored only by means of simulations. The lecture reports on simulation results obtained by means of the so-called force-biased algorithm, which can be also used to obtain super-dense packings with densities larger than 0.64. Various point-process-statistical summary characteristics can be used to describe the increasing degree of order in packings of identical spheres with increasing density. They include the pair correlation function, the bond-orientational order characteristic and various characteristics related to the Voronoi tessellation with respect to the sphere centres. Random packings of spheres with variable diameters show many interesting phenomena. Some of them are discussed in the lecture. Additionally to the point-process characteristics also set-theoretic characteristics are discussed, in particuar lineal characteristics such as chord lengths. Simulations show that for the chord length and linear contact distributions the exponential distribution is an excellent approximation, while for the spherical contact distribution the positive part of the normal distribution can be recommended as approximation.