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).