Meta note: I record here some thought processes for this code more thoroughly than usual because it is tricky logic and also because I have come to wonder how such programming ideas arise. Sometime text such as the following is written as a explanation of code that has already been debugged and perhaps put into production. Here it is a record of thought processes that lead to the code; it precedes and enables the code.

In the following the terms ‘cube’ and ‘plane’ are taken from the case when n=3. ‘Hyper-cube’ for [0, 1]n and ‘hyper plane’ for an n−1 dimensional linear space would be more accurate.

The limitation that all components of the argument to c be positive is awkward. I think it is tricky but efficient to eliminate this restriction. This program eliminates this restriction.

H = {x | a∙x < 1} and C = [0, 1]n = {(x1, … xn) | ∀j(0 < j ≤ n → 0 ≤ xj ≤ 1)} and we seek the content and center of gravity of H∩C. Linear function f(x) = a∙x is maximized on some vertex P of the cube. It is clear that a coordinate of P is 1 iff the corresponding component of a is positive, and otherwise 0. Here and in the code the sum of the positive components is p and of the negative components q. If p < 1 at P then the f(x) < 1 at all vertices and H∩C = C and the calculation is trivial. f(x) is minimized at the vertex Q which is opposite P. f(P) = p and f(Q) = q and thus f(c) = (p+q)/2 where c is the cube center and thus midway between P and Q. Note q ≤ 0 ≤ p and if p+q > 2 then the clip plane is closer to P and otherwise to Q. To minimize the number of vertices we visit we change coordinates and adopt as origin the closer of P and Q to the clip plane. The task is to find a new function g whose range is the same space which is 0 at the new origin and is some nonzero multiple of f, and such that f(x) = g(x) whenever f(x) = 1. For each i, x'i = x'i if a[i]>0 and x'i = 1 - x'i otherwise.

I think I must have a name ej for the points where the clip plane intercepts the coordinate axis j. ej = 1/aj.

If a[1] < 0 and a[j] > 0 for j > 1 then we adopt the vertex at (1, 0, …, 0) as the new origin and use primes on the variables for the new coordinate system. The new f' will be 0 at the new origin and 1 on the clip plane which is unmoved in this alias story. x'1 = 1 − x1 and x'j = xj for j > 1. e'1 = 1 − e1 and e'j = ej − a1/aj for j > 1. Finally a'j = 1/e'j. Perhaps we need not name e or e' in the code.
a'1 = 1/e'1 = 1/(1 - 1/a1) = a1/(a1 − 1)
a'j = 1/(1/aj − a1/aj) = aj/(1 − a1)

There remains the following problems:

We must now decide on notation. Is the argument to f a set of coordinates, or a point in the space?