  The flux of $\gamma$ rays is measured with unprecedented accuracy by the
$\textit{Fermi}$ Large Area Telescope from 100 MeV to almost 1 TeV. In the
future, the Cherenkov Telescope Array will have the capability to measure
photons up to 100 TeV. To accurately interpret this data, precise predictions
of the production processes, specifically the cross section for the production
of photons from the interaction of cosmic-ray protons and helium with atoms of
the ISM, are necessary. In this study, we determine new analytical functions
describing the Lorentz-invariant cross section for $\gamma$-ray production in
hadronic collisions. We utilize the limited total cross section data for
$\pi^0$ production channels and supplement this information by drawing on our
previous analyses of charged pion production to infer missing details. In this
context, we highlight the need for new data on $\pi^0$ production. Our
predictions include the cross sections for all production channels that
contribute down to the 0.5% level of the final cross section, namely $\eta$,
$K^+$, $K^-$, $K^0_S$, and $K^0_L$ mesons as well as $\Lambda$, $\Sigma$, and
$\Xi$ baryons. We determine the total differential cross section
$d\sigma(p+p\rightarrow \gamma+X)/dE_{\gamma}$ from 10 MeV to 100 TeV with an
uncertainty of 10% below 10 GeV of $\gamma$-ray energies, increasing to 20% at
the TeV energies. We provide numerical tables and a script for the community to
access our energy-differential cross sections, which are provided for incident
proton (nuclei) energies from 0.1 to $10^7$ GeV (GeV/n).
