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ölder 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 ≥ 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 trimmed/crooked Hamiltonians, for example); this does not affect the main phenomena observed already in the simplest setting.Contrary to the situation with the usual lattice Bernoulli-Anderson Hamiltonians, the proof of Anderson localization in the models with infinite-range interaction follows in a fairly standard way from the main bounds on the finite-volume IDS; taking into account a considerable size of the current text, the proof of localization is now presented in a separate work. Another distinction from the approach developed by Bourgain and Kenig for the continuous Bernoulli-Anderson Hamiltonians, and later extended by Germinet and Klein to arbitrary (locally IID) disorder, is that all nontrivial marginal distributions are treated in a unified way, via harmonic analysis and without reduction to an embedded Bernoulli model, thus keeping potential benefits from less singular forms of underlying disorder.Long-range models have an amazingly large number of connections to several classical problems of harmonic analysis, probability theory, dynamical systems and number theory.