The other source of ideas is from a paper referred to me by Don Chakerian. I did not understand their development but I evidently absorbed enough ideas from the paper sufficient for this code, for in retrospect I see similarities.
We have a cube [0, 1]n. The argument to c is a vector ‘a’. {x | a∙x = 1} is a plane about which we take moments. All the components of a must be positive. The task is to compute the volume of that part of the cube on the origin side of that plane, and also the moment of that part about that plane. The structure S holds such value pairs for a series of similar tetrahedra sums and differences of which amount to that part of the hypercube.
Id ‘m’ is the moment of BT.
In r(b, j):
r returns two values which multiplied respectively by the content and moment of BT provides the content and moment of other tents.
b is the value of 1 − a∙x at the cube vertex under consideration.
j is the coordinate index.
If a component of the arg of c is small cancelations errors are introduced. The good news is that omitting a zero component, thus shortening the vector, leads to the limiting result otherwise difficult to compute. E.g. c((2, 5, 0, 3)) = c((2, 5, 3)). More precisely c is symmetric in all its arguments and limit of c((a, b, c, … w, x)) as x goes to 0 is c((a, b, c, … w))