  Determining the sound speed $c_s$ in compact stars is an important open
question with numerous implications on the behaviour of matter at large
densities and hence on gravitational-wave emission from neutron stars. To this
scope, we construct more than $10^7$ equations of state (EOSs) with continuous
sound speed and build more than $10^8$ nonrotating stellar models consistent
not only with nuclear theory and perturbative QCD, but also with astronomical
observations. In this way, we find that EOSs with sub-conformal sound speeds,
i.e. with $c^2_s < 1/3$ within the stars, are possible in principle but very
unlikely in practice, being only $0.03\%$ of our sample. Hence, it is natural
to expect that $c^2_s > 1/3$ somewhere in the stellar interior. Using our large
sample, we obtain estimates at $95\%$ credibility of neutron-star radii for
representative stars with $1.4$ and $2.0$ solar masses,
$R_{1.4}=12.42^{+0.52}_{-0.99}\,{\rm km}$, $R_{2.0}=12.12^{+1.11}_{-1.23}\,{\rm
km}$, and for the binary tidal deformability of the GW170817 event,
$\tilde\Lambda_{1.186}=485^{+225}_{-211}$. Interestingly, our lower-bounds on
the radii are in very good agreement with the prediction derived from very
different arguments, namely, the threshold mass. Finally, we provide simple
analytic expressions to determine the minimum and maximum values of
$\tilde\Lambda$ as a function of the chirp mass.
