This paper seems to have achieved what I attempted here.

The idea here is to do octrees and call the potential function only at ‘corners’. What’s a corner? There is a 3D grid of points with integer coordinates. Each point is the corner of 8 unit cubes. You evaluate the potential function there just if the ‘alternating sum’ is not zero. Each adjacent cube contributes to the sum in an ‘alternating fashion’. (This works in n dimensions.)

The ‘corners’ determine the interior set. Define C thus: C((x, y, z)) = {(X, Y, Z) | X>x ∧ Y>y ∧ Z>z}. Then if S (a set of triples of integers) is the set of corners of body T then z∈T iff z belongs to an odd number of C(q) for q in S.

I describe the above again in other words. We want to describe a body bounded by unit squares whose corners have integral coordinates. We argue that it suffices to specify the set of corners. A corner is a point with integral coordinates where an odd number of octants of the body meet the point. (i.e. (x, y, z) is a corner if x, y and z are integers and of the eight points (x±½, y±½, z±½) an odd number of them are in the body.) The rule for whether a point (without integral coordinates) is in the body is whether it dominates an odd number of corners. (X, Y, Z) dominates (x, y, z) just if X>x ∧ Y>y ∧ Z>z. This seems to be for n dimensions with suitable replacement of the squares.

This has a simple constructive proof: Consider a 3D array C of bits, 0 except at the corners where it is 1. We will simply compute whether (3D bit array B) the unit cube whose center is at (i+½, j+½, k+½) belongs to the body.

FOR j TO m DO FOR k TO m DO B[j, k, 0] := 0; B[j, 0, k] := 0; B[0, j, k] := 0 OD OD; FOR j TO m DO FOR k TO m DO FOR l TO m DO B[j, k, l] := B[j-1, k, l] + B[j, k-1, l] + B[j, k, l-1] + C[j, k, l] MOD 2 OD OD OD;

Rely on this using the function u there for what we call the ‘potential’ here. The pragmatics of this is to compute the gravitational potential of an odd shaped body where distance 1 is the resolution unit of the body.

Beware the 2D body {(x, y) | 0 < x < 1 ∧ 0 < y < 1 ∨ 1 < x < 2 ∧ 1 < y < 2)}. To get the right answer we must assign a weight of 2 to the point (1, 1). We must also weight the corners with ±1.

Enumerating the corners does not well serve the necessity of visiting the entire surface of the blob. The octree serves that purpose and also finding corners.

It is slightly tricky to compute the surface of the blob as a set of squares in octrees. Closely related to this is computing the corners. Perhaps there are algorithms on for this the net.

When we move a set of cubies we must (probably) recompute potentials at every surface point.
We can compute the difference of potentials in order to avoid m^{4} work.
(m is roughly the number of cubies in one dimension across the blob — perhaps about 1000.)