I've had to do this operation once when I was doing research on Voronoi diagrams. Heike's method is nice. Here is one possible alternative, making use of a fold-over technique to account for periodicity:

```
cPerDist = Compile[{{v1, _Real, 1}, {v2, _Real, 1}, {size, _Real}},
Norm[size/2 - Abs[size/2 - Abs[v1 - v2]]]];
PeriodicDistance[x_?MatrixQ, size_:1] := Outer[cPerDist[#1, #2, size] &, x, x, 1]
```

(You might want to look at a plot of the function $\frac12-\left|\frac12-\left|x\right|\right|$ to get an idea of how the fold-over works.)

Here, we generate and use a compiled function, since this will seem to be used many times on inexact numbers. If you need exact expressions for the lengths, it should be straightforward to write an uncompiled version of `cPerDist`

.

Try it out:

```
PeriodicDistance[{{0.1, 0, 0}, {0.9, 0, 0}}]
{{0., 0.2}, {0.2, 0.}}
pts = BlockRandom[SeedRandom[42, Method -> "Legacy"]; RandomReal[1, {12, 3}]];
PeriodicDistance[pts]
{{0.,0.369061,0.557349,0.273739,0.580937,0.607056,0.630916,0.414241,0.149751,0.214448,0.577028,0.467806},
{0.369061,0.,0.431869,0.636656,0.611335,0.628004,0.641477,0.584125,0.465563,0.34617,0.454007,0.574909},
{0.557349,0.431869,0.,0.653288,0.526924,0.323731,0.230079,0.450194,0.538027,0.589348,0.581119,0.669271},
{0.273739,0.636656,0.653288,0.,0.521814,0.433284,0.539068,0.278055,0.201587,0.392315,0.607604,0.52001},
{0.580937,0.611335,0.526924,0.521814,0.,0.333371,0.527569,0.529642,0.484959,0.691898,0.237093,0.474143},
{0.607056,0.628004,0.323731,0.433284,0.333371,0.,0.308742,0.343472,0.590752,0.69726,0.468426,0.596956},
{0.630916,0.641477,0.230079,0.539068,0.527569,0.308742,0.,0.486608,0.564332,0.543762,0.589842,0.62605},
{0.414241,0.584125,0.450194,0.278055,0.529642,0.343472,0.486608,0.,0.396313,0.377811,0.538572,0.380586},
{0.149751,0.465563,0.538027,0.201587,0.484959,0.590752,0.564332,0.396313,0.,0.307534,0.602186,0.56417},
{0.214448,0.34617,0.589348,0.392315,0.691898,0.69726,0.543762,0.377811,0.307534,0.,0.600932,0.352568},
{0.577028,0.454007,0.581119,0.607604,0.237093,0.468426,0.589842,0.538572,0.602186,0.600932,0.,0.276919},
{0.467806,0.574909,0.669271,0.52001,0.474143,0.596956,0.62605,0.380586,0.56417,0.352568,0.276919,0.}}
```

Some limited tests I did indicate that the method is at least as fast as Heike's.

How can it be less than O(n^2) in memory or speed? That's the information-theoretic lower bound. – Daniel Lichtblau – 2012-06-04T20:48:28.800

Well I am not asking for the impossible ;-). Are you implying the code I wrote is already at the theoretical lower bound? – chris – 2012-06-04T20:53:02.887

2"Millions of points" would produce a matrix with

tens of trillionsof (nonzero) entries. I don't see scaling up to that inMathematicaany time soon :-). – whuber – 2012-06-05T15:22:03.127@chris Hi Chris, great question and very instructive answers. Quick question about a slight extension, how could we e.g. extend Heike's method to non-cubic boxes? (say a cubic box where for instance one dimension is elongated.) – None – 2019-05-02T08:31:02.503

@user929304 it seems you just need to replace Norm by something anisotropic ? As in

`Norm1 := Sqrt[3] Norm[#/{1, 2}] &`

?? I have not checked. – chris – 2019-05-02T09:01:54.123