General Approach--can 3D methods be generalized?

[student5]

Hmmm... Even though a lot of them turned out to be EKF, I'd say that bilbiro-ing and thawro-ing are still techniques that could produce further families of CRFpolychora. Like creating two "eggshells" that a brute force program could further try to connect. But then again, this will go into ex-family territory and thus get a combinatorial explosion.
Anyways, Quickfur, could you maybe share your program? I'd be interested to take a look at it and maybe make it more user-friendly (I'm into GUI's these days )

[quickfur]

Hmmm... Even though a lot of them turned out to be EKF, I'd say that bilbiro-ing and thawro-ing are still techniques that could produce further families of CRFpolychora. Like creating two "eggshells" that a brute force program could further try to connect. But then again, this will go into ex-family territory and thus get a combinatorial explosion.
Anyways, Quickfur, could you maybe share your program? I'd be interested to take a look at it and maybe make it more user-friendly (I'm into GUI's these days )

Which program are you talking about? Polyview? That's the only program I have right now, and the only thing it does is to make (convex) polytopes from input coordinates and render them using POVRay. I don't have any working prototype of a program that can do brute-force searching of CRFs, or even do EKF operations. In the past I had to manually construct the CD diagrams and lace cities and then manually construct explicit coordinates for them, before Polyview could render them.

I believe Keiji actually has a working instance of Polyview (an older version) running somewhere on this site, or at least on a personal server. I have since rewritten Polyview in another language and added some features (and perhaps introduced some new bugs ), but the core functionality remains the same. You need to have explicit coordinates to construct a polytope before it can render it. Though I did add a lot of syntactic sugar for the input: you can use permutation operators like ±, apacs, epacs, etc., and you can even delete vertices from the current set of coordinates (to make it easier to construct e.g. diminishings of the 600-cell). But there's currently no way of generating a polytope from a CD diagram or lace tower / lace city directly, you have to expand it by hand first.

[student5]

D4.3.1 and D4.3.2 are really cool: I couldn't find anything EKF-like from them.

I see it as an EKF (i.e., modified Stott expansion) of some faceted prism of an icosahedral family polyhedron.

In this lace tower configuration here, I found it difficult to "see" the bilibiro's, but they are on the sides:

Code: Select all

x2x    x2o       x2o
o2F    o2f       o2f
f2o <- f2(-x) <- f2o
o2F    o2f       o2f
x2x    x2o       x2o


Ah yes, now I see it: The bilbs are expanded along the square, so the original shape is a dodecahedral-prism kind of thing:

Code: Select all

x5o3x    x5o3o                   x5o3o
o5o3F    o5o3f                   o5o3f
f5o3o <- f5o3(-x) <- f5(-x)3x <- o5x3o
o5o3F    o5o3f                   o5o3f
x5o3x    x5o3o                   x5o3o


Notice that we have to swap twice to get back to a non-faceted something. However, this is not anything as far as I can see... The ikes don't show up, because the faceting doesn't work out... I can't find it as an ekf now...

Made me think that the stott-expanded version may actually be an EKF that then gets stott-contracted:

Code: Select all

x5x3x    x5x3o       x5x3o
o5x3F    o5x3f       o5x3f
f5x3o <- f5x3(-x) <- f5o3x
o5x3F    o5x3f       o5x3f
x5x3x    x5x3o       x5x3o


This does look like it may be CRF.. Indeed it's part of rox:

Code: Select all

x3o5o
o3x5o
F3o5o
o3f5o
f3o5x
o3x5x <-
f3x5o <-
V3o5o (V=2f)
x3o5f <-
f3x5o <-
o3x5x <-
...


[student5]

And D.4.16: Could it be derived from an EKF?
from the wiki: xfox-2-oxfo-3-ooox&#xt

Code: Select all

x2o3o   x2o3x
f2x3o    f2x3x
o2f3o <- o2f3x
x2o3x    x2o3u


nope...

The mibdies are in this case on directly visible in the lace tower:

Code: Select all

mibdi:
x2o
f2x
o2f
x2o
lace city:
    o o
   x   x
     f
    o o


compare to ike:

Code: Select all

x2o
o2f
f2x
o2f
x2o


trying to expand into lace city:

Code: Select all

   o3o o3o
  x3o   x3o
     f3o
   o3x o3x <- where does this x come from?


Original thread describing the process which indeed needs some brute-forcing to find the mysterious x, which quickfur did by hand. Also, it is not EKF-derived! (not proven)
I still have this hunch that it is some contraction and diminishing of trigonally-expanded ex.../rhombochoron, trying to expand:

Code: Select all

   o3x o3x
  x3x   x3x
     f3x
   x3o x3o <- what to do here?


It doesn't seem to be part of this: lace city source

Code: Select all

                x3o               
                                 
        o3x   x3f F3o   x3x       
                                 
                                 
x3o x3f   F3x       o3F
                                 
   F3o       x3F   f3x o3x
                                 
                                 
x3x   o3F f3x   o3x

        o3x


[Klitzing]

D4.7: Can be augmented undocumented but probably some EKF?

Cf. https://bendwavy.org/klitzing/incmats/id=srid=F-ike=ti.htm.

--- rk

[Klitzing]

However, just found an other old one (http://hi.gher.space/forum/viewtopic.php?p=21746#p21746) which wasn't so far on my CRF or EKF pages.
In fact the true EKF (starting from ex = 600-cell) within A2 x A2 subsymmetry here is:

Code: Select all

ABCDEFGHIJ   ABCDEFGHIJ   ABCDEFGHIJ   ABCDEFGHIJ
-------------------------------------------------
fFoxffooxo 3 foFfxofxoo 2 oxofofxFof 3 ooxofxfoFf &#zx  (ex)
fFo(-x)ffoo(-x)o 3 foFFxofxxo 2 oxofofxFof 3 ooxofxfoFf &#zx  (D1+I1 rewrite)
FAxoFFxxox 3 foFFxofxxo 2 oxofofxFof 3 ooxofxfoFf &#zx  (first node expanded -> already contained EKF)

FAxoFFxxo. 3 foFFxofxx. 2 oxofofxFo. 3 ooxofxfoF. &#zx  (J-level diminished = student5's/Quickfur's above mentioned, but so far missing find)

.AxoFFxxo. 3 .oFFxofxx. 2 .xofofxFo. 3 .oxofxfoF. &#zx  (A+J-levels diminished = todays new find!)


Btw. did you already know of the facet totals of that old one?
Those were:

• 18 ikes
• 9 mibdies (J62)
• 54 squippies (J1)
• 84 tets
• 6 thawroes (J92)
• 12 trips

And the facet totals of my new one then are:

• 18 gyepips (J11)
• 9 mibdies (J62)
• 54 squippies (J1)
• 6 teddies (J63)
• 54 tets
• 6 thawroes (J92)
• 12 trips

--- rk

[Klitzing]

Just found a further CRF:
FxFox 5 xFoFx 2 xofox 5 oxofx &#zx
Its facet totals are:

• 50 bilbiroes (J91)
• 10 dips
• 10 paps
• 10 pips
• 50 squippies (J1)

Code: Select all

                 x5x           x5x                 
                        F5o                       
                                                   
            o5F                     o5F           
     x5x                                   x5x     
                        F5x                       
                 x5F           x5F                 
                                                   
     F5o                                   F5o     
            F5x                     F5x           
                                                   
x5x                                             x5x
                                                   
            x5F                     x5F           
     o5F                                   o5F     
                                                   
                 F5x           F5x                 
                        x5F                       
     x5x                                   x5x     
            F5o                     F5o           
                                                   
                        o5F                       
                 x5x           x5x                 


The decagonal periodic -dip-{10}-dip- ring clearly is obvious here. Then there is an orthogonal, alternating -pip-{5}-pap-gyro{5}-gyropip-gyro{5}-pap-{5}-pip- ring. And, what is quite notable, the pentagons of the bilbiroes just connect two of those each. In fact the latter ones get squeezed between those 2 orthogonal rings, while the squippies just fill in the remaining gaps. Note that it does not use any tetrahedra (nor triangular prisms or any other "usual" small stuff).

Code: Select all

FxFox5xFoFx xofox5oxofx&#zx

o....5o.... o....5o....     & | 100  *   * |  1   2   2   1   0   0   0 |  2  1   3  1   2   0   0  0  0  0 |  1  1  3  0  0
..o..5..o.. ..o..5..o..     & |   * 50   * |  0   0   0   2   4   0   0 |  0  0   0  1   4   2   2  0  0  0 |  0  0  4  1  0
....o5....o ....o5....o       |   *  * 100 |  0   0   0   0   2   2   2 |  0  0   0  0   1   2   2  1  2  2 |  0  0  2  2  2
------------------------------+------------+----------------------------+-----------------------------------+---------------
..... x.... ..... .....     & |   2  0   0 | 50   *   *   *   *   *   * |  2  0   0  1   0   0   0  0  0  0 |  0  1  2  0  0
..... ..... x.... .....     & |   2  0   0 |  * 100   *   *   *   *   * |  1  1   1  0   0   0   0  0  0  0 |  1  1  1  0  0
oo...5oo... oo...5oo...&#x    |   2  0   0 |  *   * 100   *   *   *   * |  0  0   2  0   1   0   0  0  0  0 |  1  0  2  0  0
o.o..5o.o.. o.o..5o.o..&#x  & |   1  1   0 |  *   *   * 100   *   *   * |  0  0   0  1   2   0   0  0  0  0 |  0  0  3  0  0
..o.o5..o.o ..o.o5..o.o&#x  & |   0  1   1 |  *   *   *   * 200   *   * |  0  0   0  0   1   1   1  0  0  0 |  0  0  2  1  0
....x ..... ..... .....     & |   0  0   2 |  *   *   *   *   * 100   * |  0  0   0  0   0   1   0  1  1  1 |  0  0  1  1  2
..... ..... ....x .....     & |   0  0   2 |  *   *   *   *   *   * 100 |  0  0   0  0   0   0   1  0  1  1 |  0  0  1  1  1
------------------------------+------------+----------------------------+-----------------------------------+---------------
..... x.... x.... .....     & |   4  0   0 |  2   2   0   0   0   0   0 | 50  *   *  *   *   *   *  *  *  * |  0  1  1  0  0
..... ..... x....5o....     & |   5  0   0 |  0   5   0   0   0   0   0 |  * 20   *  *   *   *   *  *  *  * |  1  1  0  0  0
..... ..... xo... .....&#x  & |   3  0   0 |  0   1   2   0   0   0   0 |  *  * 100  *   *   *   *  *  *  * |  1  0  1  0  0
..... x.o.. ..... .....&#x  & |   2  1   0 |  1   0   0   2   0   0   0 |  *  *   * 50   *   *   *  *  *  * |  0  0  2  0  0
ooooo5ooooo ooooo5ooooo&#xr   |   2  2   1 |  0   0   1   2   2   0   0 |  *  *   *  * 100   *   *  *  *  * |  0  0  2  0  0
..... ..o.x ..... .....&#x  & |   0  1   2 |  0   0   0   0   2   1   0 |  *  *   *  *   * 100   *  *  *  * |  0  0  1  1  0
..... ..... ..... ..o.x&#x  & |   0  1   2 |  0   0   0   0   2   0   1 |  *  *   *  *   *   * 100  *  *  * |  0  0  1  1  0
....x5....x ..... .....       |   0  0  10 |  0   0   0   0   0  10   0 |  *  *   *  *   *   *   * 10  *  * |  0  0  0  0  2
....x ..... ....x .....     & |   0  0   4 |  0   0   0   0   0   2   2 |  *  *   *  *   *   *   *  * 50  * |  0  0  0  1  1
....x ..... ..... ....x     & |   0  0   4 |  0   0   0   0   0   2   2 |  *  *   *  *   *   *   *  *  * 50 |  0  0  1  0  1
------------------------------+------------+----------------------------+-----------------------------------+---------------
..... ..... xo...5ox...&#x    |  10  0   0 |  0  10  10   0   0   0   0 |  0  2  10  0   0   0   0  0  0  0 | 10  *  *  *  *  pap
..... x.... x....5o....     & |  10  0   0 |  5  10   0   0   0   0   0 |  5  2   0  0   0   0   0  0  0  0 |  * 10  *  *  *  pip
FxFox ..... ..... oxofx&#xt & |   6  4   4 |  2   2   4   6   8   2   2 |  1  0   2  2   4   2   2  0  0  1 |  *  * 50  *  *  bilbiro (tower: 21435)
..... ..o.x ..... ..o.x&#x  & |   0  1   4 |  0   0   0   0   4   2   2 |  0  0   0  0   0   2   2  0  1  0 |  *  *  * 50  *  squippy
....x5....x ....x .....     & |   0  0  20 |  0   0   0   0   0  20  10 |  0  0   0  0   0   0   0  2  5  5 |  *  *  *  * 10  dip


--- rk

[mr_e_man]

Nice find, Klitzing!
It is interesting, how the bilbiros fit together.

I encountered that vertex figure (2 dip + 2 squippy + 2 bilbiro/thawro/pocuro) in my search, but I didn't dwell on it, as I was focusing on 7-gons etc. rather than 10-gons.

[mr_e_man]

Here are my thoughts on the general approach for n>4 dimensions, as it doesn't appear to be explicitly spelled out anywhere.


([j, k]-subpolytopes are defined here.)

Assuming we've enumerated all CRH polytopes in the previous dimension, these are the possible [-1, n-1]-spyts (i.e. faces, using Wendy's terminology) of any CRH polytope in n dimensions.

Of course this also determines the possible [j, n-1]-spyts for all j (the vertex figures and edge figures etc. of the faces), including the [n-4, n-1]-spyts, which can be considered as spheric polygons.

Then any [n-4, n]-spyt can be considered as a spheric polyhedron, with faces from that known set of polygons. So we try to build polyhedra in 3D spheric space. By Cauchy's rigidity theorem, the polyhedron's dihedral angles are determined by its arrangement of faces. These angles are the [n-2, n]-spyts.

So we know the possible [-1, n-1] and [n-4, n]-spyts. Then we need to fit these subpolytopes together, matching the [n-4, n-1]-spyts between the two types, to build the whole polytope. This task could be broken down further: First enumerate the [n-4, n]-spyts, then enumerate the [n-5, n]-spyts using the known [n-5, n-1] and [n-4, n]-spyts, then enumerate the [n-6, n]-spyts using the known [n-6, n-1] and [n-5, n]-spyts, and so on. In other words, at each step (expanding dimensionally downward), consider faces and verfs to determine the whole subpolytope.


The point is that all the flexibility is confined to 3D. The rest of the problem is discrete/combinatorial.

[mr_e_man]

Assuming we've enumerated all CRH polytopes in the previous dimension

Of course we need to understand "enumerate" loosely, since we've seen that there are more than 10^64 different polychora derived from the o5o3o3o family. For example we could say "classify" instead.

Or we could just enumerate the vertex figures, not bothering with the whole polychora, and use this data to enumerate the edge figures of polytera.

[mr_e_man]

Hello again. If I can remember it well, one of the steps of enumerating johnson solids involved the vertex angle defect.
In 2 dimensions, the exterior angles of every non-self-intersecting polygon sum up to 2π because it seems like they can be deformed into a circle. The case in 3 dimensions sums up to 4π, equivalent to a sphere. My question is: is there an equivalent of that in 4 dimensions? The complete theorem is the Gauss-Bonnet theorem, but when I tried to find such relation, I drowned in heavy calculus. So I was curious to know if someone found this result.

If it's equivalent to hypersphere, then it would be 2π^2.

But if that's the case, there's some kind of angle defect which has this sum. What would it be? Dihedral angles and vertex solid angles are not comensurable with π, except for the cube (which can tile the space alone) and, maybe, truncated octahedron. https://en.wikipedia.org/wiki/Platonic_solid#Geometric_properties
That led me to question if this invariant even exists, which it does because of that theorem.

The duals of the verfs must exactly cover the sphere. (The dual of the whole polytope is a discrete version of the Gauss map of a smooth surface.) In particular, their (n-1)-dimensional measures must add up to exactly the measure of the sphere. That is the generalization of polyhedron vertex angle deficits adding up to 4π.

But the Gauss-Bonnet theorem generalizes in a different way, only to even dimensions. I don't know much about that.


This fact can be used to show (see here) that the number of CRH n-polytopes is finite, if the (Euclidean, regular) polygons have limited size.

[mr_e_man]

Also, I have a conjecture:

If a polytope, in n-dimensional Euclidean or spheric space with n>2, has convex faces (the [-1, n-1]-spyts) and convex vertex figures (the [0, n]-spyts), then the polytope (the [-1, n]-spyt) is convex.

Assuming this is true, we would also get the following, by induction on dimensions:

If a polytope, in n-dimensional Euclidean or spheric space with n>2, has convex [k, k+3]-subpolytopes for all k (including hedra with k = -1), then it is convex.

Thus, given an n-polytope with CRH faces, we only need to ensure that the [n-3, n]-spyts are convex (spheric) polygons, and then it should follow that the whole polytope is convex.


(But this depends on what a "polytope" even is, if it's not convex. That is a different topic.)

[student5]

I've read a bit of thingies here and started writing a program in Rust (Rust is much faster than Python @mr_e_man). So far it can:

parse ASCII coxeter-dynkin notation into a PetGraph (a graph datastructure)
Generate the first point for 3-dimensional coxgraphs based on a paper I cannot currently find back (I lost my laptop and didn't upload the pdf cuz copyright)
do the initial reflections

Currently open questions:

What is a nice sequence to do all the reflections and put them in a datastructure (actually another graph with coordinates at all vertices, but depending on whether there is an x or o on an edge, some coincide, so this involves some edge-cases, including what edges are there. Maybe a convex-hull algorithm could work, but knowing which edges lace within a polygon is useful for constructing lace towers
Further develop exxact crate so it can:

Do division, multiplication etc
do cos(pi/n) easily using DeMoivre's theorem, cuz hard-coding is boring and ugly

I have to go now, please look at the stuff and/or add your own programs/snippets to it in a different folder so we can cooperate!

[mr_e_man]

Thanks for the hint. I don't do a lot of programming. I know Python is slow, but it hasn't been a problem for me yet. Maybe I'll try Rust in the future.

Further develop exxact crate so it can:

Do division, multiplication etc
do cos(pi/n) easily using DeMoivre's theorem, cuz hard-coding is boring and ugly

See Minimal polynomial of 2cos(2π/n)

I have to go now, please look at the stuff and/or add your own programs/snippets to it in a different folder so we can cooperate!

I am focusing on arbitrary asymmetric things. CD diagrams are too specialized.

[mr_e_man]

Okay, here's the algorithm I'm using for the first part of my program.

(My data structures are based on flags. Each flag in an n-polytope has references to n other flags, and n measures (lengths or angles) associated with these connections. An edge in a multihedron is just a set of four flags, or two flags if the edge is incomplete.)

• We have a list of polygons, which are the allowed types of faces of the polyhedra we're trying to build. From this, we can make a sorted list of the angles in the polygons, etc.


• Given a multihedron (an incomplete polyhedron), find the incomplete vertex with the greatest angle sum. By convexity, any vertex angle sum must be less than 360°. The multihedron is likely to close up into a polyhedron more quickly if we work on this vertex. Of the two incomplete edges at this vertex, take the edge whose other vertex has the greater angle sum. We'll try to complete this edge, either by adding a new face there, or by connecting it to another existing incomplete edge. (It's tempting to only look at the other edge at the same vertex, but we must consider the possibility of other existing edges intruding between them. In other words, two incomplete vertices may combine.)


• Try to connect this edge to each existing edge:
• Make sure to connect them with the right orientation. I.e. don't accidentally make a Möbius band.
• Make sure the edge lengths match.
• Don't connect them if they already share a face. I.e. a single polygon shouldn't be twisted back onto itself.
• Don't connect them if a vertex will be complete but have only 2 faces around it.
• Add the angle sums at the corresponding vertices (unless they were already the same vertex). The resulting sum(s) must be less than 360°.
• Optionally, calculate some dihedral angles. If there are several ways to get an angle, the results must agree with each other.
• If it passes all these tests, accept the new multihedron, with the two edges connected. (Or perhaps I should say the two edges become one edge.) Of course there are other possible tests, but they're more complicated.
• (Don't replace the original multihedron. Continue trying its other edges, to make different new multihedra.)
• Try to add a new face at this edge. Go through the list of all orientations of all polygons, in order of the angles in the polygons:
• Make sure the edge lengths match.
• At the vertex with the greatest angle sum, add the corresponding angle in the polygon. The result must be less than 360°. If it's not, then anything further in the list (since it's sorted) will also make the sum exceed 360°, so we can immediately stop trying new faces here.
• At the other vertex of this edge, add the corresponding angle in the polygon. The result must be less than 360°. (Note, the list can't be sorted in two ways at the same time.)
• Accept the new multihedron, with the new face.
• (Don't replace the original multihedron. Continue trying other polygons.)
• Now we have a list of new multihedra, each containing the original multihedron as a subset. I guess I'll call this list an "elaboration" of the multihedron. (Sometimes I think of it as opening a box and taking several smaller boxes out of it.)


• Given a list of multihedra, replace each one with its elaboration, unless it's already a polyhedron, which cannot be elaborated. This generally makes the list longer, but it may make it shorter, if some multihedron elaborates to nothing, i.e. if we find out that that multihedron just doesn't work. For a depth-first search, go through the list elaborating in front of you. For a breadth-first search, elaborate behind you and keep moving forward, and if you reach the end of the list, start again from the beginning. Either way, it's supposed to eventually boil down to a list of polyhedra. (Of course this could be done using trees instead of lists. The root of the tree is the original multihedron, and the leaves of the tree are the possible polyhedra.)

[mr_e_man]

• Try to add a new face at this edge. Go through the list of all orientations of all polygons, in order of the angles in the polygons:
• Make sure the edge lengths match.
[...]

To make this part more efficient, I suppose we could further refine the polygon data, and make a separate list for each edge length.


There's something I forgot to mention.

• Add the angle sums at the corresponding vertices (unless they were already the same vertex). The resulting sum(s) must be less than 360°.

Each vertex keeps a record of its angle sum. So when an incomplete vertex with m incident faces is combined with an incomplete vertex with n incident faces, you only need to add 2 numbers. If you don't keep these records, you need to add m+n numbers. But it probably doesn't matter much (unless you have hundreds of faces around a single vertex!).

[mr_e_man]

We need to be able to solve a polygon; to calculate all of its angles, given some of its angles and all of its edge lengths.

Thus, suppose we have an n-gon (possibly irregular or non-Euclidean), and we know all n of the edge lengths and at least n-3 of the angles. If n=3 then the trigon is easily solved by the law of cosines (after checking the triangle inequality), so suppose n>3. Then, since we know at least one vertex angle, and we also know the lengths of the two edges at this vertex, we can use the law of cosines to find the distance between the two vertices adjacent to this one. The three vertices form a trigon which we can solve. Cutting off the trigon, we get an (n-1)-gon, with all n-1 edge lengths known, and at least n-4 angles known. (If an angle in the n-gon adjacent to the original angle is known, then we can subtract an angle in the trigon, to get an angle in the (n-1)-gon.) Solve this smaller polygon, using recursion. Then the larger polygon is also solved. (If an angle in the n-gon adjacent to the original angle is unknown, then we can find it by adding an angle in the trigon and an angle in the (n-1)-gon which was just calculated.)

Also, whenever you add or subtract two angles, check that the result is between 0° and 180°. Everything is supposed to be convex.


Now let's apply this to a multihedron's verfs. Expanding on this step:

• Optionally, calculate some dihedral angles. If there are several ways to get an angle, the results must agree with each other.

We were connecting two edges. If one of the two resulting vertices is complete, and n-valent, and we know at least n-3 of the dihedral angles there, then solve the verf to get the other dihedral angles. For each edge that was updated with a dihedral angle, go to its other vertex, and try to solve that verf also. Keep doing this until no more verfs can be solved.

Whenever you solve a verf where more than n-3 of the angles were already known, you can use just n-3 of them in the calculations, and afterward compare the results with the remaining known angles. If they're unequal then the multihedron can be discarded. (Floating-point arithmetic generally would say that they're unequal, even when they're actually equal. Of course the floating-point results may be close to equal, but how close is close enough? If you use interval arithmetic instead, the test is whether the two intervals overlap.)

This can only be applied if at least one vertex somewhere is trivalent.

[mr_e_man]

• [...] We'll try to complete this edge, either by adding a new face there, or by connecting it to another existing incomplete edge. (It's tempting to only look at the other edge at the same vertex, but we must consider the possibility of other existing edges intruding between them. In other words, two incomplete vertices may combine.)

Actually the temptation is harmless. We can just connect the two edges at the same vertex. Let's prove the following statement:

Any convex polyhedron (more generally, any polyhedron with spherical topology, i.e. with Euler characteristic 2) can be built, starting from a single face (more generally, from any connected subset of the polyhedron's faces), using two operations: adding a face and connecting it by just one edge; and connecting two opposing incomplete edges to form one complete edge.

(Note, it's false with toroidal topology; eventually you'd have to join two non-opposing edges.)

Here "opposing" means the edges already share a vertex before we connect them. They're at opposite ends of the incomplete vertex figure.
(I tried to think of a better word, but I couldn't.)


Let me clarify the above statement using the language of polyflags.
We have a polyhedron P, and a multihedron M consisting of a subset of the flags in P and a subset of the adjacencies between flags.
Each face of the multihedron is complete; i.e. every flag x in M has x⟨0⟩ and x⟨1⟩ defined; but some flags may have x⟨2⟩ undefined. (Even if two flags exist in M and are adjacent in P, they may not be adjacent in M.)
The multihedron is connected; any flag in M can be reached from any other flag in M by applying adjacencies in M.
An incomplete edge in M is a pair of flags x, x⟨0⟩, with x⟨2⟩ undefined. A complete edge is a quadruple of flags x, x⟨0⟩, x⟨2⟩, x⟨0⟩⟨2⟩.
You can add a face to M by adding its set of flags, all connected by ⟨0⟩ and ⟨1⟩, and adding a pair of ⟨2⟩ connections to make a complete edge.
You can join two incomplete edges to make a complete edge, if they are already connected by a chain of ⟨1⟩ and ⟨2⟩ in M (and are not actually two different edges in P, of course).

The claim is that it's possible to construct P from M by repeatedly applying those two operations to M.


Let me further clarify using pictures. Consider this arbitrary polyhedron (unrelated to CRFs) made of 2 hexagons, 2 pentagons, 4 tetragons, and 4 trigons.
Flag adjacencies ⟨0⟩ are represented by double lines, ⟨1⟩ by single lines, and ⟨2⟩ by dashed lines.

Multihedron 1.png (145.28 KiB) Viewed 23902 times

In the next images, the multihedron M is black and yellow, and the polyhedron P is light grey. A new face is green. Completing an edge is red. Opposition of two incomplete edges is blue.

Multihedron 2.png (110.76 KiB) Viewed 23902 times

Multihedron 3.png (127.02 KiB) Viewed 23902 times

Multihedron 4.png (102.15 KiB) Viewed 23902 times

Now for the proof.

Notice that if you join two edges and then add a face, you could have instead added the face first and then joined the edges. Thus, we can just add all the faces first, so M has all the flags of P, and is only missing some ⟨2⟩ connections. All that needs to be shown, is that M has at least one pair of opposing edges, so that we can continue building toward P. (We assume P is finite, so the process will end eventually.)

Consider the graph G formed by the edges of P which are incomplete in M. A vertex of G is a vertex of P which is incomplete in M. An edge of G is one edge of P but two edges of M. (In other words, you can cut P along G to make M.) A pair of opposing edges in M is just an edge incident with a monovalent vertex in G. So we ask, why must G have a monovalent vertex?

Here comes the spherical topology. G cannot have any cycles, because a cycle would make M disconnected, but we assumed that M is connected (and we maintained connectivity while adding faces to M). Hence, G must be a tree, or more generally a forest.

G has at least one edge (if we assume M ≠ P); and any vertex of G is incident with at least one edge of G (by definition of G). As any forest with V vertices and E edges consists of V-E many trees, we have V - E ≥ 1. Now let's count the vertex-edge pairs in G. On one hand, this number is just 2*E, since an edge is incident with exactly 2 vertices. On the other hand, this number is V1+2*V2+3*V3+..., where Vk is the number of k-valent vertices. Also note that V = V1+V2+V3+... . Putting all this together, we get

2*V1 + (2*V2 + 2*V3 + 2*V4 + ...) - 2 = 2*(V - 1) ≥ 2*E = V1 + (2*V2 + 3*V3 + 4*V4 + ...) ≥ V1 + (2*V2 + 2*V3 + 2*V4 + ...)

and we conclude that V1 ≥ 2; there are at least two monovalent vertices.

[mr_e_man]

Also, it's important to note that you can freely choose which edges to join and which faces to add. At first glance, the theorem in the previous post only says that there is some sequence of operations you can apply to M to make P. However, any sequence of operations will work; the operations preserve connectedness of M, so the theorem still applies, no matter what you do to M. (In particular, you don't need to add all the faces first, as the proof might suggest.)

And in practice, at least in the context of the CRH polytope search, we wouldn't know P in advance; we'd be building M blindly, expecting to get some polyhedron P containing M as a subset, or get some contradiction and thus prove such polyhedra don't exist.

[mr_e_man]

There is more room for improvement in the algorithm described above:

• Now we have a list of new multihedra, each containing the original multihedron as a subset. I guess I'll call this list an "elaboration" of the multihedron. (Sometimes I think of it as opening a box and taking several smaller boxes out of it.)
• Given a list of multihedra, replace each one with its elaboration, unless it's already a polyhedron, which cannot be elaborated. This generally makes the list longer, but it may make it shorter, if some multihedron elaborates to nothing, i.e. if we find out that that multihedron just doesn't work. For a depth-first search, go through the list elaborating in front of you. For a breadth-first search, elaborate behind you and keep moving forward, and if you reach the end of the list, start again from the beginning. Either way, it's supposed to eventually boil down to a list of polyhedra. (Of course this could be done using trees instead of lists. The root of the tree is the original multihedron, and the leaves of the tree are the possible polyhedra.)

Regarding memory usage, basically we're keeping records of a horizontal section of the tree (with BFS), or maybe a diagonal section (with DFS). It would be better to keep a vertical section (because a tree's width is exponentially greater than its height).

The recursion functionality built into most programming languages could do this automatically. I.e. we needn't explicitly keep any list of multihedra (a multihedron being a node in the tree), though they will still be there in the computer's memory. I mean, the whole tree isn't in memory, but the part of the tree between the root and the current node will be in memory.

https://en.wikipedia.org/wiki/Tree_traversal

[mr_e_man]

• We have a list of polygons, which are the allowed types of faces of the polyhedra we're trying to build. From this, we can make a sorted list of the angles in the polygons, etc.

Hmm, I'll have to say something about sorting intervals:
They should be sorted by their infima, i.e. lower bounds, i.e. left endpoints.

The primary reason for sorting the angles is that, if x ≤ y, and s is some partial sum, and s + x ≥ 360°, then it's guaranteed that s + y ≥ 360° also.

Now suppose X=[a,b] and Y=[c,d] are intervals, sorted so that a ≤ c. (Whether b ≤ c doesn't matter.)
Let Z be an interval (possibly wider than necessary) containing s+x for all x in X, so in particular Z contains s+a.
Suppose Z ≥ 360°, meaning that z ≥ 360° for all z in Z. Then in particular s+a ≥ 360°.
Take any y in Y, meaning that c ≤ y ≤ d.
Then we get s+y ≥ s+c ≥ s+a ≥ 360°.
In summary, if a calculation in interval arithmetic says s+x ≥ 360°, then it's guaranteed that s+y ≥ 360° also, without even knowing whether x ≤ y.


However... I'm implementing a compound data type, with both exact values and approximate values (intervals) allowing different kinds of calculations. I guess I'll need to be careful about mixing them; an exact s+x≥360° doesn't imply an approximate s+x≥360°. But of course s+y≥360° would still be guaranteed on the condition that x≤y. So we just need to avoid the following situation:

• x is known exactly (such as x = arccos(1/√3)), and also approximately, by an interval X containing x (such as X = [54.7300°, 54.7400°]).
• inf(X) ≤ inf(Y).
• x > y for some relevant y in Y.

I say "relevant" because Y may be wider than necessary, containing some stuff we don't care about.
(For example, for large n, an n-gonal biantiprismatic wedge has a dichoral angle y near x = arccos(1/√3). Let's say Y = [54.0100°, 54.7400°]. We don't care about the subinterval [54.7357°, 54.7400°], since we know y < x.)

To simplify the third condition, we should split Y if necessary (in the example it isn't necessary), defining a relevant subinterval left of x and a different relevant subinterval right of x. So either x ≤ y for all relevant y, or x ≥ y for all relevant y. The former case is good. In the latter case, to avoid the bad situation, we must have inf(X) > inf(Y).
(Continuing the example, a bad situation would be Y = [54.7310°, 54.7400°], narrower than X = [54.7300°, 54.7400°]. Then infimum-sorting would put X before Y, and a sum might be exactly s+x = 360°, but s+y < 360°.)

EDIT: In the program provided below, the exact and approximate calculations are separated. So infimum-sorting is sufficient. Don't worry about that "bad situation".

[mr_e_man]

All right! I've been working on my program, and I finally got some useful results out of it.
(I do intend to post it here. Next week perhaps. First I want to write more comments and docstrings, rename some variables, delete obsolete things, etc.)


I considered the non-composite pentagonal verfs: 5.3.3.3.3 A; 4.3.3.3.3 A,E,G,H,I; 3.3.3.3.3 D,E,F,G,I,J,K,L,M.

5.3.3.3.3 A (from the snub dodecahedron) cannot be used as a face in any polyhedron, except the two polyhedra which were already known: the pyramid with edge length 90° based on that pentagon (which is the verf of the snub dodecahedron prism); and an augmentation of that (since the polychoron can be augmented with a pentagonal prism pyramid).

4.3.3.3.3 A (from snic) cannot be used as a face in any polyhedron, except the one which was already known (the verf of the snic prism).

4.3.3.3.3 E,G,H,I similarly can only be used in verfs of prisms.

3.3.3.3.3 D,I,K,M similarly can only be used in verfs of prisms. And 3.3.3.3.3 E,F,J otherwise can only be used in combination with n.3.3.3 for n > 12. (I'm pretty sure those polyhedra actually don't exist. We will have to consider n.3.3.3 separately, to rule them out. The program has trouble with some calculations because of the angles near 180°. At least now it gives less than 100 results for a pair of n.3.3.3 tetragons.)

3.3.3.3.3 G (from the sphenocorona) then produced a large number of polyhedra: 146, counting the trivial prism verf. But all 145 others were just combinations of polygons of four types: 3.3.3.3.3 G and 5.5.3 and 5.3.3 and 3.3.3 . (It couldn't calculate dihedral angles because of high vertex valence. I'm not sure if all of these polyhedra actually exist. But certainly nothing else exists.)

3.3.3.3.3 L (from the hebesphenomegacorona) similarly produced 19 polyhedra, 18 of which were combinations of 3.3.3.3.3 L and 5.5.3 and 5.3.3 and 3.3.3 .

[quickfur]

Wow, this sounds useful! Are these the only possibilities, or only what the program has found so far? If the former, it would go a long way towards proving the non-existence of non-trivial crown jewels besides prisms of the 3D crown jewels.

[mr_e_man]

It won't let me upload the 'multihedron.txt' file.

It won't let me upload pictures either.

Oh well. I'll just copy the file into a code block. You'll have to figure things out without the pictures.

[mr_e_man]

I still have some doubts about my program's design, e.g. whether a function should be moved to a different class, or be split into several smaller functions (particularly 'solve_verf' is too long), or whether the 2D geometry should be in a separate module from the 3D geometry.

Also, the output of the main function ('build_polyhedra') is just a list of polyhedra, but it ought to include the intermediate multihedra as well, so that the results can be more easily understood. In other words, the output ought to be a tree, not just its leaves. What I said before, about exponential memory usage, shouldn't be a problem if we cut off all the dead branches (the multihedra that don't produce any polyhedra). The total number of nodes in the tree (after pruning) shouldn't be much greater than the number of leaves.


Well, here's my code (Python as usual), implementing the general 3D algorithm. (It depends on decball.py, my interval arithmetic module. I'm being careful not to use the standard floating-point numbers in any way.)

exact_angle.py

Code: Select all

# Copyright (c) 2025 Anonymous.
#
# Permission to use, copy, modify, and/or distribute this software for
# any purpose with or without fee is hereby granted.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
# WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS.  IN NO EVENT SHALL THE
# AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
# DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
# PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
# TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
# PERFORMANCE OF THIS SOFTWARE.


class Angle:
    """An angle may be specified both approximately and exactly.

    Supported operations include addition, subtraction, and comparisons.
    (Simplicity is intended.)

    An Angle has a 'value', which may be an integer (or some other exact
    type), or an interval known to contain the true value.  The units
    can be degrees or radians or anything; it has no effect on the
    functionality of this class.

    An Angle also has 'coefs', a list of integer coefficients.  These
    are used to express the angle exactly as a linear combination of
    some other known angles.  The first "basis angle" is assumed to be
    1, for compatibility with integers, and the other "basis angles" are
    assumed to be positive, for comparisons.  They have no other effect
    on the functionality of this class.  In particular, the "basis
    angles" are not assumed to be linearly independent over the
    integers.  An angle may be 0 even if the coefficients are not all 0.

    For example, your "basis angles" could be 1 degree, and the regular
    heptagon's 128.5714, and the square pyramid's 54.7356 .  Then the
    regular octahedron's 109.4712 can be given coefs=[0,0,2], and the
    regular tetrahedron's 70.5288 can be given coefs=[180,0,-2].  Note
    that coefs=[-900,7,0] is equivalent to coefs=[0,0,0], but this class
    has no way to know that.

    >>> from decball import B
    >>> x = Angle(B("109.4~1"), [0,0,2])
    >>> y = Angle(B("70.5~1"), [180,0,-2])
    >>> x + y == 180  # This is determined by 'coefs'
    True
    >>> x < 2*y  # This is determined by 'value'
    True
    >>> z = Angle(coefs=[-900,7,0])
    >>> z == 0
    False
    >>> z != 0
    False
    >>> z > -900
    True

    Comparisons have the same meaning as in interval arithmetic: 'True'
    means "definitely true", and 'False' means "possibly false", or
    "unknown".

    An Angle can also have a 'cosine' and a 'sine'.  If two angles have
    these, then the sum of the angles will have them calculated
    automatically, using the angle addition laws.  This doesn't actually
    use any trigonometric functions.  This class doesn't enforce any
    relation between the 'value' and the 'cosine' and 'sine'.

    An Angle can also have a 'name', as another way to determine exact
    equality.  (Angles with the same name are equal.  Angles with
    different names may or may not be equal.)  This is useful if
    coefs=None; that is, if the angle is not any known combination of
    the "basis angles".
    """

    def __init__(self, value=None, coefs=None,
                 cosine=None, sine=None, name=None):
        self.value = value
        self.coefs = coefs
        self.cosine = cosine
        self.sine = sine
        self.name = name

    def __repr__(self):
        return repr(self.value)

    def _eq0(x):
        a = x.value
        if a is not None and a == 0:
            return True
        l = x.coefs
        if l is not None and all(n == 0 for n in l):
            return True
        return False

    def _ne0(x):
        a = x.value
        if a is not None and a != 0:
            return True
        ##l = x.coefs
        ##if l is not None and any(n != 0 for n in l):
        ##    return True  # Assuming linear independence
        return False

    def _lt0(x):
        a = x.value
        if a is not None and a < 0:
            return True
        l = x.coefs
        if l is not None and all(n <= 0 for n in l) and any(n < 0 for n in l):
            return True
        return False

    def _le0(x):
        a = x.value
        if a is not None and a <= 0:
            return True
        l = x.coefs
        if l is not None and all(n <= 0 for n in l):
            return True
        return False

    def _gt0(x):
        a = x.value
        if a is not None and a > 0:
            return True
        l = x.coefs
        if l is not None and all(n >= 0 for n in l) and any(n > 0 for n in l):
            return True
        return False

    def _ge0(x):
        a = x.value
        if a is not None and a >= 0:
            return True
        l = x.coefs
        if l is not None and all(n >= 0 for n in l):
            return True
        return False

    def __pos__(x):
        return x

    def __neg__(x):
        (a, l, c, s) = (x.value, x.coefs, x.cosine, x.sine)
        if a is not None:
            a = -a
        if l is not None:
            l = [-n for n in l]
        if s is not None:
            s = -s
        return Angle(a, l, c, s)

    def __mul__(x, y):
        if not isinstance(y, int):
            return NotImplemented
        (a, l) = (x.value, x.coefs)
        if a is not None:
            a = a*y
        if l is not None:
            l = [n*y for n in l]
        return Angle(a, l)

    __rmul__ = __mul__

    def __add__(x, y):
        (ax, listx, cosx, sinx) = (x.value, x.coefs, x.cosine, x.sine)
        if isinstance(y, int):
            # An integer n, converted to an Angle, should have
            # coefficients [n, 0, 0, ... , 0].
            ay = y
            if listx is None:
                listy = None
            else:
                listy = [y] + [0]*(len(listx) - 1)
            (cosy, siny) = (None, None)
            ##(q, r) = divmod(y, 90)  # Assuming degrees
            ##if r == 0:
            ##    q %= 4
            ##    if q == 0:
            ##        (cosy, siny) = (1, 0)
            ##    elif q == 1:
            ##        (cosy, siny) = (0, 1)
            ##    elif q == 2:
            ##        (cosy, siny) = (-1, 0)
            ##    elif q == 3:
            ##        (cosy, siny) = (0, -1)
        elif isinstance(y, Angle):
            (ay, listy, cosy, siny) = (y.value, y.coefs, y.cosine, y.sine)
        else:
            (ay, listy, cosy, siny) = (y, None, None, None)
        if ax is None or ay is None:
            az = None
        else:
            az = ax + ay
        if listx is None or listy is None:
            listz = None
        else:
            if len(listx) != len(listy):
                raise ValueError
            listz = [m + n for (m, n) in zip(listx, listy)]
        if cosx is None or sinx is None or cosy is None or siny is None:
            cosz = None
            sinz = None
        else:
            cosz = cosx*cosy - sinx*siny
            sinz = sinx*cosy + cosx*siny
        return Angle(az, listz, cosz, sinz)

    __radd__ = __add__

    def __sub__(x, y):  return x + -y
    def __rsub__(x, y):  return -x + y

    def __eq__(x, y):
        if isinstance(y, Angle):
            if x.name is not None and y.name is not None and x.name == y.name:
                return True
        return (x + -y)._eq0()

    def __ne__(x, y):  return (x + -y)._ne0()
    def __lt__(x, y):  return (x + -y)._lt0()
    def __le__(x, y):  return (x + -y)._le0()
    def __gt__(x, y):  return (x + -y)._gt0()
    def __ge__(x, y):  return (x + -y)._ge0()



multihedron.py

Code: Select all

# Copyright (c) 2025 Anonymous.
#
# Permission to use, copy, modify, and/or distribute this software for
# any purpose with or without fee is hereby granted.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
# WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS.  IN NO EVENT SHALL THE
# AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
# DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
# PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
# TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
# PERFORMANCE OF THIS SOFTWARE.

"""Building convex polyhedra with given types of faces."""

from copy import deepcopy

import decball
from decball import sqrt, sind, cosd, cisd, acosd
from exact_angle import Angle


class Nonconvexity(Exception):
    pass

class PrimaryNonconvexity(Nonconvexity):
    pass

class SecondaryNonconvexity(Nonconvexity):
    pass


class Flag:
    """A triple (vertex, edge, face), all three being incident.

    An adjacent flag, self.next[i], has a different i-dimensional
    surtope, and all other surtopes the same as in self.

    A dyad measure, self.dyme[i], is the length or angle separating self
    from an adjacent flag.  This includes edge lengths (between two
    vertices), digrammal angles (between two edges), and dihedral angles
    (between two faces).

    The surtopes (i.e. vertices, edges, faces) usually are not directly
    referred to.  Instead, a surtope can be considered as the set of all
    flags that involve the surtope.
    """

    def __init__(self):
        self.next = [None]*3
        self.dyme = [None]*3

    def face_circuit(self, reverse=False, include_self=True):
        """Go around the polygon using next[0] and next[1].

        The first flag is self, the second flag is self.next[0], the
        third flag is self.next[0].next[1], the fourth flag is
        self.next[0].next[1].next[0], and so on.  It stops when it
        reaches self again (in which case the last flag is self.next[1])
        or it reaches None (in which case the polygon is incomplete, and
        self.next[1] either is never reached or is None).

        The previous paragraph assumes default arguments.  If
        reverse=True, then next[1] is done before next[0].

        This is a generator function, used to make 'for' loops.
        """
        if include_self:
            yield self
        (i, j) = (0, 1)
        if reverse:
            (i, j) = (j, i)
        flag = self
        while True:
            flag = flag.next[i]
            if flag is None:  # (flag can't be self, by orientability)
                return
            yield flag
            flag = flag.next[j]
            if flag is None or flag is self:
                return
            yield flag

    def verf_circuit(self, reverse=False, include_self=True):
        """Go around the vertex using next[1] and next[2].

        The first flag is self, the second flag is self.next[1], the
        third flag is self.next[1].next[2], the fourth flag is
        self.next[1].next[2].next[1], and so on.  It stops when it
        reaches self again (in which case the last flag is self.next[2])
        or it reaches None (in which case the vertex figure is
        incomplete, and self.next[2] either is never reached or is
        None).

        The previous paragraph assumes default arguments.  If
        reverse=True, then next[2] is done before next[1].

        This is a generator function, used to make 'for' loops.
        """
        if include_self:
            yield self
        (i, j) = (1, 2)
        if reverse:
            (i, j) = (j, i)
        flag = self
        while True:
            flag = flag.next[i]
            if flag is None:  # (flag can't be self, by orientability)
                return
            yield flag
            flag = flag.next[j]
            if flag is None or flag is self:
                return
            yield flag

    def polygon_is_isomorphic(self, flag):
        """Return True if the two flags are in equivalent positions in
        congruent polygons.

        If both flags are in the same polygon, this is a way to find its
        symmetries.

        Both polygons are assumed to be complete, and to have known edge
        lengths and vertex angles.  If they're known only by intervals,
        this function will probably return False, because it uses '=='.
        """
        # Check for equal lengths and angles.
        for (flag1, flag2) in zip(self.face_circuit(), flag.face_circuit()):
            for (m1, m2) in zip(flag1.dyme[:2], flag2.dyme[:2]):
                if m1 is None or m2 is None:
                    raise ValueError("polygon dyad measures must be specified")
                if not m1 == m2:  # (Don't replace this with 'm1 != m2')
                    return False
        # Check for equal valence.  If they're unequal, then 'zip'
        # finishes with the smaller polygon, and the larger polygon
        # isn't finished.
        if flag1.next[1] is not self or flag2.next[1] is not flag:
            return False
        return True

    def solve_verf(self):
        """Calculate all the dihedral angles at this vertex.

        Return True if successful, False if the vertex figure is
        incomplete or underdetermined (with more than 3 unknown (None)
        angles).

        This modifies the dihedrals (only at this vertex, not the
        adjacent vertices), replacing all occurrences of None with
        values calculated by trigonometry, assuming convexity.

        The digrammals and dihedrals are assumed to be of Angle type.
        Furthermore, the digrammals are assumed to have cosine and sine
        specified (not None) and of Interval type.

        Warning:  If this function is interrupted by an Exception, the
        vertex may be "broken"; the flag adjacencies may be changed.  If
        you want to preserve the original structure, make a deepcopy
        first.
        """
        unknown_count = 0
        known_flag = None
        for (i, flag) in enumerate(self.verf_circuit()):
            if flag.dyme[1] is None:
                raise ValueError("digrammals must be specified")
            if flag.dyme[2] is None:
                unknown_count += 1
            elif known_flag is None and i%2 == 0:
                known_flag = flag
        if flag.next[2] is None:
            # The vertex figure is incomplete.
            return False
        valence = (i+1)//2
        unknown_count //= 2
        if valence < 3:
            raise ValueError("vertex should have at least 3 edges")
        elif valence == 3:
            flag0 = self
            flag1 = flag0.next[1]
            flag2 = flag1.next[2]
            flag3 = flag2.next[1]
            flag4 = flag3.next[2]
            flag5 = flag4.next[1]
            assert flag5.next[2] is flag0
            # Spherical triangle's edge lengths and vertex angles:
            l0 = flag0.dyme[1]
            a0 = flag2.dyme[2]
            l1 = flag2.dyme[1]
            a1 = flag4.dyme[2]
            l2 = flag4.dyme[1]
            a2 = flag0.dyme[2]
            if l0 + l1 < l2 or l1 + l2 < l0 or l2 + l0 < l1:
                raise Nonconvexity("triangle inequality is violated")
            # Solve the triangle using the law of cosines.
            (cos0, sin0) = (l0.cosine, l0.sine)
            (cos1, sin1) = (l1.cosine, l1.sine)
            (cos2, sin2) = (l2.cosine, l2.sine)
            if None in (cos0, sin0, cos1, sin1, cos2, sin2):
                raise ValueError("Angle must have a specified cosine and sine")
            cos0_ = (cos2 - cos0*cos1)/(sin0*sin1)  # (int/int = float
            cos1_ = (cos0 - cos1*cos2)/(sin1*sin2)  #  must be avoided)
            cos2_ = (cos1 - cos2*cos0)/(sin2*sin0)
            a0_ = Angle(value=acosd(cos0_), cosine=cos0_)
            a1_ = Angle(value=acosd(cos1_), cosine=cos1_)
            a2_ = Angle(value=acosd(cos2_), cosine=cos2_)
            # Compare the calculated angles to any previously known
            # angles.  Keep the previous ones if they're more precise.
            if a0 is not None:
                if a0 != a0_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a0.value in a0_.value:  # (Interval containment)
                    a0_ = a0
            if a1 is not None:
                if a1 != a1_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a1.value in a1_.value:
                    a1_ = a1
            if a2 is not None:
                if a2 != a2_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a2.value in a2_.value:
                    a2_ = a2
            flag1.dyme[2] = flag2.dyme[2] = a0_
            flag3.dyme[2] = flag4.dyme[2] = a1_
            flag5.dyme[2] = flag0.dyme[2] = a2_
            return True
        # Solve the spherical polygon by recursion, cutting off a
        # triangle to make a simpler polygon.
        # Find some known angle.
        if unknown_count > 3:
            return False
        a = known_flag.dyme[2]
        cos_a = a.cosine
        if cos_a is None:
            (cos_a, sin_a) = cisd(a.value)
            a.cosine = cos_a
            a.sine = sin_a
        flag5 = known_flag
        flag4 = flag5.next[2]
        flag3 = flag4.next[1]
        flag0 = flag3.next[2]
        flag6 = flag5.next[1]
        flag9 = flag6.next[2]
        prev_angle = flag3.dyme[2]
        prev_length = flag4.dyme[1]
        next_length = flag5.dyme[1]
        next_angle = flag6.dyme[2]
        # Use the law of cosines to find the length of the cutting line.
        (cos_prev, sin_prev) = (prev_length.cosine, prev_length.sine)
        (cos_next, sin_next) = (next_length.cosine, next_length.sine)
        if None in (cos_prev, sin_prev, cos_next, sin_next):
            raise ValueError("Angle must have a specified cosine and sine")
        cos_l = cos_prev*cos_next + sin_prev*sin_next*cos_a
        sin_l = sqrt(1 - cos_l**2)
        l = Angle(value=acosd(cos_l), cosine=cos_l, sine=sin_l)
        # The cutting line becomes an edge, with four new flags.
        flag1 = Flag()
        flag2 = Flag()
        flag7 = Flag()
        flag8 = Flag()
        flag1.next[1] = flag8;  flag8.next[1] = flag1
        flag2.next[1] = flag7;  flag7.next[1] = flag2
        flag1.dyme[1] = flag8.dyme[1] = l
        flag2.dyme[1] = flag7.dyme[1] = l
        # Cut.
        flag0.next[2] = flag1;  flag1.next[2] = flag0
        flag3.next[2] = flag2;  flag2.next[2] = flag3
        flag6.next[2] = flag7;  flag7.next[2] = flag6
        flag9.next[2] = flag8;  flag8.next[2] = flag9
        flag0.dyme[2] = flag3.dyme[2] = None
        flag6.dyme[2] = flag9.dyme[2] = None
        # Solve the triangle.
        cos2 = (cos_next - cos_prev*cos_l)/(sin_prev*sin_l)
        cos7 = (cos_prev - cos_next*cos_l)/(sin_next*sin_l)
        a2 = Angle(value=acosd(cos2), cosine=cos2)
        a7 = Angle(value=acosd(cos7), cosine=cos7)
        flag3.dyme[2] = flag2.dyme[2] = a2
        flag6.dyme[2] = flag7.dyme[2] = a7
        # Subtract the angles in the triangle from those in the original
        # polygon (if known).
        if prev_angle is not None:
            a1 = prev_angle - a2
            if a1 < 0:
                raise Nonconvexity
            flag0.dyme[2] = flag1.dyme[2] = a1
        if next_angle is not None:
            a8 = next_angle - a7
            if a8 < 0:
                raise Nonconvexity
            flag9.dyme[2] = flag8.dyme[2] = a8
        # Solve the remaining polygon.  (Certainly it can be solved.)
        success = flag1.solve_verf()
        assert success
        # Uncut.
        flag1.next[2] = flag8.next[2] = None
        flag2.next[2] = flag7.next[2] = None
        flag0.next[2] = flag3;  flag3.next[2] = flag0
        flag9.next[2] = flag6;  flag6.next[2] = flag9
        # Add the angles in the triangle to get those in the original
        # polygon (if unknown).
        if prev_angle is None:
            a1 = flag1.dyme[2]
            prev_angle = a1 + a2
            if prev_angle > 180:
                raise Nonconvexity
        flag0.dyme[2] = flag3.dyme[2] = prev_angle
        if next_angle is None:
            a8 = flag8.dyme[2]
            next_angle = a8 + a7
            if next_angle > 180:
                raise Nonconvexity
        flag9.dyme[2] = flag6.dyme[2] = next_angle
        return True

    def solve_all(self):
        """Calculate all dihedral angles, or as many as possible.

        Start by solving this vertex figure, and travel to adjacent
        vertices via edges with newly calculated dihedrals, and solve
        those vertex figures also, and so on.

        This modifies the dihedrals.  Nothing is returned.

        The digrammals and dihedrals are assumed to be of Angle type.
        Furthermore, the digrammals are assumed to have cosine and sine
        specified (not None) and of Interval type.

        Warning:  If this function is interrupted by an Exception, a
        vertex may be "broken"; the flag adjacencies may be changed.  If
        you want to preserve the original structure, make a deepcopy
        first.
        """
        success = self.solve_verf()
        if not success:
            return
        newly_solved = [self]  # (Use 'collections.deque' instead of 'list'?)
        while newly_solved:
            flag1 = newly_solved.pop(0)
            # Try all the edges around the vertex that was solved.
            for (i, flag2) in enumerate(flag1.verf_circuit()):
                # (Don't consider each edge twice.)
                if i%2 == 1:
                    continue
                a2 = flag2.dyme[2]
                assert a2 is not None  # (solve_verf was successful)
                # Go to the adjacent vertex.
                flag3 = flag2.next[0]
                a3 = flag3.dyme[2]
                if a3 is not None:
                    # The dihedral angle has been calculated before.
                    if a3 != a2:
                        raise Nonconvexity("dihedrals are inconsistent")
                    continue
                # The dihedral angle is new.  Transfer it to the
                # adjacent vertex.
                # But the interval may be uselessly wide, e.g. with a
                # radius greater than 10 degrees.
                b = a2.value
                if isinstance(b, int) or (isinstance(b, decball.Ball)
                        and b.get_radius_precision() + b.exp <= 1):
                    flag3.dyme[2] = flag3.next[2].dyme[2] = a2
                    success = flag3.solve_verf()
                    if success:
                        newly_solved.append(flag3)


class Polygon:
    """Needs no introduction.

    The (cyclic) order of the surtopes is:
    vertex0, edge0, vertex1, edge1, vertex2, edge2, ...

    The order of the flags is:
    flag0 = (vertex0, edge0), flag1 = (vertex1, edge0),
    flag2 = (vertex1, edge1), flag3 = (vertex2, edge1), ...

    Edge length 0 is between vertices 0 and 1, and vertex angle 0 is
    between edges 0 and 1, so length i is at edge i, but angle i is at
    vertex i+1.
    """

    def __init__(self, edge_lengths, vertex_angles, name=None):
        self.edge_lengths = list(edge_lengths)
        self.vertex_angles = list(vertex_angles)
        self.valence = len(edge_lengths)
        if self.valence != len(vertex_angles):
            raise ValueError("polygon must have equal numbers of "
                             + "vertices and edges")
        flags = []
        for i in range(2*self.valence):
            flags.append(Flag())
        for i in range(self.valence):
            j = 2*i
            flags[j].next[0] = flags[j+1];  flags[j+1].next[0] = flags[j]
            flags[j].next[1] = flags[j-1];  flags[j-1].next[1] = flags[j]
            flags[j].dyme[0] = flags[j+1].dyme[0] = edge_lengths[i]
            flags[j].dyme[1] = flags[j-1].dyme[1] = vertex_angles[i-1]
        self.flags = flags
        self.name = name

    def __str__(self):
        return ("Polygon(edge_lengths=" + str(self.edge_lengths)
                + ", vertex_angles=" + str(self.vertex_angles)
                + ")")


class Multihedron:
    """An incomplete polyhedron, made of complete polygons.

    In a complete polyhedron, every flag has flag.next[2] defined (not
    None).  In a multihedron, some flags may have flag.next[2] undefined
    (None).

    A Multihedron is initialized with a single face.  The given Polygon
    is copied (using deepcopy), so it can be re-used for other faces.

    When a Multihedron is printed, the format for each line is:
    flag_index :  adjacent_flag_indices ;  dyad_measures
    """

    def __init__(self, polygon, name=None):
        if not isinstance(polygon, Polygon):
            raise TypeError
        face = deepcopy(polygon)
        self.faces = [face]
        # We'll need three different lists to update separately.
        self.flags = list(face.flags)  # Copy
        self.incomplete_flags = list(face.flags)  # Another copy
        for flag in self.flags:
            angle = flag.dyme[1]
            # Sines and cosines won't be needed for simple addition
            # (though they will be needed to calculate dihedrals).
            s = Angle(angle.value, angle.coefs, None, None)
            flag.vertex_digrammal_sum = s
        self.name = name

    def __str__(self):
        all_flags = self.flags
        lines = []  # (Lines of text)
        if self.name is not None:
            lines.append(str(self.name))
        lines.append("=" * 72)
        i = 0
        for face in self.faces:
            if face.name is not None:
                lines.append(str(face.name))
            last = len(face.flags)-1
            for (j, flag1) in enumerate(face.flags):
                space = "__" if j == last else "  "
                assert i == all_flags.index(flag1)
                index_next = [None if flag2 is None else all_flags.index(flag2)
                              for flag2 in flag1.next]
                lines.append(space + str(i) + ":" + space
                             + str(index_next) + ";" + space
                             + str(flag1.dyme) + space)
                i += 1
        lines.append("")
        return "\n".join(lines)

    def join_edges(self, flag1):
        """Return multihedron with flag1 edge joined to opposing edge.

        It's assumed that flag1 is incomplete (so flag1.next[2] is
        None).  Let flag2 be the other incomplete flag at the same
        vertex.

        self is copied, not modified.  The new multihedron has (the
        copies of) flag1 and flag2 adjacent, sharing a vertex and an
        edge.

        The "flag" argument, if it's actually an integer, is interpreted
        as the index of a flag in the multihedron.
        """
        if isinstance(flag1, int):
            flag1 = self.flags[flag1]
        if flag1.next[2] is not None:
            raise ValueError("edge is already complete")
        for (i, flag2) in enumerate(flag1.verf_circuit()):
            pass
        valence = (i+1)//2
        if valence < 3:
            raise Nonconvexity("vertex would be monovalent or divalent")
        if flag2.dyme[0] != flag1.dyme[0]:
            raise Nonconvexity("edge lengths don't match")
        # The other ends of the edges would also be joined.
        flag1_0 = flag1.next[0]
        flag2_0 = flag2.next[0]
        # (flag1_0 is not flag2_0, by orientability.)
        for (i, flag3) in enumerate(flag1_0.verf_circuit()):
            pass
        adjacent_vertices_coincide = flag3 is flag2_0
        if adjacent_vertices_coincide:
            # The other ends of the edges are already a single vertex,
            # which will be complete.
            valence = (i+1)//2
            if valence < 3:
                raise Nonconvexity("vertex would be monovalent or divalent")
        else:
            new_sum = (flag1_0.vertex_digrammal_sum
                       + flag2_0.vertex_digrammal_sum)
            if new_sum >= 360:
                raise Nonconvexity("vertex angle deficit isn't positive")
            # (Don't join the two vertices if they're at different
            # locations in a single face?  That test is too complex.
            # And it may be redundant somehow.)
        # Join the edges; make a new multihedron.
        (multihedron, flag1_, flag2_) = deepcopy((self, flag1, flag2))
        flag1_0 = flag1_.next[0]
        flag2_0 = flag2_.next[0]
        flag1_.next[2] = flag2_;  flag2_.next[2] = flag1_
        flag1_0.next[2] = flag2_0;  flag2_0.next[2] = flag1_0
        for flag in (flag1_, flag2_, flag1_0, flag2_0):
            multihedron.incomplete_flags.remove(flag)
        if not adjacent_vertices_coincide:
            # Two vertices are joined.  Update the angle sums.
            for flag in flag1_0.verf_circuit():
                flag.vertex_digrammal_sum = new_sum
            for flag in flag2_0.verf_circuit():
                flag.vertex_digrammal_sum = new_sum
        # Calculate dihedral angles, if possible.
        flag1_.solve_all()
        if adjacent_vertices_coincide:
            flag1_0.solve_all()
        return multihedron

    def add_face(self, flag1, polygon, flag2=0):
        """Return multihedron with new face at flag1.

        flag1 is in self, and flag2 is in polygon.  All four arguments
        are copied, not modified.  The new multihedron has (the copies
        of) flag1 and flag2 adjacent, sharing a vertex and an edge.

        The polygon is assumed to be disconnected from other polygons
        (so each flag.next[2] is None); it shouldn't be a face of a
        previously existing multihedron.

        A "flag" argument, if it's actually an integer, is interpreted
        as the index of a flag in the multihedron or the polygon.

        PrimaryNonconvexity is raised, if the sum of digrammal angles
        around the vertex is at least 360, according to interval
        arithmetic.  SecondaryNonconvexity is for the adjacent vertex at
        the same edge.  If one of these sums is at least 360 according
        to exact arithmetic, only Nonconvexity is raised.
        """
        if not isinstance(polygon, Polygon):
            raise TypeError
        if isinstance(flag1, int):
            flag1 = self.flags[flag1]
        if isinstance(flag2, int):
            flag2 = polygon.flags[flag2]
        if flag1.next[2] is not None:
            raise ValueError("edge is already complete")
        # Check the angle sum at this vertex.
        new_sum = flag1.vertex_digrammal_sum + flag2.dyme[1]
        if new_sum.value >= 360:  # (Interval comparison)
            raise PrimaryNonconvexity("vertex angle deficit isn't positive")
        elif new_sum == 360:
            raise Nonconvexity("vertex angle deficit isn't positive")
        # Check the angle sum at the adjacent vertex.
        flag1_0 = flag1.next[0]
        flag2_0 = flag2.next[0]
        new_sum_0 = flag1_0.vertex_digrammal_sum + flag2_0.dyme[1]
        if new_sum_0.value >= 360:
            raise SecondaryNonconvexity("vertex angle deficit isn't positive")
        elif new_sum_0 == 360:
            raise Nonconvexity("vertex angle deficit isn't positive")
        # Add the face; make a new multihedron.
        (multihedron, flag1_) = deepcopy((self, flag1))
        (face, flag2_) = deepcopy((polygon, flag2))
        multihedron.faces.append(face)
        multihedron.flags.extend(face.flags)
        multihedron.incomplete_flags.extend(face.flags)
        for flag in face.flags:
            angle = flag.dyme[1]
            s = Angle(angle.value, angle.coefs, None, None)
            flag.vertex_digrammal_sum = s
        flag1_0 = flag1_.next[0]
        flag2_0 = flag2_.next[0]
        # Connect the face.
        flag1_.next[2] = flag2_;  flag2_.next[2] = flag1_
        flag1_0.next[2] = flag2_0;  flag2_0.next[2] = flag1_0
        for flag in (flag1_, flag2_, flag1_0, flag2_0):
            multihedron.incomplete_flags.remove(flag)
        # Update the angle sums.
        for flag in flag1_.verf_circuit():
            flag.vertex_digrammal_sum = new_sum
        for flag in flag2_.verf_circuit():  # (Only two flags here)
            flag.vertex_digrammal_sum = new_sum
        for flag in flag1_0.verf_circuit():
            flag.vertex_digrammal_sum = new_sum_0
        for flag in flag2_0.verf_circuit():  # (Only two flags here)
            flag.vertex_digrammal_sum = new_sum_0
        return multihedron

    def build_polyhedra(self, oriented_polygons, verbosity=1):
        """Return list of all possible completions of self.

        The new faces being added to the multihedron are taken from
        oriented_polygons.  This argument should be a kind of 2D array:
        a dictionary indexed by edge length, where each entry is a list
        sorted by angle, where each entry is an oriented polygon.  And
        an oriented polygon is a tuple (polygon, flag_in_polygon).

        To be more precise, each list (in oriented_polygons) should be
        sorted by inf(angle.value), assuming each angle.value is an
        interval.  (The infimum of an interval [a,b] is a.)

        Each returned polyhedron may be printed as soon as it's found,
        not waiting for this function to finish.  This is determined by
        verbosity, as follows:

        verbosity=0:  Don't print anything
        verbosity=1:  Print a dot, as part of an extended ellipsis
        verbosity=2:  Print the number of faces in the polyhedron
        verbosity=3:  Print the names of all faces in the polyhedron
        verbosity=4:  Print the whole polyhedron, with flag data
        """
        if not self.incomplete_flags:
            # The multihedron is complete, i.e. a polyhedron.
            if verbosity == 1:
                print(".", end="")
            elif verbosity == 2:
                print(len(self.faces), end=" ")
            elif verbosity == 3:
                for face in self.faces:
                    print(face.name, end="  ")
                print("")
            elif verbosity == 4:
                print(self)
            return [self]
        sorted_flags = sorted(  # (Use 'heapq' module instead?)
            self.incomplete_flags,
            key=(lambda flag: (flag.vertex_digrammal_sum,
                               flag.next[0].vertex_digrammal_sum)),
            reverse=True)
        flag1 = sorted_flags[0]
        polyhedra = []
        try:
            multihedron = self.join_edges(flag1)
        except Nonconvexity:
            pass
        else:
            polyhedra += multihedron.build_polyhedra(oriented_polygons,
                                                     verbosity)
        edge_length = flag1.dyme[0]
        for (polygon, flag2) in oriented_polygons[edge_length]:
            try:
                multihedron = self.add_face(flag1, polygon, flag2)
            except PrimaryNonconvexity:
                # The polygon list is sorted; if one angle is too large,
                # then all later angles are also too large.
                break
            except Nonconvexity:
                continue
            polyhedra += multihedron.build_polyhedra(oriented_polygons,
                                                     verbosity)
        return polyhedra



It can tell that a vertex is improper (i.e. flat, in the middle of a face or edge of the convex hull, with a 360° digrammal sum), at least in some cases. But it can't tell that an edge is improper (with a 180° dihedral), and thus it can't exclude such edges. We might as well allow improper edges in general.

In the context of 4D, this means that a hedron with a 180° dichoral angle is allowed, though an edge with a 360° dihedral sum is not allowed. Thus we needn't explicitly include composite Johnson solids (such as J87) as polychoron faces; they will appear naturally anyway.

[mr_e_man]

This next module is just defining and processing the polygon data. Compare to List of verfs of CRF polyhedra.

crf_polyhedron_verfs.py

Code: Select all

from decball import Ball, B, cisd

from multihedron import deepcopy, Angle, Polygon, Multihedron


# This is a workaround for decball not supporting singleton intervals.
# p is the number of digits after the decimal point.
def integer_infimum(x, p):
    """Given interval x, return inf(x)*10**p, if that's an integer.

    >>> integer_infimum(B("12.345~6"), 3)
    12339
    >>> integer_infimum(B("12.345~6"), 6)
    12339000
    """
    if not (isinstance(p, int) and p >= 0):
        raise TypeError
    if isinstance(x, int):
        return x * 10**p
    elif not isinstance(x, Ball):
        raise TypeError
    if x.exp < -p:
        raise ValueError
    return (x.mid - x.rad) * 10**(x.exp + p)
# (Use standard 'decimal' module for singletons instead?)


# An angle, if known exactly, is given as a combination of four basic
# angles: 1 degree, and the pentagonal cupola's 31.7175, and the
# pentagonal pyramid's 37.3774, and the square pyramid's 54.7356 .

# Dihedral angles in CRF polyhedra with faces no larger than 12-gons:
sorted_angles = angles_3_12 = [
    # Commented-out angles only appear in composite verfs (where the
    # decomposition doesn't introduce new vertices).
    Angle(  "31.7175", [  0,  1,  0,  0]),
    Angle(  "37.3774", [  0,  0,  1,  0]),
    Angle(  "45.0000", [ 45,  0,  0,  0]),
    Angle(  "54.7356", [  0,  0,  0,  1]),
    Angle(  "60.0000", [ 60,  0,  0,  0]),          # pr3
    Angle(  "63.4349", [  0,  2,  0,  0]),
  ##Angle(  "69.0948", [  0,  1,  1,  0]),
    Angle(  "70.5288", [180,  0,  0, -2]),
    Angle(  "72.9730" ),
  ##Angle(  "74.7547", [  0,  0,  2,  0]),
    Angle(  "79.1877", [180, -2, -1,  0]),
    Angle(  "86.7268" ),
    Angle(  "90.0000", [ 90,  0,  0,  0]),          # pr4
                                            #(ap13)
    Angle(  "94.3592" ),                    # ap12
    Angle(  "94.7616" ),                    # ap11
  ##Angle(  "95.1524", [  0,  3,  0,  0]),
    Angle(  "95.2466" ),                    # ap10
    Angle(  "95.8430" ),                    # ap9
    Angle(  "96.1983" ),
    Angle(  "96.5945" ),                    # ap8
    Angle(  "97.4555" ),
    Angle(  "97.5723" ),                    # ap7
    Angle(  "98.8994" ),                    # ap6
  ##Angle(  "99.7356", [ 45,  0,  0,  1]),
    Angle( "100.1939" ),
    Angle( "100.8123", [  0,  2,  1,  0]),  # ap5
    Angle( "102.5238" ),
    Angle( "103.8362" ),                    # ap4
    Angle( "108.0000", [108,  0,  0,  0]),          # pr5
    Angle( "109.4712", [  0,  0,  0,  2]),  # ap3
    Angle( "109.5240" ),
    Angle( "110.9052", [180, -1, -1,  0]),
    Angle( "111.7348" ),
    Angle( "114.6452" ),
  ##Angle( "114.7356", [ 60,  0,  0,  1]),
    Angle( "116.5651", [180, -2,  0,  0]),
    Angle( "117.0190" ),
    Angle( "117.3556" ),
    Angle( "118.8922" ),
    Angle( "120.0000", [120,  0,  0,  0]),          # pr6
    Angle( "121.7175", [ 90,  1,  0,  0]),
    Angle( "121.7432" ),
    Angle( "124.7019" ),
    Angle( "125.2644", [180,  0,  0, -1]),
  ##Angle( "126.8699", [  0,  4,  0,  0]),
  ##Angle( "126.9641" ),
  ##Angle( "127.3774", [ 90,  0,  1,  0]),
    Angle( "127.5516" ),                    # ap4
    Angle( "128.4960" ),
    Angle( "128.5714" ),                            # pr7
    Angle( "129.4446" ),
    Angle( "131.4416" ),
  ##Angle( "132.6240" ),
    Angle( "133.5912" ),
    Angle( "133.9728" ),
    Angle( "135.0000", [135,  0,  0,  0]),          # pr8
    Angle( "135.9915" ),
    Angle( "136.3359" ),
    Angle( "137.2401" ),
    Angle( "138.1897", [  0,  2,  2,  0]),  # ap5
    Angle( "140.0000" ),                            # pr9
  ##Angle( "141.0576", [360,  0,  0, -4]),
    Angle( "141.3411" ),
  ##Angle( "141.5945" ),
    Angle( "142.6226", [180,  0, -1,  0]),
    Angle( "142.9834" ),
    Angle( "143.4787" ),
    Angle( "143.7383" ),
    Angle( "144.0000", [144,  0,  0,  0]),          # pr10
    Angle( "144.1436" ),
    Angle( "144.7356", [ 90,  0,  0,  1]),
    Angle( "145.2219" ),                    # ap6
    Angle( "145.4406" ),
    Angle( "147.2727" ),                            # pr11
    Angle( "148.2825", [180, -1,  0,  0]),
    Angle( "148.4340" ),
    Angle( "149.5648" ),
    Angle( "150.0000", [150,  0,  0,  0]),          # pr12
    Angle( "150.2223" ),                    # ap7
  ##Angle( "151.3301" ),
  ##Angle( "152.1911" ),
                                                    #(pr13)
    Angle( "152.9299" ),
    Angle( "152.9756" ),
    Angle( "153.2346" ),
  ##Angle( "153.4349", [ 90,  2,  0,  0]),
  ##Angle( "153.6350" ),
  ##Angle( "153.9424", [180, -2,  1,  0]),
    Angle( "153.9624" ),                    # ap8
    Angle( "154.4188" ),
    Angle( "154.7223" ),
    Angle( "156.8662" ),                    # ap9
    Angle( "157.1481" ),
  ##Angle( "158.3754", [360, -4, -2,  0]),
  ##Angle( "158.5718" ),
  ##Angle( "158.6816" ),
                                                    #(pr17)
    Angle( "159.0948", [ 90,  1,  1,  0]),
    Angle( "159.1865" ),                    # ap10
    Angle( "159.8924" ),
                                                    #(pr18)
  ##Angle( "160.5288", [270,  0,  0, -2]),
    Angle( "161.0833" ),                    # ap11
    Angle( "161.4828" ),
    Angle( "162.6628" ),                    # ap12
  ##Angle( "162.7356", [108,  0,  0,  1]),
                                            #(ap13)
    Angle( "164.1754" ),
  ##Angle( "164.2068", [  0,  0,  0,  3]),
    Angle( "164.2574" ),
  ##Angle( "164.2596" ),
    Angle( "166.4406" ),
    Angle( "166.8114" ),
  ##Angle( "169.1877", [270, -2, -1,  0]),
  ##Angle( "169.4282" ),
  ##Angle( "169.4712", [ 60,  0,  0,  2]),
  ##Angle( "170.2644", [225,  0,  0, -1]),
  ##Angle( "171.3411", [180,  2,  1, -2]),
    Angle( "171.6457" ),
  ##Angle( "171.7546" ),
  ##Angle( "174.3401", [180,  1, -1,  0]),
  ##Angle( "174.4343" ),
  ##Angle( "174.7356", [120,  0,  0,  1]),
]
# Convert str to Ball:
for x in sorted_angles:
    x.name = x.value
    x.value = B(x.name + "0~5")  # Error bounds are 0.5 ULP

# Do the approximate values agree with the exact values?
a1 = B("31.71750~5")
a2 = B("37.37740~5")
a3 = B("54.73560~5")
for x in sorted_angles:
    l = x.coefs
    if l is not None:
        assert not (x.value != l[0] + l[1]*a1 + l[2]*a2 + l[3]*a3)
# Is the list properly sorted?
for i in range(len(sorted_angles)-1):
    x = sorted_angles[i]
    y = sorted_angles[i+1]
    assert x.value < y.value

pr13_17 = Angle(B("155.6~33"))  # prisms 13-17
pr18_ = Angle(B("170~10"))  # prisms 18+
pr_ = Angle(B("90.00000~5"), [90, 0, 0, 0])
ap13_17 = Angle(B("165.8~20"))  # antiprisms 13-17
ap18_ = Angle(B("174.2~58"))  # antiprisms 18+
ap_13_17 = Angle(B("93.54~48"))  # antiprisms 13-17
ap_18_ = Angle(B("91.4~15"))  # antiprisms 18+

interval_angles = angles_13_ = [
    pr13_17, pr18_,
    pr_,
    ap13_17, ap18_,
    ap_13_17, ap_18_,
]
# Some intervals overlap, so they can't be directly sorted.
# However, they can still be infimum-sorted.

all_angles = angles_3_12 + angles_13_

# Calculate sines and cosines, and store them with the angles.
for x in all_angles:
    (cosx, sinx) = cisd(x.value)
    x.cosine = cosx
    x.sine = sinx


verfs_3_12 = [
    # Commented-out polygons are composite.
    ([12,4,4], [90_0000, 150_0000, 90_0000], "12.4.4"),
    ([11,4,4], [90_0000, 147_2727, 90_0000], "11.4.4"),
    ([ 9,4,4], [90_0000, 140_0000, 90_0000], "9.4.4"),
    ([ 7,4,4], [90_0000, 128_5714, 90_0000], "7.4.4"),
    ([12,3,3,3], [94_3592, 162_6628, 162_6628, 94_3592], "12.3.3.3"),
    ([11,3,3,3], [94_7616, 161_0833, 161_0833, 94_7616], "11.3.3.3"),
    ([ 9,3,3,3], [95_8430, 156_8662, 156_8662, 95_8430], "9.3.3.3"),
    ([ 7,3,3,3], [97_5723, 150_2223, 150_2223, 97_5723], "7.3.3.3"),
    ([10,10,3], [116_5651, 142_6226, 142_6226], "10.10.3"),
    ([10, 6,4], [142_6226, 159_0948, 148_2825], "10.6.4"),
    ([10, 5,4], [116_5651, 148_2825, 121_7175], "10.5.4"),
    ([10, 5,3], [ 63_4349, 142_6226,  79_1877], "10.5.3"),
    ([10, 4,4], [ 90_0000, 144_0000,  90_0000], "10.4.4"),
    ([10, 4,3], [ 31_7175, 159_0948,  37_3774], "10.4.3"),
    ([ 8, 8,3], [ 90_0000, 125_2644, 125_2644], "8.8.3"),
    ([ 8, 6,4], [125_2644, 144_7356, 135_0000], "8.6.4"),
    ([ 8, 4,4], [ 90_0000, 135_0000,  90_0000], "8.4.4"),
    ([ 8, 4,3], [ 45_0000, 144_7356,  54_7356], "8.4.3"),
    ([ 6, 6,5], [138_1897, 142_6226, 142_6226], "6.6.5"),
    ([ 6, 6,4], [109_4712, 125_2644, 125_2644], "6.6.4"),
    ([ 6, 6,3], [ 70_5288, 109_4712, 109_4712], "6.6.3"),
    ([ 6, 4,4], [ 90_0000, 120_0000,  90_0000], "6.4.4"),
    ([ 6, 4,3], [ 54_7356, 125_2644,  70_5288], "6.4.3"),
    ([ 5, 5,5], [116_5651, 116_5651, 116_5651], "5.5.5"),
    ([ 5, 5,3], [ 63_4349, 100_8123, 100_8123], "5.5.3"),
    ([ 5, 4,4], [ 90_0000, 108_0000,  90_0000], "5.4.4"),
    ([ 5, 3,3], [ 37_3774, 138_1897,  37_3774], "5.3.3"),
    ([ 4, 4,4], [ 90_0000,  90_0000,  90_0000], "4.4.4"),
    ([ 4, 4,3], [ 60_0000,  90_0000,  90_0000], "4.4.3"),
    ([ 4, 3,3], [ 54_7356, 109_4712,  54_7356], "4.3.3"),
    ([ 3, 3,3], [ 70_5288,  70_5288,  70_5288], "3.3.3"),
  ##([10,3,4,3], [142_6226, 174_3401, 159_0948, 153_9424], "10.3.4.3"),
    ([10,3,3,3], [ 95_2466, 159_1865, 159_1865,  95_2466], "10.3.3.3"),
  ##([ 8,3,4,3], [125_2644, 170_2644, 144_7356, 144_7356], "8.3.4.3"),
    ([ 8,3,3,3], [ 96_5945, 153_9624, 153_9624,  96_5945], "8.3.3.3"),
  ##([ 6,4,3,3], [ 90_0000, 174_7356, 109_4712, 144_7356], "6.4.3.3 A"),
    ([ 6,4,3,3], [110_9052, 159_0948, 138_1897, 138_1897], "6.4.3.3 B"),
  ##([ 6,3,4,3], [109_4712, 164_2068, 125_2644, 141_0576], "6.3.4.3"),
    ([ 6,3,3,3], [ 98_8994, 145_2219, 145_2219,  98_8994], "6.3.3.3"),
  ##([ 5,5,3,3], [126_8699, 142_6226, 158_3754, 142_6226], "5.5.3.3 A"),
  ##([ 5,5,3,3], [116_5651, 153_9424, 138_1897, 153_9424], "5.5.3.3 B"),
  ##([ 5,5,3,3], [ 63_4349, 171_3411,  70_5288, 171_3411], "5.5.3.3 C"),
  ##([ 5,3,5,3], [142_6226, 142_6226, 142_6226, 142_6226], "5.3.5.3"),
  ##([ 5,4,4,3], [153_4349, 144_0000, 169_1877, 142_6226], "5.4.4.3 A"),
  ##([ 5,4,4,3], [148_2825, 153_4349, 159_0948, 153_9424], "5.4.4.3 B"),
  ##([ 5,4,3,4], [148_2825, 159_0948, 159_0948, 148_2825], "5.4.3.4"),
  ##([ 5,4,3,3], [ 95_1524, 159_0948, 116_5651, 142_6226], "5.4.3.3 A"),
  ##([ 5,4,3,3], [ 90_0000, 162_7356, 109_4712, 144_7356], "5.4.3.3 B"),
  ##([ 5,3,4,3], [142_6226, 110_9052, 159_0948, 100_8123], "5.3.4.3"),
  ##([ 5,3,3,3], [100_8123, 138_1897, 138_1897, 100_8123], "5.3.3.3"),
  ##([ 4,4,4,3], [135_0000, 135_0000, 144_7356, 144_7356], "4.4.4.3 A"),
  ##([ 4,4,4,3], [120_0000, 144_7356, 125_2644, 160_5288], "4.4.4.3 B"),
  ##([ 4,4,4,3], [144_0000, 121_7175, 159_0948, 127_3774], "4.4.4.3 C"),
  ##([ 4,4,3,3], [ 60_0000, 160_5288,  70_5288, 160_5288], "4.4.3.3 A"),
  ##([ 4,4,3,3], [ 90_0000, 144_7356, 109_4712, 144_7356], "4.4.3.3 B"),
  ##([ 4,4,3,3], [108_0000, 127_3774, 138_1897, 127_3774], "4.4.3.3 C"),
  ##([ 4,4,3,3], [109_4712, 125_2644, 141_0576, 125_2644], "4.4.3.3 D"),
  ##([ 4,4,3,3], [ 63_4349, 159_0948,  74_7547, 159_0948], "4.4.3.3 E"),
    ([ 4,4,3,3], [117_0190, 109_5240, 159_8924, 109_5240], "4.4.3.3 F"),
    ([ 4,4,3,3], [ 72_9730, 154_7223,  86_7268, 154_7223], "4.4.3.3 G"),
    ([ 4,4,3,3], [102_5238, 133_9728, 128_4960, 133_9728], "4.4.3.3 H"),
    ([ 4,4,3,3], [100_1939, 136_3359, 124_7019, 136_3359], "4.4.3.3 I"),
  ##([ 4,3,4,3], [125_2644, 125_2644, 125_2644, 125_2644], "4.3.4.3 A"),
  ##([ 4,3,4,3], [ 90_0000, 150_0000,  90_0000, 150_0000], "4.3.4.3 B"),
  ##([ 4,3,4,3], [144_7356,  99_7356, 144_7356,  99_7356], "4.3.4.3 C"),
  ##([ 4,3,4,3], [159_0948,  69_0948, 159_0948,  69_0948], "4.3.4.3 D"),
    ([ 4,3,3,3], [103_8362, 127_5516, 127_5516, 103_8362], "4.3.3.3 A"),
  ##([ 4,3,3,3], [ 90_0000, 144_7356, 109_4712, 114_7356], "4.3.3.3 B"),
    ([ 4,3,3,3], [ 97_4555, 135_9915, 118_8922, 109_5240], "4.3.3.3 C"),
  ##([ 3,3,3,3], [109_4712, 109_4712, 109_4712, 109_4712], "3.3.3.3 A"),
  ##([ 3,3,3,3], [ 70_5288, 141_0576,  70_5288, 141_0576], "3.3.3.3 B"),
  ##([ 3,3,3,3], [138_1897,  74_7547, 138_1897,  74_7547], "3.3.3.3 C"),
    ([ 3,3,3,3], [ 96_1983, 121_7432,  96_1983, 121_7432], "3.3.3.3 D"),
    ([ 3,3,3,3], [ 86_7268, 129_4446,  86_7268, 129_4446], "3.3.3.3 E"),
    ([5,3,3,3,3], [152_9299, 164_1754, 164_1754, 164_1754, 152_9299], "5.3.3.3.3 A"),
  ##([5,3,3,3,3], [142_6226, 174_4343, 159_1865, 159_1865, 158_6816], "5.3.3.3.3 B"),
    ([4,3,3,3,3], [142_9834, 153_2346, 153_2346, 153_2346, 142_9834], "4.3.3.3.3 A"),
  ##([4,3,3,3,3], [125_2644, 169_4282, 145_2219, 145_2219, 153_6350], "4.3.3.3.3 B"),
  ##([4,3,3,3,3], [144_7356, 151_3301, 153_9624, 153_9624, 141_5945], "4.3.3.3.3 C"),
  ##([4,3,3,3,3], [159_0948, 132_6240, 159_1865, 159_1865, 126_9641], "4.3.3.3.3 D"),
    ([4,3,3,3,3], [145_4406, 144_1436, 164_2574, 144_1436, 145_4406], "4.3.3.3.3 E"),
  ##([4,3,3,3,3], [171_7546, 109_4712, 164_2596, 159_8924, 109_5240], "4.3.3.3.3 F"),
    ([4,3,3,3,3], [137_2401, 143_7383, 171_6457, 129_4446, 154_7223], "4.3.3.3.3 G"),
    ([4,3,3,3,3], [152_9756, 141_3411, 157_1481, 157_1481, 133_9728], "4.3.3.3.3 H"),
    ([4,3,3,3,3], [154_4188, 133_5912, 166_8114, 148_4340, 136_3359], "4.3.3.3.3 I"),
  ##([3,3,3,3,3], [138_1897, 138_1897, 138_1897, 138_1897, 138_1897], "3.3.3.3.3 A"),
  ##([3,3,3,3,3], [109_4712, 158_5718, 127_5516, 127_5516, 158_5718], "3.3.3.3.3 B"),
  ##([3,3,3,3,3], [109_4712, 144_7356, 144_7356, 109_4712, 169_4712], "3.3.3.3.3 C"),
    ([3,3,3,3,3], [ 96_1983, 166_4406, 121_7432, 121_7432, 166_4406], "3.3.3.3.3 D"),
    ([3,3,3,3,3], [164_2574, 114_6452, 144_1436, 144_1436, 114_6452], "3.3.3.3.3 E"),
    ([3,3,3,3,3], [159_8924, 118_8922, 143_4787, 143_4787, 118_8922], "3.3.3.3.3 F"),
    ([3,3,3,3,3], [131_4416, 143_4787, 135_9915, 135_9915, 143_4787], "3.3.3.3.3 G"),
  ##([3,3,3,3,3], [152_1911, 109_4712, 164_2596, 118_8922, 135_9915], "3.3.3.3.3 H"),
    ([3,3,3,3,3], [ 86_7268, 171_6457, 117_3556, 117_3556, 171_6457], "3.3.3.3.3 I"),
    ([3,3,3,3,3], [161_4828, 117_3556, 143_7383, 143_7383, 117_3556], "3.3.3.3.3 J"),
    ([3,3,3,3,3], [111_7348, 157_1481, 128_4960, 128_4960, 157_1481], "3.3.3.3.3 K"),
    ([3,3,3,3,3], [149_5648, 128_4960, 141_3411, 141_3411, 128_4960], "3.3.3.3.3 L"),
    ([3,3,3,3,3], [124_7019, 148_4340, 133_5912, 133_5912, 148_4340], "3.3.3.3.3 M"),
]
# The verf angles are given as integers (multiplying 0.0001 degrees),
# rather than strings, to make the above list slightly more compact and
# easier for humans to read.  Compare these integers to the angles in
# the sorted list, and replace them accordingly:
for i in range(len(verfs_3_12)):
    (edge_lengths, vertex_angles, name) = verfs_3_12[i]
    for j in range(len(vertex_angles)):
        integer_angle = vertex_angles[j]
        for x in sorted_angles:
            if integer_angle in x.value<<4:  # (Decimal shift 4 digits)
                vertex_angles[j] = x
                break
        else:
            raise AssertionError("verf " + name
                                 + " angle " + str(integer_angle)
                                 + " not found in list")
    # Convert tuple to Polygon:
    verfs_3_12[i] = Polygon(edge_lengths, vertex_angles, name)
# The "edge lengths" are also given as integers (e.g. 5 represents the
# pentagon's angle, 108 degrees).  This is no problem; they won't be
# used for any calculations, except comparisons.

for polygon in verfs_3_12:
    distinct_flags = []
    for flag1 in polygon.flags:
        # Check for duplicates, considering symmetries of the polygon.
        for flag2 in distinct_flags:
            if flag1.polygon_is_isomorphic(flag2):
                break
        else:
            distinct_flags.append(flag1)
    polygon.distinct_flags = distinct_flags

verfs_13_ = [
    Polygon([13,4,4], [pr_, pr13_17, pr_], "[13,17].4.4"),
    Polygon([18,4,4], [pr_, pr18_  , pr_],  "[18+].4.4"),
    Polygon([13,3,3,3], [ap_13_17, ap13_17, ap13_17, ap_13_17], "[13,17].3.3.3"),
    Polygon([18,3,3,3], [ap_18_  , ap18_  , ap18_  , ap_18_  ],  "[18+].3.3.3"),
]

for polygon in verfs_13_:
    flags = polygon.flags
    # The function 'polygon_is_isomorphic' doesn't work with just
    # intervals; it needs exact equality.  So we can't use it here.
    # But, for n.4.4 and n.3.3.3, we know flags 0 and 1 should be
    # equivalent by mirror symmetry.
    if len(flags) == 6:
        polygon.distinct_flags = [flags[0], flags[2], flags[3]]
    else:
        polygon.distinct_flags = [flags[0], flags[2], flags[3], flags[4]]

all_verfs = verfs_3_12 + verfs_13_
d = all_verfs_dict = {polygon.name: polygon for polygon in all_verfs}


# (An oriented polygon is a pair (polygon, flag_in_polygon).)
# Make a dictionary indexed by "edge length".  Each entry in the
# dictionary should be a list of oriented polygons, sorted by angle.
oriented_verfs = {}
for n in [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 18]:
    oriented_verfs[n] = []
for polygon in all_verfs:
    for flag in polygon.distinct_flags:
        n = flag.dyme[0]
        oriented_verfs[n].append((polygon, flag))

for n in [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 18]:
    oriented_verfs[n].sort(
        key=(lambda t: integer_infimum(t[1].dyme[1].value, 5)))
        # t == (polygon, flag)
        # flag.dyme == [edge_length, vertex_angle]
        # angle.value has no more than 5 digits after the decimal point.



After running this, we're ready to construct multihedra using the polygons. (I'm using IDLE, in case it matters.)


Let's start with the verf of snid:

Code: Select all

>>> m = Multihedron(d["5.3.3.3.3 A"])
>>> polyhedra = m.build_polyhedra(oriented_verfs)


Here's the output; in this case there are only two polyhedra:

Code: Select all

>>> for p in polyhedra:
...     print(p)
...

========================================================================
5.3.3.3.3 A
  0:  [1, 9, 34];  [5, 152.92990~5, 108.0000~28] 
  1:  [0, 2, 35];  [5, 152.92990~5, 108.0000~28] 
  2:  [3, 1, 15];  [3, 152.92990~5, 90.00000~39] 
  3:  [2, 4, 14];  [3, 164.17540~5, 90.00000~39] 
  4:  [5, 3, 20];  [3, 164.17540~5, 90.00000~39] 
  5:  [4, 6, 21];  [3, 164.17540~5, 90.00000~39] 
  6:  [7, 5, 26];  [3, 164.17540~5, 90.00000~39] 
  7:  [6, 8, 27];  [3, 164.17540~5, 90.00000~39] 
  8:  [9, 7, 33];  [3, 164.17540~5, 90.00000~39] 
__9:__[8, 0, 32];__[3, 152.92990~5, 90.00000~39]__
4.4.3
  10:  [11, 15, 40];  [4, 90.00000~5, 166.2125~15] 
  11:  [10, 12, 41];  [4, 60.00000~5, 166.2125~15] 
  12:  [13, 11, 18];  [4, 60.00000~5, 164.175398~81] 
  13:  [12, 14, 19];  [4, 90.00000~5, 164.175398~81] 
  14:  [15, 13, 3];  [3, 90.00000~5, 90.00000~39] 
__15:__[14, 10, 2];__[3, 90.00000~5, 90.00000~39]__
4.4.3
  16:  [17, 21, 25];  [4, 90.00000~5, 164.175398~81] 
  17:  [16, 18, 24];  [4, 60.00000~5, 164.175398~81] 
  18:  [19, 17, 12];  [4, 60.00000~5, 164.175398~81] 
  19:  [18, 20, 13];  [4, 90.00000~5, 164.175398~81] 
  20:  [21, 19, 4];  [3, 90.00000~5, 90.00000~39] 
__21:__[20, 16, 5];__[3, 90.00000~5, 90.00000~39]__
4.4.3
  22:  [23, 27, 28];  [4, 90.00000~5, 164.175398~81] 
  23:  [22, 24, 29];  [4, 60.00000~5, 164.175398~81] 
  24:  [25, 23, 17];  [4, 60.00000~5, 164.175398~81] 
  25:  [24, 26, 16];  [4, 90.00000~5, 164.175398~81] 
  26:  [27, 25, 6];  [3, 90.00000~5, 90.00000~39] 
__27:__[26, 22, 7];__[3, 90.00000~5, 90.00000~39]__
4.4.3
  28:  [29, 33, 22];  [4, 90.00000~5, 164.175398~81] 
  29:  [28, 30, 23];  [4, 60.00000~5, 164.175398~81] 
  30:  [31, 29, 46];  [4, 60.00000~5, 166.215~52] 
  31:  [30, 32, 47];  [4, 90.00000~5, 166.215~52] 
  32:  [33, 31, 9];  [3, 90.00000~5, 90.00000~39] 
__33:__[32, 28, 8];__[3, 90.00000~5, 90.00000~39]__
5.3.3
  34:  [35, 39, 0];  [5, 37.37740~5, 108.0000~28] 
  35:  [34, 36, 1];  [5, 37.37740~5, 108.0000~28] 
  36:  [37, 35, 45];  [3, 37.37740~5, 157.7610~29] 
  37:  [36, 38, 44];  [3, 138.18970~5, 157.7610~16] 
  38:  [39, 37, 49];  [3, 138.18970~5, 157.7610~16] 
__39:__[38, 34, 48];__[3, 37.37740~5, 157.7610~16]__
4.3.3
  40:  [41, 45, 10];  [4, 54.73560~5, 166.2125~15] 
  41:  [40, 42, 11];  [4, 54.73560~5, 166.2125~15] 
  42:  [43, 41, 51];  [3, 54.73560~5, 164.476~93] 
  43:  [42, 44, 50];  [3, 109.47120~5, 164.47743~70] 
  44:  [45, 43, 37];  [3, 109.47120~5, 157.7610~16] 
__45:__[44, 40, 36];__[3, 54.73560~5, 157.7610~29]__
4.3.3
  46:  [47, 51, 30];  [4, 54.73560~5, 166.215~52] 
  47:  [46, 48, 31];  [4, 54.73560~5, 166.215~52] 
  48:  [49, 47, 39];  [3, 54.73560~5, 157.7610~16] 
  49:  [48, 50, 38];  [3, 109.47120~5, 157.7610~16] 
  50:  [51, 49, 43];  [3, 109.47120~5, 164.47743~70] 
__51:__[50, 46, 42];__[3, 54.73560~5, 164.476~93]__

========================================================================
5.3.3.3.3 A
  0:  [1, 9, 34];  [5, 152.92990~5, 90.00000~23] 
  1:  [0, 2, 35];  [5, 152.92990~5, 90.00000~23] 
  2:  [3, 1, 15];  [3, 152.92990~5, 90.00000~23] 
  3:  [2, 4, 14];  [3, 164.17540~5, 90.00000~39] 
  4:  [5, 3, 20];  [3, 164.17540~5, 90.00000~39] 
  5:  [4, 6, 21];  [3, 164.17540~5, 90.00000~39] 
  6:  [7, 5, 26];  [3, 164.17540~5, 90.00000~39] 
  7:  [6, 8, 27];  [3, 164.17540~5, 90.00000~39] 
  8:  [9, 7, 33];  [3, 164.17540~5, 90.00000~39] 
__9:__[8, 0, 32];__[3, 152.92990~5, 90.00000~23]__
4.4.3
  10:  [11, 15, 36];  [4, 90.00000~5, 152.929899~93] 
  11:  [10, 12, 37];  [4, 60.00000~5, 152.929899~93] 
  12:  [13, 11, 18];  [4, 60.00000~5, 164.175398~81] 
  13:  [12, 14, 19];  [4, 90.00000~5, 164.175398~81] 
  14:  [15, 13, 3];  [3, 90.00000~5, 90.00000~39] 
__15:__[14, 10, 2];__[3, 90.00000~5, 90.00000~23]__
4.4.3
  16:  [17, 21, 25];  [4, 90.00000~5, 164.175398~81] 
  17:  [16, 18, 24];  [4, 60.00000~5, 164.175398~81] 
  18:  [19, 17, 12];  [4, 60.00000~5, 164.175398~81] 
  19:  [18, 20, 13];  [4, 90.00000~5, 164.175398~81] 
  20:  [21, 19, 4];  [3, 90.00000~5, 90.00000~39] 
__21:__[20, 16, 5];__[3, 90.00000~5, 90.00000~39]__
4.4.3
  22:  [23, 27, 28];  [4, 90.00000~5, 164.175398~81] 
  23:  [22, 24, 29];  [4, 60.00000~5, 164.175398~81] 
  24:  [25, 23, 17];  [4, 60.00000~5, 164.175398~81] 
  25:  [24, 26, 16];  [4, 90.00000~5, 164.175398~81] 
  26:  [27, 25, 6];  [3, 90.00000~5, 90.00000~39] 
__27:__[26, 22, 7];__[3, 90.00000~5, 90.00000~39]__
4.4.3
  28:  [29, 33, 22];  [4, 90.00000~5, 164.175398~81] 
  29:  [28, 30, 23];  [4, 60.00000~5, 164.175398~81] 
  30:  [31, 29, 38];  [4, 60.00000~5, 152.930~14] 
  31:  [30, 32, 39];  [4, 90.00000~5, 152.929899~93] 
  32:  [33, 31, 9];  [3, 90.00000~5, 90.00000~23] 
__33:__[32, 28, 8];__[3, 90.00000~5, 90.00000~39]__
5.4.4
  34:  [35, 39, 0];  [5, 90.00000~5, 90.00000~23] 
  35:  [34, 36, 1];  [5, 90.00000~5, 90.00000~23] 
  36:  [37, 35, 10];  [4, 90.00000~5, 152.929899~93] 
  37:  [36, 38, 11];  [4, 108.00000~5, 152.929899~93] 
  38:  [39, 37, 30];  [4, 108.00000~5, 152.930~14] 
__39:__[38, 34, 31];__[4, 90.00000~5, 152.929899~93]__



Now how do we interpret this jumble of numbers? Each row represents a flag. The flags are grouped according to the faces they're in (and the type of face is shown at the top of the group). The leftmost number in the row is the flag's index (just an arbitrary identifier). The next three numbers in the row are the indices of the adjacent flags. The last three numbers in the row are the edge length, and the digrammal angle, and the dihedral angle. Actually the edge length is 180° - 360°/n, where n is the number shown. And the tilde (as explained elsewhere) indicates an error bound in ULPs.

I'll draw three images below, based on the printed numbers above. The first image is the original multihedron (just one face), and of course the other two are the polyhedra. The edge lengths are encoded by tick-marks; a 60° edge is unmarked, a 90° edge has 1 mark, and a 108° edge has 2 marks.

[missing image]
[missing image]
[missing image]

Of course we already knew these two polyhedra existed. And 'multihedron.py' doesn't even guarantee that they do exist. But it does guarantee that these two are the only possibilities involving the snid verf.


Let's continue analyzing snid. We see that '5.3.3.3.3 A' must attach to four '4.4.3' faces, and these must share their '4' edges, thus meeting at a vertex not in '5.3.3.3.3 A'. This 3D result translates to a 4D result, that every trigonal face of snid must attach to a trigonal prism (or J7, J26, J49), and these must share their square faces, thus meeting at edges not in snid. Consider the far vertex Y of such an edge, and let X be the near vertex (the vertex in snid). Currently the verf at Y has four '4.4.3' faces and no other faces. We also know the dihedral angles between these faces, i.e. the dichoral angles between trigonal prisms, since they were determined at X. Now let's see what the program has to say about the verf at Y (after we construct the multihedron).

[missing image]

Code: Select all

>>> m1 = Multihedron(d["4.4.3"])
>>> m2 = m1.add_face(3, d["4.4.3"], 0)
>>> m3 = m2.add_face(9, d["4.4.3"], 0)
>>> m4 = m3.add_face(15, d["4.4.3"], 0)
>>> # p is the last polyhedron which was printed before.  (The relevant
>>> # angles are the same in the two polyhedra.  Either one will work.)
>>> flag1 = p.flags[11]  # In the verf at X
>>> flag2 = m4.flags[1]  # In the verf at Y
>>> for (flag3, flag4) in zip(flag1.verf_circuit(), flag2.verf_circuit()):
...     # Only complete edges should have dihedral angles.
...     if flag4.next[2] is not None:
...         # Transfer the angle from X to Y.
...         a = flag3.dyme[2]
...         flag4.dyme[2] = flag4.next[0].dyme[2] = a
...
>>> # Let's try a more compact form of output, one line per polyhedron.
>>> polyhedra = m4.build_polyhedra(oriented_verfs, verbosity=3)

4.4.3  4.4.3  4.4.3  4.4.3  4.3.3  4.3.3  5.3.3  5.3.3.3.3 A
4.4.3  4.4.3  4.4.3  4.4.3  5.4.4  5.3.3.3.3 A


Thus there are only two possible polyhedra. We didn't even need to see the types of faces; we already knew that the same two polyhedra from X would also work at Y, and the program says there are no more than two, so they must be exactly these two. In particular, there must be another snid at Y. So we have the outline of a snid prism; we have all of the trigonal prisms and both of the snids; only the pentagonal prisms are missing. And the output polyhedra have either '5.3.3' and '4.3.3' or '5.4.4' filling the gap. Since '5.4.4' is the verf of nothing but a pentagonal prism (or J9, J52, J53), and similarly for '5.3.3' and '4.3.3', the gap in 4D must be filled by either a pentagonal prism or a bunch of pentagonal and square pyramids.

In conclusion, there are no crown jewel polychora involving snid. There is nothing but the snid prism, and its augmentations by pentagonal prism pyramids.


I've already gone ahead along these lines, and found that there are no crown jewels involving snic or J84-90 either. And this is independent of my claim from a few years ago, that there are no crown jewel polychora with large hedra. But if we also accept that claim, it follows that any (hard) crown jewels must have interesting verfs, not interesting faces.

[quickfur]

[...]
Thus there are only two possible polyhedra. We didn't even need to see the types of faces; we already knew that the same two polyhedra from X would also work at Y, and the program says there are no more than two, so they must be exactly these two. In particular, there must be another snid at Y. So we have the outline of a snid prism; we have all of the trigonal prisms and both of the snids; only the pentagonal prisms are missing. And the output polyhedra have either '5.3.3' and '4.3.3' or '5.4.4' filling the gap. Since '5.4.4' is the verf of nothing but a pentagonal prism (or J9, J52, J53), and similarly for '5.3.3' and '4.3.3', the gap in 4D must be filled by either a pentagonal prism or a bunch of pentagonal and square pyramids.

In conclusion, there are no crown jewel polychora involving snid. There is nothing but the snid prism, and its augmentations by pentagonal prism pyramids.


I've already gone ahead along these lines, and found that there are no crown jewels involving snic or J84-90 either. And this is independent of my claim from a few years ago, that there are no crown jewel polychora with large hedra. But if we also accept that claim, it follows that any (hard) crown jewels must have interesting verfs, not interesting faces.

Wow, so now we know for sure that J84-90 do not generate any 4D CRFs besides their prisms? That's a pretty major result (if a negative one), if this is indeed the case. And here I was hoping that J86 would somehow crop up in a crown jewel. Guess that's out the window now.

So it seems the only thing that remains now is to find crown jewels consisting of "boring" facets like prisms and pyramids, put together in unusual ways. Now I'm starting to wonder if these would be limited too, because the spherical configuration of 4D verfs imposes many restrictions on what can fit around it.

[mr_e_man]

I too was hoping to find something new.

The 6,8,10-gon antiprisms could qualify as interesting faces. They are constructible (they don't require cube roots or higher polynomials), but I don't know of any polychora involving them, other than antiprismatic prisms and biantiprismatic wedges. We should look into that.


About the program's output:
Maybe each multihedron should keep a record of the operations that were done to create it. Printing that would be less verbose than the current method, but would still give enough information, more than just the types of faces.
And I imagine that outputting the whole tree wouldn't be necessary with this approach. Or, we could consider that record as a path through the tree.

...Okay, now it's implemented as an 'ancestry' attribute.

multihedron.py

Code: Select all

# Copyright (c) 2025 Anonymous.
#
# Permission to use, copy, modify, and/or distribute this software for
# any purpose with or without fee is hereby granted.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
# WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS.  IN NO EVENT SHALL THE
# AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
# DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
# PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
# TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
# PERFORMANCE OF THIS SOFTWARE.

"""Building convex polyhedra with given types of faces."""

from copy import deepcopy

import decball
from decball import sqrt, sind, cosd, cisd, acosd
from exact_angle import Angle


class Nonconvexity(Exception):
    pass

class PrimaryNonconvexity(Nonconvexity):
    pass

class SecondaryNonconvexity(Nonconvexity):
    pass


class Flag:
    """A triple (vertex, edge, face), all three being incident.

    An adjacent flag, self.next[i], has a different i-dimensional
    surtope, and all other surtopes the same as in self.

    A dyad measure, self.dyme[i], is the length or angle separating self
    from an adjacent flag.  This includes edge lengths (between two
    vertices), digrammal angles (between two edges), and dihedral angles
    (between two faces).

    The surtopes (i.e. vertices, edges, faces) usually are not directly
    referred to.  Instead, a surtope can be considered as the set of all
    flags that involve the surtope.
    """

    def __init__(self):
        self.next = [None]*3
        self.dyme = [None]*3

    def face_circuit(self, reverse=False, include_self=True):
        """Go around the polygon using next[0] and next[1].

        The first flag is self, the second flag is self.next[0], the
        third flag is self.next[0].next[1], the fourth flag is
        self.next[0].next[1].next[0], and so on.  It stops when it
        reaches self again (in which case the last flag is self.next[1])
        or it reaches None (in which case the polygon is incomplete, and
        self.next[1] either is never reached or is None).

        The previous paragraph assumes default arguments.  If
        reverse=True, then next[1] is done before next[0].

        This is a generator function, used to make 'for' loops.
        """
        if include_self:
            yield self
        (i, j) = (0, 1)
        if reverse:
            (i, j) = (j, i)
        flag = self
        while True:
            flag = flag.next[i]
            if flag is None:  # (flag can't be self, by orientability)
                return
            yield flag
            flag = flag.next[j]
            if flag is None or flag is self:
                return
            yield flag

    def verf_circuit(self, reverse=False, include_self=True):
        """Go around the vertex using next[1] and next[2].

        The first flag is self, the second flag is self.next[1], the
        third flag is self.next[1].next[2], the fourth flag is
        self.next[1].next[2].next[1], and so on.  It stops when it
        reaches self again (in which case the last flag is self.next[2])
        or it reaches None (in which case the vertex figure is
        incomplete, and self.next[2] either is never reached or is
        None).

        The previous paragraph assumes default arguments.  If
        reverse=True, then next[2] is done before next[1].

        This is a generator function, used to make 'for' loops.
        """
        if include_self:
            yield self
        (i, j) = (1, 2)
        if reverse:
            (i, j) = (j, i)
        flag = self
        while True:
            flag = flag.next[i]
            if flag is None:  # (flag can't be self, by orientability)
                return
            yield flag
            flag = flag.next[j]
            if flag is None or flag is self:
                return
            yield flag

    def polygon_is_isomorphic(self, flag):
        """Return True if the two flags are in equivalent positions in
        congruent polygons.

        If both flags are in the same polygon, this is a way to find its
        symmetries.

        Both polygons are assumed to be complete, and to have known edge
        lengths and vertex angles.  If they're known only by intervals,
        this function will probably return False, because it uses '=='.
        """
        # Check for equal lengths and angles.
        for (flag1, flag2) in zip(self.face_circuit(), flag.face_circuit()):
            for (m1, m2) in zip(flag1.dyme[:2], flag2.dyme[:2]):
                if m1 is None or m2 is None:
                    raise ValueError("polygon dyad measures must be specified")
                if not m1 == m2:  # (Don't replace this with 'm1 != m2')
                    return False
        # Check for equal valence.  If they're unequal, then 'zip'
        # finishes with the smaller polygon, and the larger polygon
        # isn't finished.
        if flag1.next[1] is not self or flag2.next[1] is not flag:
            return False
        return True

    def solve_verf(self):
        """Calculate all the dihedral angles at this vertex.

        Return True if successful, False if the vertex figure is
        incomplete or underdetermined (with more than 3 unknown (None)
        angles).

        This modifies the dihedrals (only at this vertex, not the
        adjacent vertices), replacing all occurrences of None with
        values calculated by trigonometry, assuming convexity.

        The digrammals and dihedrals are assumed to be of Angle type.
        Furthermore, the digrammals are assumed to have cosine and sine
        specified (not None) and of Interval type.

        Warning:  If this function is interrupted by an Exception, the
        vertex may be "broken"; the flag adjacencies may be changed.  If
        you want to preserve the original structure, make a deepcopy
        first.
        """
        unknown_count = 0
        known_flag = None
        for (i, flag) in enumerate(self.verf_circuit()):
            if flag.dyme[1] is None:
                raise ValueError("digrammals must be specified")
            if flag.dyme[2] is None:
                unknown_count += 1
            elif known_flag is None and i%2 == 0:
                known_flag = flag
        if flag.next[2] is None:
            # The vertex figure is incomplete.
            return False
        valence = (i+1)//2
        unknown_count //= 2
        if valence < 3:
            raise ValueError("vertex should have at least 3 edges")
        elif valence == 3:
            flag0 = self
            flag1 = flag0.next[1]
            flag2 = flag1.next[2]
            flag3 = flag2.next[1]
            flag4 = flag3.next[2]
            flag5 = flag4.next[1]
            assert flag5.next[2] is flag0
            # Spherical triangle's edge lengths and vertex angles:
            l0 = flag0.dyme[1]
            a0 = flag2.dyme[2]
            l1 = flag2.dyme[1]
            a1 = flag4.dyme[2]
            l2 = flag4.dyme[1]
            a2 = flag0.dyme[2]
            if l0 + l1 < l2 or l1 + l2 < l0 or l2 + l0 < l1:
                raise Nonconvexity("triangle inequality is violated")
            # Solve the triangle using the law of cosines.
            (cos0, sin0) = (l0.cosine, l0.sine)
            (cos1, sin1) = (l1.cosine, l1.sine)
            (cos2, sin2) = (l2.cosine, l2.sine)
            if None in (cos0, sin0, cos1, sin1, cos2, sin2):
                raise ValueError("Angle must have a specified cosine and sine")
            cos0_ = (cos2 - cos0*cos1)/(sin0*sin1)  # (int/int = float
            cos1_ = (cos0 - cos1*cos2)/(sin1*sin2)  #  must be avoided)
            cos2_ = (cos1 - cos2*cos0)/(sin2*sin0)
            a0_ = Angle(value=acosd(cos0_), cosine=cos0_)
            a1_ = Angle(value=acosd(cos1_), cosine=cos1_)
            a2_ = Angle(value=acosd(cos2_), cosine=cos2_)
            # Compare the calculated angles to any previously known
            # angles.  Keep the previous ones if they're more precise.
            if a0 is not None:
                if a0 != a0_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a0.value in a0_.value:  # (Interval containment)
                    a0_ = a0
            if a1 is not None:
                if a1 != a1_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a1.value in a1_.value:
                    a1_ = a1
            if a2 is not None:
                if a2 != a2_:
                    raise Nonconvexity("dihedrals are inconsistent")
                elif a2.value in a2_.value:
                    a2_ = a2
            flag1.dyme[2] = flag2.dyme[2] = a0_
            flag3.dyme[2] = flag4.dyme[2] = a1_
            flag5.dyme[2] = flag0.dyme[2] = a2_
            return True
        # Solve the spherical polygon by recursion, cutting off a
        # triangle to make a simpler polygon.
        # Find some known angle.
        if unknown_count > 3:
            return False
        a = known_flag.dyme[2]
        cos_a = a.cosine
        if cos_a is None:
            (cos_a, sin_a) = cisd(a.value)
            a.cosine = cos_a
            a.sine = sin_a
        flag5 = known_flag
        flag4 = flag5.next[2]
        flag3 = flag4.next[1]
        flag0 = flag3.next[2]
        flag6 = flag5.next[1]
        flag9 = flag6.next[2]
        prev_angle = flag3.dyme[2]
        prev_length = flag4.dyme[1]
        next_length = flag5.dyme[1]
        next_angle = flag6.dyme[2]
        # Use the law of cosines to find the length of the cutting line.
        (cos_prev, sin_prev) = (prev_length.cosine, prev_length.sine)
        (cos_next, sin_next) = (next_length.cosine, next_length.sine)
        if None in (cos_prev, sin_prev, cos_next, sin_next):
            raise ValueError("Angle must have a specified cosine and sine")
        cos_l = cos_prev*cos_next + sin_prev*sin_next*cos_a
        sin_l = sqrt(1 - cos_l**2)
        l = Angle(value=acosd(cos_l), cosine=cos_l, sine=sin_l)
        # The cutting line becomes an edge, with four new flags.
        flag1 = Flag()
        flag2 = Flag()
        flag7 = Flag()
        flag8 = Flag()
        flag1.next[1] = flag8;  flag8.next[1] = flag1
        flag2.next[1] = flag7;  flag7.next[1] = flag2
        flag1.dyme[1] = flag8.dyme[1] = l
        flag2.dyme[1] = flag7.dyme[1] = l
        # Cut.
        flag0.next[2] = flag1;  flag1.next[2] = flag0
        flag3.next[2] = flag2;  flag2.next[2] = flag3
        flag6.next[2] = flag7;  flag7.next[2] = flag6
        flag9.next[2] = flag8;  flag8.next[2] = flag9
        flag0.dyme[2] = flag3.dyme[2] = None
        flag6.dyme[2] = flag9.dyme[2] = None
        # Solve the triangle.
        cos2 = (cos_next - cos_prev*cos_l)/(sin_prev*sin_l)
        cos7 = (cos_prev - cos_next*cos_l)/(sin_next*sin_l)
        a2 = Angle(value=acosd(cos2), cosine=cos2)
        a7 = Angle(value=acosd(cos7), cosine=cos7)
        flag3.dyme[2] = flag2.dyme[2] = a2
        flag6.dyme[2] = flag7.dyme[2] = a7
        # Subtract the angles in the triangle from those in the original
        # polygon (if known).
        if prev_angle is not None:
            a1 = prev_angle - a2
            if a1 < 0:
                raise Nonconvexity
            flag0.dyme[2] = flag1.dyme[2] = a1
        if next_angle is not None:
            a8 = next_angle - a7
            if a8 < 0:
                raise Nonconvexity
            flag9.dyme[2] = flag8.dyme[2] = a8
        # Solve the remaining polygon.  (Certainly it can be solved.)
        success = flag1.solve_verf()
        assert success
        # Uncut.
        flag1.next[2] = flag8.next[2] = None
        flag2.next[2] = flag7.next[2] = None
        flag0.next[2] = flag3;  flag3.next[2] = flag0
        flag9.next[2] = flag6;  flag6.next[2] = flag9
        # Add the angles in the triangle to get those in the original
        # polygon (if unknown).
        if prev_angle is None:
            a1 = flag1.dyme[2]
            prev_angle = a1 + a2
            if prev_angle > 180:
                raise Nonconvexity
        flag0.dyme[2] = flag3.dyme[2] = prev_angle
        if next_angle is None:
            a8 = flag8.dyme[2]
            next_angle = a8 + a7
            if next_angle > 180:
                raise Nonconvexity
        flag9.dyme[2] = flag6.dyme[2] = next_angle
        return True

    def solve_all(self):
        """Calculate all dihedral angles, or as many as possible.

        Start by solving this vertex figure, and travel to adjacent
        vertices via edges with newly calculated dihedrals, and solve
        those vertex figures also, and so on.

        This modifies the dihedrals.  Nothing is returned.

        The digrammals and dihedrals are assumed to be of Angle type.
        Furthermore, the digrammals are assumed to have cosine and sine
        specified (not None) and of Interval type.

        Warning:  If this function is interrupted by an Exception, a
        vertex may be "broken"; the flag adjacencies may be changed.  If
        you want to preserve the original structure, make a deepcopy
        first.
        """
        success = self.solve_verf()
        if not success:
            return
        newly_solved = [self]  # (Use 'collections.deque' instead of 'list'?)
        while newly_solved:
            flag1 = newly_solved.pop(0)
            # Try all the edges around the vertex that was solved.
            for (i, flag2) in enumerate(flag1.verf_circuit()):
                # (Don't consider each edge twice.)
                if i%2 == 1:
                    continue
                a2 = flag2.dyme[2]
                assert a2 is not None  # (solve_verf was successful)
                # Go to the adjacent vertex.
                flag3 = flag2.next[0]
                a3 = flag3.dyme[2]
                if a3 is not None:
                    # The dihedral angle has been calculated before.
                    if a3 != a2:
                        raise Nonconvexity("dihedrals are inconsistent")
                    continue
                # The dihedral angle is new.  Transfer it to the
                # adjacent vertex.
                # But the interval may be uselessly wide, e.g. with a
                # radius greater than 10 degrees.
                b = a2.value
                if isinstance(b, int) or (isinstance(b, decball.Ball)
                        and b.get_radius_precision() + b.exp <= 1):
                    flag3.dyme[2] = flag3.next[2].dyme[2] = a2
                    success = flag3.solve_verf()
                    if success:
                        newly_solved.append(flag3)


class Polygon:
    """Needs no introduction.

    The (cyclic) order of the surtopes is:
    vertex0, edge0, vertex1, edge1, vertex2, edge2, ...

    The order of the flags is:
    flag0 = (vertex0, edge0), flag1 = (vertex1, edge0),
    flag2 = (vertex1, edge1), flag3 = (vertex2, edge1), ...

    Edge length 0 is between vertices 0 and 1, and vertex angle 0 is
    between edges 0 and 1, so length i is at edge i, but angle i is at
    vertex i+1.
    """

    def __init__(self, edge_lengths, vertex_angles, name=None):
        self.edge_lengths = list(edge_lengths)
        self.vertex_angles = list(vertex_angles)
        self.valence = len(edge_lengths)
        if self.valence != len(vertex_angles):
            raise ValueError("polygon must have equal numbers of "
                             + "vertices and edges")
        flags = []
        for i in range(2*self.valence):
            flags.append(Flag())
        for i in range(self.valence):
            j = 2*i
            flags[j].next[0] = flags[j+1];  flags[j+1].next[0] = flags[j]
            flags[j].next[1] = flags[j-1];  flags[j-1].next[1] = flags[j]
            flags[j].dyme[0] = flags[j+1].dyme[0] = edge_lengths[i]
            flags[j].dyme[1] = flags[j-1].dyme[1] = vertex_angles[i-1]
        self.flags = flags
        self.name = name

    def __str__(self):
        return ("Polygon(edge_lengths=" + str(self.edge_lengths)
                + ", vertex_angles=" + str(self.vertex_angles)
                + ")")


class Multihedron:
    """An incomplete polyhedron, made of complete polygons.

    In a complete polyhedron, every flag has flag.next[2] defined (not
    None).  In a multihedron, some flags may have flag.next[2] undefined
    (None).

    A Multihedron is initialized with a single face.  The given Polygon
    is copied (using deepcopy), so it can be re-used for other faces.

    When a Multihedron is printed, the format for each line is:
    flag_index :  adjacent_flag_indices ;  dyad_measures

    A Multihedron's ancestry is a string describing how to build the
    Multihedron.  This includes the face names, the operations 'je' and
    'af', and the relevant flag indices.  (However, evaluating the
    string as code is not valid, because a polygon's name is not the
    polygon itself.)
    """

    def __init__(self, polygon, name=None):
        if not isinstance(polygon, Polygon):
            raise TypeError
        face = deepcopy(polygon)
        self.faces = [face]
        # We'll need three different lists to update separately.
        self.flags = list(face.flags)  # Copy
        self.incomplete_flags = list(face.flags)  # Another copy
        for flag in self.flags:
            angle = flag.dyme[1]
            # Sines and cosines won't be needed for simple addition
            # (though they will be needed to calculate dihedrals).
            s = Angle(angle.value, angle.coefs, None, None)
            flag.vertex_digrammal_sum = s
        self.name = name
        self.ancestry = repr(face.name) + " "

    def __str__(self):
        all_flags = self.flags
        lines = []  # (Lines of text)
        if self.name is not None:
            lines.append(str(self.name))
        lines.append("=" * 72)
        i = 0
        for face in self.faces:
            if face.name is not None:
                lines.append(str(face.name))
            last = len(face.flags)-1
            for (j, flag1) in enumerate(face.flags):
                space = "__" if j == last else "  "
                assert i == all_flags.index(flag1)
                index_next = [None if flag2 is None else all_flags.index(flag2)
                              for flag2 in flag1.next]
                lines.append(space + str(i) + ":" + space
                             + str(index_next) + ";" + space
                             + str(flag1.dyme) + space)
                i += 1
        lines.append("")
        return "\n".join(lines)

    def join_edges(self, flag1):
        """Return multihedron with flag1 edge joined to opposing edge.

        It's assumed that flag1 is incomplete (so flag1.next[2] is
        None).  Let flag2 be the other incomplete flag at the same
        vertex.

        self is copied, not modified.  The new multihedron has (the
        copies of) flag1 and flag2 adjacent, sharing a vertex and an
        edge.

        The "flag" argument, if it's actually an integer, is interpreted
        as the index of a flag in the multihedron.

        "join_edges" may be abbreviated as "je".
        """
        if isinstance(flag1, int):
            index1 = flag1
            flag1 = self.flags[index1]
        else:
            index1 = self.flags.index(flag1)
        if flag1.next[2] is not None:
            raise ValueError("edge is already complete")
        for (i, flag2) in enumerate(flag1.verf_circuit()):
            pass
        valence = (i+1)//2
        if valence < 3:
            raise Nonconvexity("vertex would be monovalent or divalent")
        if flag2.dyme[0] != flag1.dyme[0]:
            raise Nonconvexity("edge lengths don't match")
        # The other ends of the edges would also be joined.
        flag1_0 = flag1.next[0]
        flag2_0 = flag2.next[0]
        # (flag1_0 is not flag2_0, by orientability.)
        for (i, flag3) in enumerate(flag1_0.verf_circuit()):
            pass
        adjacent_vertices_coincide = flag3 is flag2_0
        if adjacent_vertices_coincide:
            # The other ends of the edges are already a single vertex,
            # which will be complete.
            valence = (i+1)//2
            if valence < 3:
                raise Nonconvexity("vertex would be monovalent or divalent")
        else:
            new_sum = (flag1_0.vertex_digrammal_sum
                       + flag2_0.vertex_digrammal_sum)
            if new_sum >= 360:
                raise Nonconvexity("vertex angle deficit isn't positive")
            # (Don't join the two vertices if they're at different
            # locations in a single face?  That test is too complex.
            # And it may be redundant somehow.)
        # Join the edges; make a new multihedron.
        (multihedron, flag1_, flag2_) = deepcopy((self, flag1, flag2))
        flag1_0 = flag1_.next[0]
        flag2_0 = flag2_.next[0]
        flag1_.next[2] = flag2_;  flag2_.next[2] = flag1_
        flag1_0.next[2] = flag2_0;  flag2_0.next[2] = flag1_0
        for flag in (flag1_, flag2_, flag1_0, flag2_0):
            multihedron.incomplete_flags.remove(flag)
        if not adjacent_vertices_coincide:
            # Two vertices are joined.  Update the angle sums.
            for flag in flag1_0.verf_circuit():
                flag.vertex_digrammal_sum = new_sum
            for flag in flag2_0.verf_circuit():
                flag.vertex_digrammal_sum = new_sum
        # Calculate dihedral angles, if possible.
        flag1_.solve_all()
        if adjacent_vertices_coincide:
            flag1_0.solve_all()
        multihedron.ancestry += ".je(" + str(index1) + ")"
        return multihedron

    je = join_edges

    def add_face(self, flag1, polygon, flag2=0):
        """Return multihedron with new face at flag1.

        flag1 is in self, and flag2 is in polygon.  All four arguments
        are copied, not modified.  The new multihedron has (the copies
        of) flag1 and flag2 adjacent, sharing a vertex and an edge.

        The polygon is assumed to be disconnected from other polygons
        (so each flag.next[2] is None); it shouldn't be a face of a
        previously existing multihedron.

        A "flag" argument, if it's actually an integer, is interpreted
        as the index of a flag in the multihedron or the polygon.

        PrimaryNonconvexity is raised, if the sum of digrammal angles
        around the vertex is at least 360, according to interval
        arithmetic.  SecondaryNonconvexity is for the adjacent vertex at
        the same edge.  If one of these sums is at least 360 according
        to exact arithmetic, only Nonconvexity is raised.

        "add_face" may be abbreviated as "af".
        """
        if not isinstance(polygon, Polygon):
            raise TypeError
        if isinstance(flag1, int):
            index1 = flag1
            flag1 = self.flags[index1]
        else:
            index1 = self.flags.index(flag1)
        if isinstance(flag2, int):
            index2 = flag2
            flag2 = polygon.flags[index2]
        else:
            index2 = polygon.flags.index(flag2)
        if flag1.next[2] is not None:
            raise ValueError("edge is already complete")
        # Check the angle sum at this vertex.
        new_sum = flag1.vertex_digrammal_sum + flag2.dyme[1]
        if new_sum.value >= 360:  # (Interval comparison)
            raise PrimaryNonconvexity("vertex angle deficit isn't positive")
        elif new_sum == 360:
            raise Nonconvexity("vertex angle deficit isn't positive")
        # Check the angle sum at the adjacent vertex.
        flag1_0 = flag1.next[0]
        flag2_0 = flag2.next[0]
        new_sum_0 = flag1_0.vertex_digrammal_sum + flag2_0.dyme[1]
        if new_sum_0.value >= 360:
            raise SecondaryNonconvexity("vertex angle deficit isn't positive")
        elif new_sum_0 == 360:
            raise Nonconvexity("vertex angle deficit isn't positive")
        # Add the face; make a new multihedron.
        (multihedron, flag1_) = deepcopy((self, flag1))
        (face, flag2_) = deepcopy((polygon, flag2))
        multihedron.faces.append(face)
        multihedron.flags.extend(face.flags)
        multihedron.incomplete_flags.extend(face.flags)
        for flag in face.flags:
            angle = flag.dyme[1]
            s = Angle(angle.value, angle.coefs, None, None)
            flag.vertex_digrammal_sum = s
        flag1_0 = flag1_.next[0]
        flag2_0 = flag2_.next[0]
        # Connect the face.
        flag1_.next[2] = flag2_;  flag2_.next[2] = flag1_
        flag1_0.next[2] = flag2_0;  flag2_0.next[2] = flag1_0
        for flag in (flag1_, flag2_, flag1_0, flag2_0):
            multihedron.incomplete_flags.remove(flag)
        # Update the angle sums.
        for flag in flag1_.verf_circuit():
            flag.vertex_digrammal_sum = new_sum
        for flag in flag2_.verf_circuit():  # (Only two flags here)
            flag.vertex_digrammal_sum = new_sum
        for flag in flag1_0.verf_circuit():
            flag.vertex_digrammal_sum = new_sum_0
        for flag in flag2_0.verf_circuit():  # (Only two flags here)
            flag.vertex_digrammal_sum = new_sum_0
        multihedron.ancestry += (".af(" + str(index1) + ", "
                                 + repr(face.name) + ", " + str(index2) + ")")
        return multihedron

    af = add_face

    def build_polyhedra(self, oriented_polygons, verbosity=1):
        """Return list of all possible completions of self.

        The new faces being added to the multihedron are taken from
        oriented_polygons.  This argument should be a kind of 2D array:
        a dictionary indexed by edge length, where each entry is a list
        sorted by angle, where each entry is an oriented polygon.  And
        an oriented polygon is a tuple (polygon, flag_in_polygon).

        To be more precise, each list (in oriented_polygons) should be
        sorted by inf(angle.value), assuming each angle.value is an
        interval.  (The infimum of an interval [a,b] is a.)

        Each returned polyhedron may be printed as soon as it's found,
        not waiting for this function to finish.  This is determined by
        verbosity, as follows:

        verbosity=0:  Don't print anything
        verbosity=1:  Print a dot, as part of an extended ellipsis
        verbosity=2:  Print the number of faces in the polyhedron
        verbosity=3:  Print the names of all faces in the polyhedron
        verbosity=4:  Print the ancestry of the polyhedron
        verbosity=5:  Print the whole polyhedron, with flag data
        """
        if not self.incomplete_flags:
            # The multihedron is complete, i.e. a polyhedron.
            if verbosity == 1:
                print(".", end="")
            elif verbosity == 2:
                print(len(self.faces), end=" ")
            elif verbosity == 3:
                for face in self.faces:
                    print(face.name, end="  ")
                print("")
            elif verbosity == 4:
                print(self.ancestry, end="\n\n")
            elif verbosity == 5:
                print(self)
            return [self]
        sorted_flags = sorted(  # (Use 'heapq' module instead?)
            self.incomplete_flags,
            key=(lambda flag: (flag.vertex_digrammal_sum,
                               flag.next[0].vertex_digrammal_sum)),
            reverse=True)
        flag1 = sorted_flags[0]
        polyhedra = []
        try:
            multihedron = self.join_edges(flag1)
        except Nonconvexity:
            pass
        else:
            polyhedra += multihedron.build_polyhedra(oriented_polygons,
                                                     verbosity)
        edge_length = flag1.dyme[0]
        for (polygon, flag2) in oriented_polygons[edge_length]:
            try:
                multihedron = self.add_face(flag1, polygon, flag2)
            except PrimaryNonconvexity:
                # The polygon list is sorted; if one angle is too large,
                # then all later angles are also too large.
                break
            except Nonconvexity:
                continue
            polyhedra += multihedron.build_polyhedra(oriented_polygons,
                                                     verbosity)
        return polyhedra