Obtaining the charge from the charge density using distribution theory In electrostatics, for several reasons, it seems that the correct way to understand the charge density $\rho$ isn't as a function $\rho : \mathbb{R}^3\to \mathbb{R}$, but rather as a distribution $\rho \in \mathcal{S}'(\mathbb{R}^3)$.
For example, denoting $\delta_{a}$ the Dirac delta centered at $a\in \mathbb{R}^3$ and defined by $(\delta_a,f)=f(a)$, a collection of point charges is then given by
$$\rho=\sum q_i \delta_{x_i}.$$
Here I'm using rigorous distribution theory, where tempered distributions are continuous linear functionals defined on the Schwartz space. So there isn't any $\rho(x)$, because $\rho$ isn't a function on $\mathbb{R}^3$, nor any $\delta(x-x_i)$, because $\delta_{x_i}$ also isn't a function on $\mathbb{R}^3$.
Now, the point is that since $\rho$ isn't a function on $\mathbb{R}^3$, differently than what is done in Physics, it cannot be integrated. Indeed,
$$Q=\int \rho(x)d^3x$$
is meaningless because $\rho$ isn't a function. Actually, even if one tried to say that this would be solved with a measure on the space of functions, this wouldn't help much. If that were the case, we would be able to integrate functions with distribution values, but $\rho$ isn't that. It is a single distribution generaling a function on $\mathbb{R}^3$.
The only thing we can do with $\rho$ is to apply to functions. Interestingly if we pick $f = 1$, we have $(\rho,f)=Q$ in the point charge case, however, this is not correct since $1\notin \mathcal{S}(\mathbb{R}^3)$ and hence $\rho$ can't act over it.
So my question is: given $\rho\in \mathcal{S}'(\mathbb{R}^3)$ the charge density considered as a tempered distribution, and considering we want to do things with mathematical rigor, how do we recover the charge $Q$ from $\rho$ since we can't integrate it?
 A: There's no need to go to distribution theory: as far as the total charge in any arbitrary volume goes,* you can just define the charge density to be a measure (or, more precisely, a signed measure) over the usual Lebesgue $\sigma$-algebra $\mathscr A$ of sets in 3D space.
Thus, you can see $q:\mathscr A\to \Bbb R$ as a measure which assigns to any Lebesgue-measurable set $S\in\mathscr A$ the charge
$$
S\mapsto q(S)=\int_S\mathrm dq.
$$
Here $q$ can be a Dirac measure (essentially, giving $1$ if the chosen point is in $S$, and zero otherwise), giving you point charges, or you can also model volumetric charge densities by asking that $q$ be absolutely continuous with respect to the usual Lebesgue measure $\mu$, in which case the charge density is the Radon-Nikodym derivative of $q$ with respect to $\mu$.
If you then want to extend the measure to a distribution, for whatever reason, it's perfectly easy to do - just use the functional that will take any function $f:\mathbb R^3\to\mathbb R$ to its integral with respect to $q$,
$$
f\mapsto \int f(\mathbf r) \:\mathrm dq(\mathbf r).
$$
However, just because the charge distribution can give rise to a distribution doesn't mean that that is what it really is at its most fundamental $-$ what it really is is a measure.

*As noted in the comments, the Lebesgue-measure approach does miss out on being able to describe singular multipolar charge distributions such as the point-dipole charge density $\rho(\mathbf r) = p\delta'(x)
\delta(y)\delta(z)$ and its relatives. Those do matter, and if you want to describe them, then you need to up the charge density from signed measure to distribution. However, the charge densities that you gain by doing that (or at least the multipolar ones) integrate to zero charge over any well-behaved volume, so they don't affect the question as posed. 
A: $\mathcal{S}(\mathbb{R}^3)$ does not include the function $f(\mathbf{r}) = 1$, but does contain a sequence of function which tends to it in the relevant sense. 
Explicitly, $\mathcal{S}(\mathbb{R}^3)$ includes the bump functions on $\mathbb{R}^3$ so I can construct a function, $f_n(\mathbf{r})$, in $\mathcal{S}(\mathbb{R}^3)$ such that.
\begin{align}
f_n(\mathbf{r}) = 1\quad &\forall\mathbf{r}\in\bar{B}(nl,0)\\
f_n(\mathbf{r})= 0\quad &\forall\mathbf{r}\not\in B\left((n+1)l,0\right)
\end{align}
where $B\left(nl,0\right)$ and $\bar{B}(nl,0)$ are respectively the open and closed balls of radius $nl$ centred on the origin. Clearly we can then write
$$
Q = \lim_{n\rightarrow\infty}(\rho,f_n) = \lim_{n\rightarrow\infty}\int\mathrm{d}^3\mathbf{r}\;\rho(\mathbf{r})f_n(\mathbf{r})
$$
A: There are two issues: local regularity and 
decay at infinity which are in duality via the Fourier transform. Here one needs to introduce distributions primarily in order to allow very bad local regularity as arises for instance in the idealized situation of a point charge or a surface charge distribution etc.
In 1D, take a unit positive charge and concentrate it at the origin. The object you need is then the delta function. But if you take two such things of opposite charge and take a suitable limit where you put them both at the origin you get things like dipoles which require the derivative in the sense of distributions of the delta function which is not a signed measure anymore as in Emilio's answer. You can of course play this game in higher dimension and generate multipoles etc.
Now your problem has to do with decay at infinity rather than local regularity.
You need a restricted class of distributions which can be fed the constant test function equal to one. Although you will not find that in most textbooks on distributions (especially the "math for physicists" ones) Laurent Schwartz considered lots of spaces other than $\mathcal{D},\mathcal{D}',\mathcal{S},\mathcal{S}'$. there are also $\mathcal{E},\mathcal{E}',\mathcal{O}_M,\mathcal{O}_M',\mathcal{O}_C,\mathcal{O}_C',\ldots$
The one you need here, I think, is $\mathcal{O}_C'$ also called the space of rapidly decaying distributions. 
