  It is shown that in a large class of disordered systems with non-degenerate
disorder, in presence of non-local interactions, the Integrated Density of
States (IDS) is at least H\"older continuous in one dimension and universally
infinitely differentiable in higher dimensions. This result applies also to the
IDS in any finite volume subject to the random potential induced by an ambient,
infinitely extended disordered media. Dimension one is critical: in the
Bernoulli case, within the class of exponential interactions, the IDS measure
undergoes continuity phase transitions, from absolutely continuous to singular
continuous behaviour (the singularity in the latter case was known before). The
continuity transitions do not occur for sub-exponential or slower decaying
interactions, nor for $d\ge 2$. Technically, the case of polynomial decay is
the simplest one.
  The proposed approach provides a complement to the classical Wegner estimate
which says, essentially, that the IDS in the short-range models is at least as
regular as the marginal distribution of the disorder. In the models with
non-local interaction the IDS is actually much more regular than the underlying
disorder, which can even be discrete, due to the smoothing effect of multiple
convolutions. In turn, smoothness of the IDS is responsable for a mechanism
complementing the usual Lifshitz tails phenomenon.
  It is also shown that the disorder can take various forms (e.g., substitution
or random displacements) and need not be stochastically stationary (as in
Delone--Anderson or trim\-med/crooked Hamiltonians, for example); this does not
affect the main phenomena observed already in the simplest setting.
  Long-range models have an amazingly large number of connections to several
classical problems of harmonic analysis, probability theory, dynamical systems
and number theory.
