We use recent progress in simulating the production of magnetohydrodynamic jets around black holes to derive the cosmic spin history of the most massive black holes, with masses >~10^8 Msol. Assuming the jet efficiency depends on spin a, we can approximately reproduce the observed `radio loudness of quasars and the local radio luminosity function. Using the X-ray luminosity function and the local mass function of supermassive black holes, SMBHs we can reproduce the individual radio luminosity functions of radio sources showing high- and low-excitation narrow emission lines. The data favour spin distributions that are bimodal, with one component around spin zero and the other close to maximal spin. In the low-excitation galaxies, the two components have similar amplitudes. For the high-excitation galaxies, the amplitude of the high-spin peak is typically much smaller than that of the low-spin peak. A bimodality should be seen in the radio loudness of quasars. We predict that the low-excitation galaxies are dominated by SMBHs with masses >~10^8 Msol, down to radio luminosity densities ~10^21 W Hz-1 sr-1 at 1.4~GHz. Our model is also able to predict the radio luminosity function at z=1, and predicts it to be dominated by high-excitation galaxies above luminosity densities >~10^26 W Hz-1 sr-1, in full agreement with the observations. From our parametrisation and using the best fitting jet efficiencies there is marginal evidence for evolution in spin: the mean spin increases slightly from <a>~0.25 at z=1 to <a>~0.35 at z=0, and the fraction of SMBHs with a>=0.5 increases from 0.16+-0.03 at z=1 to 0.24+-0.09 at z=0. Our results are in excellent agreement with the mean radiative efficiency of quasars, as well as recent cosmological simulations. We discuss the implications in terms of accretion and SMBH mergers, and galactic black holes (Abridged).