I) Right, the differential form of Gauss's law
$$\tag{1} {\bf\nabla} \cdot{\bf E}~=~ \frac{\rho}{\varepsilon_0} $$
uses the relatively advanced mathematical concept of Dirac delta distributions in case of point charges
$$\tag{2} \rho({\bf r})~=~\sum_{i=1}^n q_i\delta^3({\bf r}-{\bf r}_i).$$
Note in particular, that it is technically wrong to claim (as OP seems to do) that the Dirac delta distribution $\delta^3({\bf r})$ is merely a function $f:\mathbb{R}^3\to [0,\infty]$ that takes the value zero everywhere except at the origin where the value is infinity:
$$\tag{3} f({\bf r})~:=~\left\{ \begin{array}{rcl} \infty& {\rm for}& {\bf r}={\bf 0}, \cr 0& {\rm for}& {\bf r}\neq {\bf 0}.\end{array}\right. $$
For starters, for an arbitrary test function $g:\mathbb{R}^3\to [0,\infty[$, the Lebesgue integral$^1$
$$\tag{4} \int_{\mathbb{R}^3} \! d^3r~f({\bf r})g({\bf r}) ~=~0 $$
vanishes, in contrast to the defining property of the Dirac delta distribution
$$\tag{5} \int_{\mathbb{R}^3} \! d^3r~\delta^3({\bf r})g({\bf r}) ~=~g({\bf 0}). $$
The Dirac delta distribution $\delta^3({\bf r})$ is not a function. It is instead a generalized function. It is possible to give a mathematically consistent treatment of the Dirac delta distribution. However, it should be stressed that the analysis does not reduce to the investigation of two separate cases ${\bf r}= {\bf 0}$ and ${\bf r}\neq {\bf 0}$, but instead (typically) involves (smeared) test functions. To get a flavor of the various intricacies that can arise with distributions, the reader might find thisthis Phys.SE post interesting.
II) To avoid the notion of distributions, it is more safe (and probably more intuitive) to work with the equivalentequivalent integral form of Gauss's law
$$\tag{6} \Phi_{\bf E}~=~ \frac{Q_e}{\varepsilon_0}. $$
The corresponding Gauss's law for magnetism
$$\tag{7} \Phi_{\bf{B}}~=~ 0 $$
expresses (without employing double standards) the fact that there is no magnetic charge $Q_m$.
--
$^1$ Eq. (4) relies crucially on the fact that in integration theory for non-negative functions, one defines multiplication $\cdot: [0,\infty]\times[0,\infty]\to[0,\infty]$ on the extended real halfline $[0,\infty]$ so that $0\cdot\infty:=0$. Eq. (4) is essentially caused by the fact that $f$ is zero almost everywhere. Also we should mention the well-known fact that integration theory can be appropriately generalized from non-negative functions to complex-valued functions.