Is particle number a problem for formulating statistical physics in a mathematically rigorous manner? Quantities like the chemical potential can be expressed as something like 
$$\mu=-T\left(\tfrac{\partial S}{\partial N}\right)_{E,V}.$$
Now the entropy is the log some volume, which depends on the particle number $N$. As in this definition, we sum over natural numbers of particles, is there any good way of actually evaluating the derivative? 
What one practically does, i.e. when dealing with an ideal gas, is computing the quantity $\Gamma_N$, which might turn out to be $\frac{\pi^{N/2}}{(n/2)!}$, and then one will get an expression $S(E,N,V)$ which can of course be treated as if it was a function over $\mathbb{R}^3$. Even if that assumes that one has a closed expression which is a function $N$.
In principle I'd be fine with that - if one has a given function (Or at least the bunch of values for all $N$) over a grid and a procedure to introduce more and more grid points to get a finer mash, then there is a notion of convergence to a derivative. But here the N's are clearly always at least 1 value apart - no matter how many N s there are (thermodynamics limit), the mesh doesn't get finer between any two given points. 
You might define the derivative as computing the average rate of change between two partcle numbers $n$ and $n+d$ and say $\tfrac{\partial S}{\partial N}$ evaluated at $N'$ gives the same value for all $N'$ in one of the $d$-Intervals, but then you would have to postulate how to come up with $d$ in every new situation. This might be overcome in very specific situations in coming up with a "reasonable" fraction of Avogadros number, but this is not quite mathematical and the values of different finite difference approximation schemes are always different.
In full generality, I feel there is no categorical understanding of what the fractional dimensional space (phase space in this case) has to be and so the procedure of evaluation of the derivative should be explicitly postulate.
 A: No, it's not a problem. The reason is that, in order for expressions like
$$\mu=-T\left(\tfrac{\partial S}{\partial N}\right)_{E,V}.$$
to be meaningful, you have to be using the grand canonical ensemble (or a generalisation thereof), in which particles are able to enter and leave the system. Consequently, $N$ stands not for an integer number of particles, but for an expectation value: $N= \sum_i p_i N_i$, where the sum is over all possible states of the system, and $N_i$ is the number of particles when the system is in state $i$.  Since the probabilities $p_i$ are continuous, $N$ is continuous also, and there is no mathematical problem with considering an arbitrarily small change in it.
The way to do this is to note that, for the grand-canonical ensemble with a single chemical species,
$$
p_i = \frac{1}{Z(\mathbf{\lambda})} \exp\left( -\lambda_U U_i -\lambda_V V_i - \lambda_N N_i \right),
$$
where $\lambda_U$, often denoted $\beta$, is $1/T$; $\lambda_V = P/T$ and $\lambda_N = -\mu/T$, (I have suppressed Boltzmann's constant for readability), and $\mathbf{\lambda}$ indicates a vector of all three $\lambda$'s. This allows us to write the entropy as
$$
S = -\sum_i p_i \log p_i = \sum_i p_i (\log Z(\mathbf{\lambda}) + \lambda_U U_i + \lambda_V V_i + \lambda_N N_i)\\
= \log Z(\mathbf{\lambda}) + \lambda_U U + \lambda_V V + \lambda_N N,
$$
where again the $U$, $V$ and $N$ without subscripts indicate the expectation values. $S$ may now be considered a function of these expectation values, allowing us to write
$$
\frac{\partial S(U,V,N)}{\partial N} = \lambda_N = -\mu/T,
$$
which gives us your original expression.
Considering $S$ as a function of the extensive variables $U$, $V$ and $N$ rather than the intensive $\lambda$'s might seem a bit of an odd move, but it makes sense if you take the MaxEnt approach to statistical mechanics. In this formalism the above expression for $p_i$ is derived by maximising $-\sum_i p_i \log p_i$ subject to the constraints that the expectation values $U$, $V$ and $N$ must take particular values. In this case it is very natural to consider $S$ a function of their values, and $Z$ a function only of the $\lambda$'s.
A: The formula you write down is one of thermodynamics. In the statistical mechanics version it is valid in the grand canonical ensemble only if you interpret the extensive variables as expectation values. (See, e.g., Chapter 9 of my online book Classical and Quantum Mechanics via Lie algebras, arXiv:0810.1019.) But expectation values are continuous even when the corresponding operator has a discrete spectrum.
If you work with instead with the canonical or microcanonical ensemble, the formula is valid only in the thermodynamic limit. In this limit, the thermodynamic variable that makes sense in the limit is not $N$ but $\overline N = N\overline V/V$, where $\overline V$ is a finite reference volume that in the thermodynamic formulas takes the place of $V$. Again, $\overline N$ takes in the limit a continuum of values.
Thus in each of the versions typically covered in statistical mechanics, everything is consistent, as there is no theoretical problem in taking the derivative. Of course, finding closed formulas for the thermodynamic limit is the hard part, tractable only in a few cases, but in these cases one can verify differentiability except at points of phase transition. 
A: Gibbs' thought on this was (Elementary principles, page 204 footnote) "Strictly speaking, $\psi_{\rm gen}$ is not determined as function of $\nu_1,\ldots\nu_h$, except for integral values of these variables. Yet we may suppose it to be determined as a continuous function by any suitable process of interpolation." Here $\psi_{\rm gen}$ is the free energy of the canonical ensemble including correct Boltzmann counting for indistinguishable particles, and $\nu_i$ are numbers of different particle species.
However, "any suitable process" is indeed a bit vague. I could construct a hypothetical system where even particle numbers were vastly favoured over odd particle number (say, making an energy penalty for an odd particle number). If I tried to define chemical potential by a simple finite difference then I would get into trouble.
The answer (as noted by Nathaniel and Arnold) is that when particle number fluctuations are allowed, i.e. in grand ensembles, then the expectation value of particle number varies continuously and the partial derivative is good. With the grand canonical ensemble, in particular, Gibbs showed that we have a well defined and useful value of chemical potential that works all the way down to microscopic systems.
(Actually the discrete particle number problem is no different from the discrete energy value problem, in quantum systems. Again the canonical ensemble, either grand or petit, provides a perfectly functional value of temperature even in microscopic systems. To have a well defined temperature in a microscopic system, the price is that we must allow the ensemble to involve systems of a variety of energy values.)
