I'm a bit slow in understanding the computation of darktable from the code, so for now I played a bit in MATLAB. But I definitely want to challenge myself with real code.
I don't know if I'm missing something fundamental so I'll write down few notes. Criticate freely.
Adding noise using the response curve of the photographic paper
I started slightly modifing a logarithmic sigmoidal function f(x) which represents the typical shape of a response curve of photographic paper,
here x is the logarithm of the exposure of the photographic paper, gamma is the slope in the midtones, and delta is a small shift. The values of f(x) represent an hypothetical density impressed on the paper.
Here an example of response curve.
Let's assume that the grain is present uniformly on the photographic film. The light filtering trough the photographic film will be affected uniformly and we need to add the noise to the logarithm of the exposure.
In order to add the noise we must convert the lightness channel of the photo (that I assume proportional to the density impressed on the photographic paper) in the originating log-exposure values. We can use the inverse of f(x) defined as
Taking the lightness channel (L) defined as between 0 and 1, gamma equal to 1 and delta equal to 0.005, we have
Now we can add the grain-noise (n) to this function,
and then convert the noisy exposure into the noisy density on the photographic paper using the function f
We have obtained the grain with a midtones bias and highlights and shadows skewed noise distributions. The average values of lightness are unaffected.
The lightness channel (L) of the photograph should be edited accordingly to the following equation
where k is a parameter to tune the amount of grain. The procedure leaves unaffected the average picture but distort and stretches the noise adapting its characteristics to the specific lightness value.
A quick look at the grain
Using the final equation, we can apply some grain and see how it looks. At this stage I'm using Gaussian random noise:
To me it looks quite similar to the snapseed result.
The highligths/shadows skewed distributions are more evident when watching only the grain subtracting the linear trend:
On the left of the image there are more white dots and on the right of the image there are more black dots.
The statistical properties of the noise are compelling.
Example with Fourier filtered grain
We can Fourier filter the Gaussian noise in order to mimic the spatial properties of real grain and then apply the photographic response curve. In ref.  a simple model for the power spectrum in the spatial frequency frequency is used: a sum of two Gaussian function associated at two typical size of film grains.
Here an example using a real image and Fourier filtered grain.
@houz: With Fourier filtering is easy to achieve large blotches but I have some concerns about performances. Am I right?
 Stephenson I, Saunders A., Simulating Film Grain using the Noise-Power Spectrum, EG UK Theory and Practice of Computer Graphics (2007).