Here are some preliminary functions. If you just want to read the document, then you can ignore this first cell.

:

#preliminary functions

def hyperplane_arrangement_to_fan(hyperplane_arrangement):
'''
Input: Linear hyperplane arrangement
Output: Maximal cones of the induced fan. Each maximal cone is the closure of a connected component
of R^d after removal of the hyperplane arrangement
'''
n = len(hyperplane_arrangement)
p = list(powerset(range(n)))
d = hyperplane_arrangement.ambient_dim()
cones = []
for entry in p:
signs = [1 if i in entry else -1 for i in range(n)]
inequalities = []
for i in range(n):
hyperplane = hyperplane_arrangement[i]
A = hyperplane.Hrepresentation().A()
ieq = +[signs[i]*a for a in A]
inequalities.append(ieq)
C = Polyhedron(ieqs = inequalities)
if C.dimension() == d:
cones.append(C)
return cones

def intersecting_edges(polytope, u):
'''
Input: A polytope and a normal vector u
Output: Set of edges of polytope that has nonempty intersection with the hyperplane orthogonal to u
'''
u_perp = Polyhedron(eqns= [ +list(u) ])
intersecting_edges = []
for edge in polytope.faces(1):
if edge.as_polyhedron().intersection(u_perp).dimension()>-1:
intersecting_edges.append(edge)
return intersecting_edges


# Computing the intersection body of the cube.¶

We consider the $$3$$-dimensional cube $$[-1,1]^3$$ and compute the radial function of its intersection body step by step.

:

P = polytopes.cube()


The combinatorial type of the intersection of $$P$$ with a hyperplane through the origin is fixed for each region of the normal fan of the following zonotope Z. You can visualize this fan by looking at the dual polytope.

:

Z = sum([Polyhedron(vertices = [-vector(v),vector(v)]) for v in P.vertices()])
dualZ = Z.polar()
dualZ.show(frame=False, point = False)


We now construct the associated hyperplane arrangement $$H$$.

:

H = list(set([Polyhedron(eqns = [+list(vertex)]) for vertex in P.vertices()]))
Fan = hyperplane_arrangement_to_fan(H)

Let us focus now on one of the open chambers of this hyperplane arrangement, and compute the radial function in here.
We sample a vector in the interior of this region.
:

C = Fan
sample_vector = sum([vector(ray) for ray in C.rays()])
sample_vector

:

(-2, 2, 2)


We now collect all facets of $$P$$ that intersect the hyperplane perp orthogonal to sample_vector. Note that in this case $$Q = P\cap$$perp is a hexagon: indeed we are intersecting $$6$$ facets of $$P$$.

:

perp = Polyhedron(eqns= [ +list(sample_vector) ])

intersecting_facets = []
if P.intersection(perp).dimension()==2:
for facet in P.facets():
if facet.as_polyhedron().intersection(perp).dimension() == 1:
intersecting_facets.append(facet)
intersecting_facets

:

[A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices,
A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices]


The vertices of such a hexagon arise as the intersection of perp with some edges of the above facets of $$P$$. We now compute these edges.

:

intersecting_edges_of_facets = []
for facet in intersecting_facets:
edges = intersecting_edges(facet.as_polyhedron(),sample_vector)
intersecting_edges_of_facets.append(edges)
intersecting_edges_of_facets

:

[[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices],
[A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices,
A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices]]


Let $$x \in C$$ arbitrary. By the theory we know that $$x^\perp$$ intersects $$P$$ in the set of edges in intersecting_edges_of_facets. Therefore, for each such edge $$\hbox{conv}(a,b)$$ we write the corresponding vertex $$v$$ of $$P\cap x^\perp$$ as:

$v = \frac{(b\cdot x)\,a - (a\cdot x)\,b}{(b-a)\cdot x} \,.$
:

x = vector(var('x1,x2,x3'))
x1, x2, x3 = x
v = []
for edges_of_facet in intersecting_edges_of_facets:
intersecting_points_of_facets = []
for edge in edges_of_facet:
[a,b] = list(edge.vertices())
a = vector(a)
b = vector(b)
intersection_point = ((b.inner_product(x)*a - a.inner_product(x)*b )/ (b-a).inner_product(x))
intersecting_points_of_facets.append(intersection_point)
v.append(intersecting_points_of_facets)
v

:

[[((x2 - x3)/x1, -1, 1), (-1, -1, (x1 + x2)/x3)],
[(-(x2 - x3)/x1, 1, -1), (-1, (x1 + x3)/x2, -1)],
[(-1, (x1 + x3)/x2, -1), (-1, -1, (x1 + x2)/x3)],
[(1, -(x1 + x3)/x2, 1), ((x2 - x3)/x1, -1, 1)],
[(1, 1, -(x1 + x2)/x3), (-(x2 - x3)/x1, 1, -1)],
[(1, -(x1 + x3)/x2, 1), (1, 1, -(x1 + x2)/x3)]]


We are now ready to compute the volume of the hexagon. This will be done by subdividing the polygon into simplices, and by computing the volume of each simplex as a determinant.

:

determinants = []
for simplex in v:
M = matrix(
simplex + [x]
)
det_M = M.determinant()
if det:
sgn = sign(det_M(x1 = sample_vector, x2 = sample_vector, x3 = sample_vector))
determinants.append([sgn,M.determinant()])
determinants

:

[[-1, -x1*((x1 + x2)/x3 - 1) - x2 - (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 - x3],
[1, x1*((x1 + x3)/x2 - 1) - (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 + x2 + x3],
[1,
x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3],
[-1, -x1*((x1 + x3)/x2 - 1) + (x2 + (x1 + x3)*x3/x2)*(x2 - x3)/x1 - x2 - x3],
[1, x1*((x1 + x2)/x3 - 1) + x2 + (x2 - x3)*((x1 + x2)*x2/x3 + x3)/x1 + x3],
[1,
x1*((x1 + x2)*(x1 + x3)/(x2*x3) - 1) + x2 + (x1 + x2)*x2/x3 + (x1 + x3)*x3/x2 + x3]]


Finally, by adding all the determinants with a correction factor, we obtain the radial function of $$IP$$ restricted to the region $$C$$.

:

correcting_factor = 1/(factorial(P.dimension()-1) * sum([x_i^2 for x_i in x]))
rho = correcting_factor*sum([sng*det for [sng,det] in determinants])
rho.factor()

:

(x1^2 + 2*x1*x2 + x2^2 + 2*x1*x3 - 2*x2*x3 + x3^2)/(x1*x2*x3)


From the radial function we can compute the algebraic boundary of $$IP$$. In $$C$$ this is the zero set of the following polynomial.

:

p, q = rho.numerator_denominator()
alg_boundary = q-p
alg_boundary

:

x1*x2*x3 - x1^2 - 2*x1*(x2 + x3) - x2^2 + 2*x2*x3 - x3^2