  We present a lattice determination of the vector and scalar form factors of
the $D \to \pi(K) \ell \nu$ semileptonic decays, which are relevant for the
extraction of the CKM matrix elements $|V_{cd}|$ and $|V_{cs}|$ from
experimental data. Our analysis is based on the gauge configurations produced
by the European Twisted Mass Collaboration with $N_f = 2+1+1$ flavors of
dynamical quarks, at three different values of the lattice spacing and with
pion masses as small as 210 MeV. The matrix elements of both vector and scalar
currents are determined for a plenty of kinematical conditions in which parent
and child mesons are either moving or at rest. Lorentz symmetry breaking due to
hypercubic effects is clearly observed in the data and included in the
decomposition of the current matrix elements in terms of additional form
factors. After the extrapolations to the physical pion mass and to the
continuum limit we determine the vector and scalar form factors in the whole
kinematical region from $q^2 = 0$ up to $q^2_{max} = (M_D - M_{\pi(K)})^2$
accessible in the experiments, obtaining a good overall agreement with
experiments, except in the region at high values of $q^2$ where some deviations
are visible. A set of synthetic data points, representing our results for
$f_+^{D \pi(K)}(q^2)$ and $f_0^{D \pi(K)}(q^2)$ for several selected values of
$q^2$, is provided and also the corresponding covariance matrix is available.
At zero 4-momentum transfer we get: $f_+^{D \to \pi}(0) = 0.612 ~ (35)$ and
$f_+^{D \to K}(0) = 0.765 ~ (31)$. Using the experimental averages for
$|V_{cd}| f_+^{D \to \pi}(0)$ and $|V_{cs}| f_+^{D \to K}(0)$, we extract
$|V_{cd}| = 0.2330 ~ (137)$ and $|V_{cs}| = 0.945 ~ (38)$, respectively. The
second-row of the CKM matrix is found to be in agreement with unitarity within
the current uncertainties: $|V_{cd}|^2 + |V_{cs}|^2 + |V_{cb}|^2 = 0.949 ~
(78)$.
