In a comment to Recursive randomness of integers Rahul pointed out a continuous version of the same problem: pick uniformly in
, then
uniformly in
, then
uniformly in
, etc. What is the distribution of the sum
?
The continuous version turns out to be easier to analyse. To begin with, it’s equivalent to picking uniformly distributed, independent and letting
. Then the sum is
which can be written as
So, is a stationary point of the random process
where
is uniformly distributed in
. Simply put,
and
have the same distribution. This yields the value of
in a much simpler way than in the previous post:
hence .
We also get an equation for the cumulative distribution function . Indeed,
The latter probability is . Conclusion:
. Differentiate to get an equation for the probability density function
, namely
. It’s convenient to change the variable of integration to
, which leads to
Another differentiation turns the integral equation into a delay differential equation,
Looks pretty simple, doesn’t it? Since the density is zero for negative arguments, it is constant on . This constant, which I’ll denote
, is
, or simply
. I couldn’t get an analytic formula for
. My attempt was
where
are the moments of
. The moments can be computed recursively using
, which yields
The first few moments, starting with , are 1, 1, 3/2, 17/6, 19/3, 81/5, 8351/80… Unfortunately the series
diverges, so this approach seems doomed. Numerically
which is not far from the Euler-Mascheroni constant, hence the choice of notation.
On the interval (1,2) we have , hence
for
.
The DDE gets harder to integrate after that… on the interval the solution already involves the dilogarithm (Spence’s function):
following SciPy’s convention for Spence. This is as far as I went… but here is an experimental confirmation of the formulas obtained so far (link to full size image).

To generate a sample from distribution S, I begin with a bunch of zeros and repeat “add 1, multiply by U[0,1]” many times. That’s it.
import numpy as np import matplotlib.pyplot as plt trials = 10000000 terms = 10000 x = np.zeros(shape=(trials,)) for _ in range(terms): np.multiply(x+1, np.random.uniform(size=(trials,)), out=x) _ = plt.hist(x, bins=1000, normed=True) plt.show()
I still want to know the exact value of … after all, it’s also the probability that the sum of our random decreasing sequence is less than 1.
Update
The constant I called “” is in fact
where
is indeed Euler’s constant… This is what I learned from the Inverse Symbolic Calculator after solving the DDE (with initial value 1) numerically, and calculating the integral of the solution. From there, it did not take long to find that
- The p.d.f. of S is the Dickman function, normalized to have integral 1. Wikipedia plots it on log scale, but Wolfram Mathworld has a plot identical to the above. Neither one mention the sum of random decreasing sequence.
- This exact problem was discussed on Mathematics Stack Exchange in February 2017.
- The problem and its solution go back to 1973 paper “A probabilistic approach to differential-difference equations arising in analytic number theory” by J.-M-F. Chamayou, cited in the survey Euler’s constant: Euler’s work and modern developments by J. C. Lagarias.
Oh well. At least I practiced solving delay differential equations in Python. There is no built-in method in SciPy for that, and although there are some modules for DDE out there, I decided to roll my own. The logic is straightforward: solve the ODE on an interval of length 1, then build an interpolating spline out of the numeric solution and use it as the right hand side in the ODE, repeat. I used Romberg’s method for integrating the solution; the integration is done separately on each interval [k, k+1] because of the lack of smoothness at the integers.
import numpy as np from scipy.integrate import odeint, romb from scipy.interpolate import interp1d numpoints = 2**12 + 1 solution = [lambda x: 1] integrals = [1] for k in range(1, 15): y0 = solution[k-1](k) t = np.linspace(k, k+1, numpoints) rhs = lambda y, x: -solution[k-1](np.clip(x-1, k-1, k))/x y = odeint(rhs, y0, t, atol=1e-15, rtol=1e-13).squeeze() solution.append(interp1d(t, y, kind='cubic', assume_sorted=True)) integrals.append(romb(y, dx=1/(numpoints-1))) total_integral = sum(integrals) print("{:.15f}".format(1/total_integral))
As a byproduct, the program found the probabilities of the random sum being in each integer interval:
- 56.15% in [0,1]
- 34.46% in [1,2]
- 8.19% in [2,3]
- 1.1% in [3,4]
- 0.1% in [4,5]
- less than 0.01% chance of being greater than 5