How Far Apart Are Two Random Points in a Hypercube? 3

Posted by Cleve Moler,

Two days ago I wrote about random points in a square. At the last minute I added the paragraph asking about the generalization to random points in a cube. I have to admit that I didn't check the Web to see what was known about the question.

Shortly after my post appeared, Pieterjan Robbe of KU Leuven in Belgium submitted a comment pointing out that the problem has been studied extensively. He provided several pointers. The solution even has a name; it's called Robbins Constant. I have listed a few of the references below.


Robbins Constant

Robbins Constant is the value of this six-fold integral.

$$\delta = \int_0^1 \int_0^1 \int_0^1 \int_0^1 \int_0^1 \int_0^1 \! \sqrt{(x_1-y_1)^2 + (x_2-y_2)^2 + (x_3-y_3)^2} \ \mathrm{d}x_1 \mathrm{d}y_1 \mathrm{d}x_2 \mathrm{d}y_2 \mathrm{d}x_3 \mathrm{d}y_3 $$

A change of variables makes it a triple integral.

$$\delta = 8 \int_0^1\int_0^1 \int_0^1 \! \sqrt{x^2 + y^2 + z^2} \ (1-x)(1-y)(1-z) \ \mathrm{d}x \mathrm{d}y \mathrm{d}z$$

Numerical Integration

Numerical integration with the default tolerances gets about nine decimal places.

   format long
   F = @(x,y,z) sqrt(x.^2+y.^2+z.^2).*(1-x).*(1-y).*(1-z)
   delta = 8*integral3(F,0,1,0,1,0,1)
F =
  function_handle with value:
delta =

Analytic Solution

The exact analytic solution is impressive. I can't explain where it comes from.

$$\delta = \frac{1}{105}(4 + 17\sqrt{2} - 6\sqrt{3} + 21\log{(1+\sqrt{2})} + 84\log{(1+\sqrt{3})} - 42\log{(2)} - 7\pi)$$

Here is the numeric value.

   delta = (1/105)*(4 + 17*sqrt(2) - 6*sqrt(3) + 21*log(1+sqrt(2)) + ...
            84*log(1+sqrt(3)) - 42*log(2) - 7*pi)
delta =

Hypercubes and Vectorize

I also received email from MathWorks colleague Matt Tearle informing me about vecnorm, a new function that is in MATLAB Release R2017b. Unlike the norm function that computes matrix norms, vecnorm treats the elements along a specified dimension as vectors and calculates the norm of each vector. Here is the documentation page.

I don't yet have R2017b on this computer, but Matt's email prompted me to think about vectorizing the simulation, even without vecnorm. While we're at it, let's generalize to d dimensional hypercubes. And, I'll make it a one-liner.

   n = 10000000;  % n samples.
   for d = 1:10   % d dimensions.
      delta = mean(sqrt(sum((rand(n,d)-rand(n,d)).^2,2)));
      fprintf('%6d %7.4f\n',d,delta)
     1  0.3335
     2  0.5215
     3  0.6617
     4  0.7778
     5  0.8786
     6  0.9692
     7  1.0517
     8  1.1282
     9  1.1997
    10  1.2675


Thanks to Pieterjan Robbe and Matt Tearle.


D. P. Robbins, Problem 2629, Average distance between two points in a box, Amer. Mathematical Monthly, 85.4, 1978, p. 278.

D. H. Bailey, J. M. Borwein and R. E. Crandall, Advances in the Theory of Box Integrals, <>.

Johan Philip, The Probability Distribution of the Distance Between Two Random Points in a Box,

Eric W. Weisstein, Hypercube Line Picking. MathWorld, <>.

N. J. A. Sloane, editor, The On-Line Encyclopedia of Integer Sequences, Decimal Expansion of Robbins Constant, <>.

Get the MATLAB code

Published with MATLAB® R2017a


Comments are closed.

3 CommentsOldest to Newest

Paul replied on : 1 of 3
The exact form for Robbins Constant shown here: is: 1/105*(4 + 17*sqrt(2) - 6*sqrt(3) + 21*log(1+sqrt(2)) + 42*log(2+sqrt(3)) -7*pi) This form is equivalent to what you show, but I'm kind of curious how both came to be. Did someone look at the form you have and realize that it could be reduced by a single term?
denis replied on : 3 of 3
Different metrics -- L1, L2, Lmax -- are surprisingly close after scaling:
    All-pair distances D = L1/dim       ~ .33 +- .25/sqrt dim
    All-pair distances D = L2/sqrt dim  ~ .4  +- .25/sqrt dim
    All-pair distances D = Lmax/2       ~ .3 .. .4 +- .2/sqrt dim

in 3d 8d 16d 32d. Code (python, sorry about that) and log might be under