/*****************************************************************************
* #ident "Id: main.c,v 3.27 2002-01-06 16:23:01+02 rl Exp "
* kaleido
*
* Kaleidoscopic construction of uniform polyhedra
* Copyright (c) 1991-2002 Dr. Zvi Har'El <rl@math.technion.ac.il>
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
*
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in
* the documentation and/or other materials provided with the
* distribution.
*
* 3. The end-user documentation included with the redistribution,
* if any, must include the following acknowledgment:
* "This product includes software developed by
* Dr. Zvi Har'El (http://www.math.technion.ac.il/~rl/)."
* Alternately, this acknowledgment may appear in the software itself,
* if and wherever such third-party acknowledgments normally appear.
*
* This software is provided 'as-is', without any express or implied
* warranty. In no event will the author be held liable for any
* damages arising from the use of this software.
*
* Author:
* Dr. Zvi Har'El,
* Deptartment of Mathematics,
* Technion, Israel Institue of Technology,
* Haifa 32000, Israel.
* E-Mail: rl@math.technion.ac.il
*
* ftp://ftp.math.technion.ac.il/kaleido/
* http://www.mathconsult.ch/showroom/unipoly/
*
* Adapted for xscreensaver by Jamie Zawinski <jwz@jwz.org> 25-Apr-2004
*
*****************************************************************************
*/
#ifdef HAVE_CONFIG_H
# include "config.h"
#endif
#include <math.h>
#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <stdlib.h>
#include <errno.h>
#include "polyhedra.h"
extern const char *progname;
#ifndef MAXLONG
#define MAXLONG 0x7FFFFFFF
#endif
#ifndef MAXDIGITS
#define MAXDIGITS 10 /* (int)log10((double)MAXLONG) + 1 */
#endif
#ifndef DBL_EPSILON
#define DBL_EPSILON 2.2204460492503131e-16
#endif
#define BIG_EPSILON 3e-2
#define AZ M_PI/7 /* axis azimuth */
#define EL M_PI/17 /* axis elevation */
#define Err(x) do {\
fprintf (stderr, "%s: %s\n", progname, (x)); \
exit (1); \
} while(0)
#define Free(lvalue) do {\
if (lvalue) {\
free((char*) lvalue);\
lvalue=0;\
}\
} while(0)
#define Matfree(lvalue,n) do {\
if (lvalue) \
matfree((char*) lvalue, n);\
lvalue=0;\
} while(0)
#define Malloc(lvalue,n,type) do {\
if (!(lvalue = (type*) calloc((n), sizeof(type)))) \
abort();\
} while(0)
#define Realloc(lvalue,n,type) do {\
if (!(lvalue = (type*) realloc(lvalue, (n) * sizeof(type)))) \
abort();\
} while(0)
#define Calloc(lvalue,n,type) do {\
if (!(lvalue = (type*) calloc(n, sizeof(type))))\
abort();\
} while(0)
#define Matalloc(lvalue,n,m,type) do {\
if (!(lvalue = (type**) matalloc(n, (m) * sizeof(type))))\
abort();\
} while(0)
#define Sprintfrac(lvalue,x) do {\
if (!(lvalue=sprintfrac(x)))\
return 0;\
} while(0)
#define numerator(x) (frac(x), frax.n)
#define denominator(x) (frac(x), frax.d)
#define compl(x) (frac(x), (double) frax.n / (frax.n-frax.d))
typedef struct {
double x, y, z;
} Vector;
typedef struct {
/* NOTE: some of the int's can be replaced by short's, char's,
or even bit fields, at the expense of readability!!!*/
int index; /* index to the standard list, the array uniform[] */
int N; /* number of faces types (atmost 5)*/
int M; /* vertex valency (may be big for dihedral polyhedra) */
int V; /* vertex count */
int E; /* edge count */
int F; /* face count */
int D; /* density */
int chi; /* Euler characteristic */
int g; /* order of symmetry group */
int K; /* symmetry type: D=2, T=3, O=4, I=5 */
int hemi;/* flag hemi polyhedron */
int onesided;/* flag onesided polyhedron */
int even; /* removed face in pqr| */
int *Fi; /* face counts by type (array N)*/
int *rot; /* vertex configuration (array M of 0..N-1) */
int *snub; /* snub triangle configuration (array M of 0..1) */
int *firstrot; /* temporary for vertex generation (array V) */
int *anti; /* temporary for direction of ideal vertices (array E) */
int *ftype; /* face types (array F) */
int **e; /* edges (matrix 2 x E of 0..V-1)*/
int **dual_e; /* dual edges (matrix 2 x E of 0..F-1)*/
int **incid; /* vertex-face incidence (matrix M x V of 0..F-1)*/
int **adj; /* vertex-vertex adjacency (matrix M x V of 0..V-1)*/
double p[4]; /* p, q and r; |=0 */
double minr; /* smallest nonzero inradius */
double gon; /* basis type for dihedral polyhedra */
double *n; /* number of side of a face of each type (array N) */
double *m; /* number of faces at a vertex of each type (array N) */
double *gamma; /* fundamental angles in radians (array N) */
char *polyform; /* printable Wythoff symbol */
char *config; /* printable vertex configuration */
char *group; /* printable group name */
char *name; /* name, standard or manifuctured */
char *dual_name; /* dual name, standard or manifuctured */
char *class;
char *dual_class;
Vector *v; /* vertex coordinates (array V) */
Vector *f; /* face coordinates (array F)*/
} Polyhedron;
typedef struct {
long n,d;
} Fraction;
static Polyhedron *polyalloc(void);
static Vector rotate(Vector vertex, Vector axis, double angle);
static Vector sum3(Vector a, Vector b, Vector c);
static Vector scale(double k, Vector a);
static Vector sum(Vector a, Vector b);
static Vector diff(Vector a, Vector b);
static Vector pole (double r, Vector a, Vector b, Vector c);
static Vector cross(Vector a, Vector b);
static double dot(Vector a, Vector b);
static int same(Vector a, Vector b, double epsilon);
static char *sprintfrac(double x);
static void frac(double x);
static void matfree(void *mat, int rows);
static void *matalloc(int rows, int row_size);
static Fraction frax;
static const struct {
char *Wythoff, *name, *dual, *group, *class, *dual_class;
short Coxeter, Wenninger;
} uniform[] = {
/****************************************************************************
* Dihedral Schwarz Triangles (D5 only)
***************************************************************************/
/* 0 */ {"2 5|2", "Pentagonal Prism",
"Pentagonal Dipyramid",
"Dihedral (D[1/5])",
"",
"",
0, 0},
/* 2 */ {"|2 2 5", "Pentagonal Antiprism",
"Pentagonal Deltohedron",
"Dihedral (D[1/5])",
"",
"",
0, 0},
/* (2 2 5/2) (D2/5) */
/* 4 */ {"2 5/2|2", "Pentagrammic Prism",
"Pentagrammic Dipyramid",
"Dihedral (D[2/5])",
"",
"",
0, 0},
/* 6 */ {"|2 2 5/2", "Pentagrammic Antiprism",
"Pentagrammic Deltohedron",
"Dihedral (D[2/5])",
"",
"",
0, 0},
/* (5/3 2 2) (D3/5) */
/* 8 */ {"|2 2 5/3", "Pentagrammic Crossed Antiprism",
"Pentagrammic Concave Deltohedron",
"Dihedral (D[3/5])",
"",
"",
0, 0},
/****************************************************************************
* Tetrahedral
***************************************************************************/
/* (2 3 3) (T1) */
/* 10 */ {"3|2 3", "Tetrahedron",
"Tetrahedron",
"Tetrahedral (T[1])",
"Platonic Solid",
"Platonic Solid",
15, 1},
/* 12 */ {"2 3|3", "Truncated Tetrahedron",
"Triakistetrahedron",
"Tetrahedral (T[1])",
"Archimedean Solid",
"Catalan Solid",
16, 6},
/* (3/2 3 3) (T2) */
/* 14 */ {"3/2 3|3", "Octahemioctahedron",
"Octahemioctacron",
"Tetrahedral (T[2])",
"",
"",
37, 68},
/* (3/2 2 3) (T3) */
/* 16 */ {"3/2 3|2", "Tetrahemihexahedron",
"Tetrahemihexacron",
"Tetrahedral (T[3])",
"",
"",
36, 67},
/****************************************************************************
* Octahedral
***************************************************************************/
/* (2 3 4) (O1) */
/* 18 */ {"4|2 3", "Octahedron",
"Cube",
"Octahedral (O[1])",
"Platonic Solid",
"Platonic Solid",
17, 2},
/* 20 */ {"3|2 4", "Cube",
"Octahedron",
"Octahedral (O[1])",
"Platonic Solid",
"Platonic Solid",
18, 3},
/* 22 */ {"2|3 4", "Cuboctahedron",
"Rhombic Dodecahedron",
"Octahedral (O[1])",
"Archimedean Solid",
"Catalan Solid",
19, 11},
/* 24 */ {"2 4|3", "Truncated Octahedron",
"Tetrakishexahedron",
"Octahedral (O[1])",
"Archimedean Solid",
"Catalan Solid",
20, 7},
/* 26 */ {"2 3|4", "Truncated Cube",
"Triakisoctahedron",
"Octahedral (O[1])",
"Archimedean Solid",
"Catalan Solid",
21, 8},
/* 28 */ {"3 4|2", "Rhombicuboctahedron",
"Deltoidal Icositetrahedron",
"Octahedral (O[1])",
"Archimedean Solid",
"Catalan Solid",
22, 13},
/* 30 */ {"2 3 4|", "Truncated Cuboctahedron",
"Disdyakisdodecahedron",
"Octahedral (O[1])",
"Archimedean Solid",
"Catalan Solid",
23, 15},
/* 32, 33, 66, and 67 are chiral, existing in both left and right handed
(enantiomeric) forms, so it would make sense to display both versions.
*/
/* 32 */ {"|2 3 4", "Snub Cube",
"Pentagonal Icositetrahedron",
"Octahedral (O[1]), Chiral",
"Archimedean Solid",
"Catalan Solid",
24, 17},
/* (3/2 4 4) (O2b) */
/* 34 */ {"3/2 4|4", "Small Cubicuboctahedron",
"Small Hexacronic Icositetrahedron",
"Octahedral (O[2b])",
"",
"",
38, 69},
/* (4/3 3 4) (O4) */
/* 36 */ {"3 4|4/3", "Great Cubicuboctahedron",
"Great Hexacronic Icositetrahedron",
"Octahedral (O[4])",
"",
"",
50, 77},
/* 38 */ {"4/3 4|3", "Cubohemioctahedron",
"Hexahemioctacron",
"Octahedral (O[4])",
"",
"",
51, 78},
/* 40 */ {"4/3 3 4|", "Cubitruncated Cuboctahedron",
"Tetradyakishexahedron",
"Octahedral (O[4])",
"",
"",
52, 79},
/* (3/2 2 4) (O5) */
/* 42 */ {"3/2 4|2", "Great Rhombicuboctahedron",
"Great Deltoidal Icositetrahedron",
"Octahedral (O[5])",
"",
"",
59, 85},
/* 44 */ {"3/2 2 4|", "Small Rhombihexahedron",
"Small Rhombihexacron",
"Octahedral (O[5])",
"",
"",
60, 86},
/* (4/3 2 3) (O7) */
/* 46 */ {"2 3|4/3", "Stellated Truncated Hexahedron",
"Great Triakisoctahedron",
"Octahedral (O[7])",
"",
"",
66, 92},
/* 48 */ {"4/3 2 3|", "Great Truncated Cuboctahedron",
"Great Disdyakisdodecahedron",
"Octahedral (O[7])",
"",
"",
67, 93},
/* (4/3 3/2 2) (O11) */
/* 50 */ {"4/3 3/2 2|", "Great Rhombihexahedron",
"Great Rhombihexacron",
"Octahedral (O[11])",
"",
"",
82, 103},
/****************************************************************************
* Icosahedral
***************************************************************************/
/* (2 3 5) (I1) */
/* 52 */ {"5|2 3", "Icosahedron",
"Dodecahedron",
"Icosahedral (I[1])",
"Platonic Solid",
"Platonic Solid",
25, 4},
/* 54 */ {"3|2 5", "Dodecahedron",
"Icosahedron",
"Icosahedral (I[1])",
"Platonic Solid",
"Platonic Solid",
26, 5},
/* 56 */ {"2|3 5", "Icosidodecahedron",
"Rhombic Triacontahedron",
"Icosahedral (I[1])",
"Archimedean Solid",
"Catalan Solid",
28, 12},
/* 58 */ {"2 5|3", "Truncated Icosahedron",
"Pentakisdodecahedron",
"Icosahedral (I[1])",
"Archimedean Solid",
"Catalan Solid",
27, 9},
/* 60 */ {"2 3|5", "Truncated Dodecahedron",
"Triakisicosahedron",
"Icosahedral (I[1])",
"Archimedean Solid",
"Catalan Solid",
29, 10},
/* 62 */ {"3 5|2", "Rhombicosidodecahedron",
"Deltoidal Hexecontahedron",
"Icosahedral (I[1])",
"Archimedean Solid",
"Catalan Solid",
30, 14},
/* 64 */ {"2 3 5|", "Truncated Icosidodecahedron",
"Disdyakistriacontahedron",
"Icosahedral (I[1])",
"Archimedean Solid",
"Catalan Solid",
31, 16},
/* 32, 33, 66, and 67 are chiral, existing in both left and right handed
(enantiomeric) forms, so it would make sense to display both versions.
*/
/* 66 */ {"|2 3 5", "Snub Dodecahedron",
"Pentagonal Hexecontahedron",
"Icosahedral (I[1]), Chiral",
"Archimedean Solid",
"Catalan Solid",
32, 18},
/* (5/2 3 3) (I2a) */
/* 68 */ {"3|5/2 3", "Small Ditrigonal Icosidodecahedron",
"Small Triambic Icosahedron",
"Icosahedral (I[2a])",
"",
"",
39, 70},
/* 70 */ {"5/2 3|3", "Small Icosicosidodecahedron",
"Small Icosacronic Hexecontahedron",
"Icosahedral (I[2a])",
"",
"",
40, 71},
/* 72 */ {"|5/2 3 3", "Small Snub Icosicosidodecahedron",
"Small Hexagonal Hexecontahedron",
"Icosahedral (I[2a])",
"",
"",
41, 110},
/* (3/2 5 5) (I2b) */
/* 74 */ {"3/2 5|5", "Small Dodecicosidodecahedron",
"Small Dodecacronic Hexecontahedron",
"Icosahedral (I[2b])",
"",
"",
42, 72},
/* (2 5/2 5) (I3) */
/* 76 */ {"5|2 5/2", "Small Stellated Dodecahedron",
"Great Dodecahedron",
"Icosahedral (I[3])",
"Kepler-Poinsot Solid",
"Kepler-Poinsot Solid",
43, 20},
/* 78 */ {"5/2|2 5", "Great Dodecahedron",
"Small Stellated Dodecahedron",
"Icosahedral (I[3])",
"Kepler-Poinsot Solid",
"Kepler-Poinsot Solid",
44, 21},
/* 80 */ {"2|5/2 5", "Great Dodecadodecahedron",
"Medial Rhombic Triacontahedron",
"Icosahedral (I[3])",
"",
"",
45, 73},
/* 82 */ {"2 5/2|5", "Truncated Great Dodecahedron",
"Small Stellapentakisdodecahedron",
"Icosahedral (I[3])",
"",
"",
47, 75},
/* 84 */ {"5/2 5|2", "Rhombidodecadodecahedron",
"Medial Deltoidal Hexecontahedron",
"Icosahedral (I[3])",
"",
"",
48, 76},
/* 86 */ {"2 5/2 5|", "Small Rhombidodecahedron",
"Small Rhombidodecacron",
"Icosahedral (I[3])",
"",
"",
46, 74},
/* 88 */ {"|2 5/2 5", "Snub Dodecadodecahedron",
"Medial Pentagonal Hexecontahedron",
"Icosahedral (I[3])",
"",
"",
49, 111},
/* (5/3 3 5) (I4) */
/* 90 */ {"3|5/3 5", "Ditrigonal Dodecadodecahedron",
"Medial Triambic Icosahedron",
"Icosahedral (I[4])",
"",
"",
53, 80},
/* 92 */ {"3 5|5/3", "Great Ditrigonal Dodecicosidodecahedron",
"Great Ditrigonal Dodecacronic Hexecontahedron",
"Icosahedral (I[4])",
"",
"",
54, 81},
/* 94 */ {"5/3 3|5", "Small Ditrigonal Dodecicosidodecahedron",
"Small Ditrigonal Dodecacronic Hexecontahedron",
"Icosahedral (I[4])",
"",
"",
55, 82},
/* 96 */ {"5/3 5|3", "Icosidodecadodecahedron",
"Medial Icosacronic Hexecontahedron",
"Icosahedral (I[4])",
"",
"",
56, 83},
/* 98 */ {"5/3 3 5|", "Icositruncated Dodecadodecahedron",
"Tridyakisicosahedron",
"Icosahedral (I[4])",
"",
"",
57, 84},
/* 100 */ {"|5/3 3 5", "Snub Icosidodecadodecahedron",
"Medial Hexagonal Hexecontahedron",
"Icosahedral (I[4])",
"",
"",
58, 112},
/* (3/2 3 5) (I6b) */
/* 102 */ {"3/2|3 5", "Great Ditrigonal Icosidodecahedron",
"Great Triambic Icosahedron",
"Icosahedral (I[6b])",
"",
"",
61, 87},
/* 104 */ {"3/2 5|3", "Great Icosicosidodecahedron",
"Great Icosacronic Hexecontahedron",
"Icosahedral (I[6b])",
"",
"",
62, 88},
/* 106 */ {"3/2 3|5", "Small Icosihemidodecahedron",
"Small Icosihemidodecacron",
"Icosahedral (I[6b])",
"",
"",
63, 89},
/* 108 */ {"3/2 3 5|", "Small Dodecicosahedron",
"Small Dodecicosacron",
"Icosahedral (I[6b])",
"",
"",
64, 90},
/* (5/4 5 5) (I6c) */
/* 110 */ {"5/4 5|5", "Small Dodecahemidodecahedron",
"Small Dodecahemidodecacron",
"Icosahedral (I[6c])",
"",
"",
65, 91},
/* (2 5/2 3) (I7) */
/* 112 */ {"3|2 5/2", "Great Stellated Dodecahedron",
"Great Icosahedron",
"Icosahedral (I[7])",
"Kepler-Poinsot Solid",
"Kepler-Poinsot Solid",
68, 22},
/* 114 */ {"5/2|2 3", "Great Icosahedron",
"Great Stellated Dodecahedron",
"Icosahedral (I[7])",
"Kepler-Poinsot Solid",
"Kepler-Poinsot Solid",
69, 41},
/* 116 */ {"2|5/2 3", "Great Icosidodecahedron",
"Great Rhombic Triacontahedron",
"Icosahedral (I[7])",
"",
"",
70, 94},
/* 118 */ {"2 5/2|3", "Great Truncated Icosahedron",
"Great Stellapentakisdodecahedron",
"Icosahedral (I[7])",
"",
"",
71, 95},
/* 120 */ {"2 5/2 3|", "Rhombicosahedron",
"Rhombicosacron",
"Icosahedral (I[7])",
"",
"",
72, 96},
/* 122 */ {"|2 5/2 3", "Great Snub Icosidodecahedron",
"Great Pentagonal Hexecontahedron",
"Icosahedral (I[7])",
"",
"",
73, 113},
/* (5/3 2 5) (I9) */
/* 124 */ {"2 5|5/3", "Small Stellated Truncated Dodecahedron",
"Great Pentakisdodecahedron",
"Icosahedral (I[9])",
"",
"",
74, 97},
/* 126 */ {"5/3 2 5|", "Truncated Dodecadodecahedron",
"Medial Disdyakistriacontahedron",
"Icosahedral (I[9])",
"",
"",
75, 98},
/* 128 */ {"|5/3 2 5", "Inverted Snub Dodecadodecahedron",
"Medial Inverted Pentagonal Hexecontahedron",
"Icosahedral (I[9])",
"",
"",
76, 114},
/* (5/3 5/2 3) (I10a) */
/* 130 */ {"5/2 3|5/3", "Great Dodecicosidodecahedron",
"Great Dodecacronic Hexecontahedron",
"Icosahedral (I[10a])",
"",
"",
77, 99},
/* 132 */ {"5/3 5/2|3", "Small Dodecahemicosahedron",
"Small Dodecahemicosacron",
"Icosahedral (I[10a])",
"",
"",
78, 100},
/* 134 */ {"5/3 5/2 3|", "Great Dodecicosahedron",
"Great Dodecicosacron",
"Icosahedral (I[10a])",
"",
"",
79, 101},
/* 136 */ {"|5/3 5/2 3", "Great Snub Dodecicosidodecahedron",
"Great Hexagonal Hexecontahedron",
"Icosahedral (I[10a])",
"",
"",
80, 115},
/* (5/4 3 5) (I10b) */
/* 138 */ {"5/4 5|3", "Great Dodecahemicosahedron",
"Great Dodecahemicosacron",
"Icosahedral (I[10b])",
"",
"",
81, 102},
/* (5/3 2 3) (I13) */
/* 140 */ {"2 3|5/3", "Great Stellated Truncated Dodecahedron",
"Great Triakisicosahedron",
"Icosahedral (I[13])",
"",
"",
83, 104},
/* 142 */ {"5/3 3|2", "Great Rhombicosidodecahedron",
"Great Deltoidal Hexecontahedron",
"Icosahedral (I[13])",
"",
"",
84, 105},
/* 144 */ {"5/3 2 3|", "Great Truncated Icosidodecahedron",
"Great Disdyakistriacontahedron",
"Icosahedral (I[13])",
"",
"",
87, 108},
/* 146 */ {"|5/3 2 3", "Great Inverted Snub Icosidodecahedron",
"Great Inverted Pentagonal Hexecontahedron",
"Icosahedral (I[13])",
"",
"",
88, 116},
/* (5/3 5/3 5/2) (I18a) */
/* 148 */ {"5/3 5/2|5/3", "Great Dodecahemidodecahedron",
"Great Dodecahemidodecacron",
"Icosahedral (I[18a])",
"",
"",
86, 107},
/* (3/2 5/3 3) (I18b) */
/* 150 */ {"3/2 3|5/3", "Great Icosihemidodecahedron",
"Great Icosihemidodecacron",
"Icosahedral (I[18b])",
"",
"",
85, 106},
/* (3/2 3/2 5/3) (I22) */
/* 152 */ {"|3/2 3/2 5/2","Small Retrosnub Icosicosidodecahedron",
"Small Hexagrammic Hexecontahedron",
"Icosahedral (I[22])",
"",
"",
91, 118},
/* (3/2 5/3 2) (I23) */
/* 154 */ {"3/2 5/3 2|", "Great Rhombidodecahedron",
"Great Rhombidodecacron",
"Icosahedral (I[23])",
"",
"",
89, 109},
/* 156 */ {"|3/2 5/3 2", "Great Retrosnub Icosidodecahedron",
"Great Pentagrammic Hexecontahedron",
"Icosahedral (I[23])",
"",
"",
90, 117},
/****************************************************************************
* Last But Not Least
***************************************************************************/
/* 158 */ {"3/2 5/3 3 5/2", "Great Dirhombicosidodecahedron",
"Great Dirhombicosidodecacron",
"Non-Wythoffian",
"",
"",
92, 119}
};
static int last_uniform = sizeof (uniform) / sizeof (uniform[0]);
static int unpacksym(const char *sym, Polyhedron *P);
static int moebius(Polyhedron *P);
static int decompose(Polyhedron *P);
static int guessname(Polyhedron *P);
static int newton(Polyhedron *P, int need_approx);
static int exceptions(Polyhedron *P);
static int count(Polyhedron *P);
static int configuration(Polyhedron *P);
static int vertices(Polyhedron *P);
static int faces(Polyhedron *P);
static int edgelist(Polyhedron *P);
static Polyhedron *
kaleido(const char *sym,
int need_coordinates, int need_edgelist, int need_approx,
int just_list)
{
Polyhedron *P;
/*
* Allocate a Polyhedron structure P.
*/
if (!(P = polyalloc()))
return 0;
/*
* Unpack input symbol into P.
*/
if (!unpacksym(sym, P))
return 0;
/*
* Find Mebius triangle, its density and Euler characteristic.
*/
if (!moebius(P))
return 0;
/*
* Decompose Schwarz triangle.
*/
if (!decompose(P))
return 0;
/*
* Find the names of the polyhedron and its dual.
*/
if (!guessname(P))
return 0;
if (just_list)
return P;
/*
* Solve Fundamental triangles, optionally printing approximations.
*/
if (!newton(P,need_approx))
return 0;
/*
* Deal with exceptional polyhedra.
*/
if (!exceptions(P))
return 0;
/*
* Count edges and faces, update density and characteristic if needed.
*/
if (!count(P))
return 0;
/*
* Generate printable vertex configuration.
*/
if (!configuration(P))
return 0;
/*
* Compute coordinates.
*/
if (!need_coordinates && !need_edgelist)
return P;
if (!vertices(P))
return 0;
if (!faces (P))
return 0;
/*
* Compute edgelist.
*/
if (!need_edgelist)
return P;
if (!edgelist(P))
return 0;
return P;
}
/*
* Allocate a blank Polyhedron structure and initialize some of its nonblank
* fields.
*
* Array and matrix field are allocated when needed.
*/
static Polyhedron *
polyalloc()
{
Polyhedron *P;
Calloc(P, 1, Polyhedron);
P->index = -1;
P->even = -1;
P->K = 2;
return P;
}
/*
* Free the struture allocated by polyalloc(), as well as all the array and
* matrix fields.
*/
static void
polyfree(Polyhedron *P)
{
Free(P->Fi);
Free(P->n);
Free(P->m);
Free(P->gamma);
Free(P->rot);
Free(P->snub);
Free(P->firstrot);
Free(P->anti);
Free(P->ftype);
Free(P->polyform);
Free(P->config);
if (P->index < 0) {
Free(P->name);
Free(P->dual_name);
}
Free(P->v);
Free(P->f);
Matfree(P->e, 2);
Matfree(P->dual_e, 2);
Matfree(P->incid, P->M);
Matfree(P->adj, P->M);
free(P);
}
static void *
matalloc(int rows, int row_size)
{
void **mat;
int i = 0;
if (!(mat = malloc(rows * sizeof (void *))))
return 0;
while ((mat[i] = malloc(row_size)) && ++i < rows)
;
if (i == rows)
return (void *)mat;
while (--i >= 0)
free(mat[i]);
free(mat);
return 0;
}
static void
matfree(void *mat, int rows)
{
while (--rows >= 0)
free(((void **)mat)[rows]);
free(mat);
}
/*
* compute the mathematical modulus function.
*/
static int
mod (int i, int j)
{
return (i%=j)>=0?i:j<0?i-j:i+j;
}
/*
* Find the numerator and the denominator using the Euclidean algorithm.
*/
static void
frac(double x)
{
static const Fraction zero = {0,1}, inf = {1,0};
Fraction r0, r;
long f;
double s = x;
r = zero;
frax = inf;
for (;;) {
if (fabs(s) > (double) MAXLONG)
return;
f = (long) floor (s);
r0 = r;
r = frax;
frax.n = frax.n * f + r0.n;
frax.d = frax.d * f + r0.d;
if (x == (double)frax.n/(double)frax.d)
return;
s = 1 / (s - f);
}
}
/*
* Unpack input symbol: Wythoff symbol or an index to uniform[]. The symbol is
* a # followed by a number, or a three fractions and a bar in some order. We
* allow no bars only if it result from the input symbol #80.
*/
static int
unpacksym(const char *sym, Polyhedron *P)
{
int i = 0, n, d, bars = 0;
char c;
while ((c = *sym++) && isspace(c))
;
if (!c) Err("no data");
if (c == '#') {
while ((c = *sym++) && isspace(c))
;
if (!c)
Err("no digit after #");
if (!isdigit(c))
Err("not a digit");
n = c - '0';
while ((c = *sym++) && isdigit(c))
n = n * 10 + c - '0';
if (!n)
Err("zero index");
if (n > last_uniform)
Err("index too big");
sym--;
while ((c = *sym++) && isspace(c))
;
if (c)
Err("data exceeded");
sym = uniform[P->index = n - 1].Wythoff;
} else
sym--;
for (;;) {
while ((c = *sym++) && isspace(c))
;
if (!c) {
if (i == 4 && (bars || P->index == last_uniform - 1))
return 1;
if (!bars)
Err("no bars");
Err("not enough fractions");
}
if (i == 4)
Err("data exceeded");
if (c == '|'){
if (++bars > 1)
Err("too many bars");
P->p[i++] = 0;
continue;
}
if (!isdigit(c))
Err("not a digit");
n = c - '0';
while ((c = *sym++) && isdigit(c))
n = n * 10 + c - '0';
if (c && isspace (c))
while ((c = *sym++) && isspace(c))
;
if (c != '/') {
sym--;
if ((P->p[i++] = n) <= 1)
Err("fraction<=1");
continue;
}
while ((c = *sym++) && isspace(c))
;
if (!c || !isdigit(c))
return 0;
d = c - '0';
while ((c = *sym++) && isdigit(c))
d = d * 10 + c - '0';
if (!d)
Err("zero denominator");
sym--;
if ((P->p[i++] = (double) n / d) <= 1)
Err("fraction<=1");
}
}
/*
* Using Wythoff symbol (p|qr, pq|r, pqr| or |pqr), find the Moebius triangle
* (2 3 K) (or (2 2 n)) of the Schwarz triangle (pqr), the order g of its
* symmetry group, its Euler characteristic chi, and its covering density D.
* g is the number of copies of (2 3 K) covering the sphere, i.e.,
*
* g * pi * (1/2 + 1/3 + 1/K - 1) = 4 * pi
*
* D is the number of times g copies of (pqr) cover the sphere, i.e.
*
* D * 4 * pi = g * pi * (1/p + 1/q + 1/r - 1)
*
* chi is V - E + F, where F = g is the number of triangles, E = 3*g/2 is the
* number of triangle edges, and V = Vp+ Vq+ Vr, with Vp = g/(2*np) being the
* number of vertices with angle pi/p (np is the numerator of p).
*/
static int
moebius(Polyhedron *P)
{
int twos = 0, j, len = 1;
/*
* Arrange Wythoff symbol in a presentable form. In the same time check the
* restrictions on the three fractions: They all have to be greater then one,
* and the numerators 4 or 5 cannot occur together. We count the ocurrences
* of 2 in `two', and save the largest numerator in `P->K', since they
* reflect on the symmetry group.
*/
P->K = 2;
if (P->index == last_uniform - 1) {
Malloc(P->polyform, ++len, char);
strcpy(P->polyform, "|");
} else
Calloc(P->polyform, len, char);
for (j = 0; j < 4; j++) {
if (P->p[j]) {
char *s;
Sprintfrac(s, P->p[j]);
if (j && P->p[j-1]) {
Realloc(P->polyform, len += strlen (s) + 1, char);
strcat(P->polyform, " ");
} else
Realloc (P->polyform, len += strlen (s), char);
strcat(P->polyform, s);
free(s);
if (P->p[j] != 2) {
int k;
if ((k = numerator (P->p[j])) > P->K) {
if (P->K == 4)
break;
P->K = k;
} else if (k < P->K && k == 4)
break;
} else
twos++;
} else {
Realloc(P->polyform, ++len, char);
strcat(P->polyform, "|");
}
}
/*
* Find the symmetry group P->K (where 2, 3, 4, 5 represent the dihedral,
* tetrahedral, octahedral and icosahedral groups, respectively), and its
* order P->g.
*/
if (twos >= 2) {/* dihedral */
P->g = 4 * P->K;
P->K = 2;
} else {
if (P->K > 5)
Err("numerator too large");
P->g = 24 * P->K / (6 - P->K);
}
/*
* Compute the nominal density P->D and Euler characteristic P->chi.
* In few exceptional cases, these values will be modified later.
*/
if (P->index != last_uniform - 1) {
int i;
P->D = P->chi = - P->g;
for (j = 0; j < 4; j++) if (P->p[j]) {
P->chi += i = P->g / numerator(P->p[j]);
P->D += i * denominator(P->p[j]);
}
P->chi /= 2;
P->D /= 4;
if (P->D <= 0)
Err("nonpositive density");
}
return 1;
}
/*
* Decompose Schwarz triangle into N right triangles and compute the vertex
* count V and the vertex valency M. V is computed from the number g of
* Schwarz triangles in the cover, divided by the number of triangles which
* share a vertex. It is halved for one-sided polyhedra, because the
* kaleidoscopic construction really produces a double orientable covering of
* such polyhedra. All q' q|r are of the "hemi" type, i.e. have equatorial {2r}
* faces, and therefore are (except 3/2 3|3 and the dihedra 2 2|r) one-sided. A
* well known example is 3/2 3|4, the "one-sided heptahedron". Also, all p q r|
* with one even denominator have a crossed parallelogram as a vertex figure,
* and thus are one-sided as well.
*/
static int
decompose(Polyhedron *P)
{
int j, J, *s, *t;
if (!P->p[1]) { /* p|q r */
P->N = 2;
P->M = 2 * numerator(P->p[0]);
P->V = P->g / P->M;
Malloc(P->n, P->N, double);
Malloc(P->m, P->N, double);
Malloc(P->rot, P->M, int);
s = P->rot;
for (j = 0; j < 2; j++) {
P->n[j] = P->p[j+2];
P->m[j] = P->p[0];
}
for (j = P->M / 2; j--;) {
*s++ = 0;
*s++ = 1;
}
} else if (!P->p[2]) { /* p q|r */
P->N = 3;
P->M = 4;
P->V = P->g / 2;
Malloc(P->n, P->N, double);
Malloc(P->m, P->N, double);
Malloc(P->rot, P->M, int);
s = P->rot;
P->n[0] = 2 * P->p[3];
P->m[0] = 2;
for (j = 1; j < 3; j++) {
P->n[j] = P->p[j-1];
P->m[j] = 1;
*s++ = 0;
*s++ = j;
}
if (fabs(P->p[0] - compl (P->p[1])) < DBL_EPSILON) {/* p = q' */
/* P->p[0]==compl(P->p[1]) should work. However, MSDOS
* yeilds a 7e-17 difference! Reported by Jim Buddenhagen
* <jb1556@daditz.sbc.com> */
P->hemi = 1;
P->D = 0;
if (P->p[0] != 2 && !(P->p[3] == 3 && (P->p[0] == 3 ||
P->p[1] == 3))) {
P->onesided = 1;
P->V /= 2;
P->chi /= 2;
}
}
} else if (!P->p[3]) { /* p q r| */
P->M = P->N = 3;
P->V = P->g;
Malloc(P->n, P->N, double);
Malloc(P->m, P->N, double);
Malloc(P->rot, P->M, int);
s = P->rot;
for (j = 0; j < 3; j++) {
if (!(denominator(P->p[j]) % 2)) {
/* what happens if there is more then one even denominator? */
if (P->p[(j+1)%3] != P->p[(j+2)%3]) { /* needs postprocessing */
P->even = j;/* memorize the removed face */
P->chi -= P->g / numerator(P->p[j]) / 2;
P->onesided = 1;
P->D = 0;
} else {/* for p = q we get a double 2 2r|p */
/* noted by Roman Maeder <maeder@inf.ethz.ch> for 4 4 3/2| */
/* Euler characteristic is still wrong */
P->D /= 2;
}
P->V /= 2;
}
P->n[j] = 2 * P->p[j];
P->m[j] = 1;
*s++ = j;
}
} else { /* |p q r - snub polyhedron */
P->N = 4;
P->M = 6;
P->V = P->g / 2;/* Only "white" triangles carry a vertex */
Malloc(P->n, P->N, double);
Malloc(P->m, P->N, double);
Malloc(P->rot, P->M, int);
Malloc(P->snub, P->M, int);
s = P->rot;
t = P->snub;
P->m[0] = P->n[0] = 3;
for (j = 1; j < 4; j++) {
P->n[j] = P->p[j];
P->m[j] = 1;
*s++ = 0;
*s++ = j;
*t++ = 1;
*t++ = 0;
}
}
/*
* Sort the fundamental triangles (using bubble sort) according to decreasing
* n[i], while pushing the trivial triangles (n[i] = 2) to the end.
*/
J = P->N - 1;
while (J) {
int last;
last = J;
J = 0;
for (j = 0; j < last; j++) {
if ((P->n[j] < P->n[j+1] || P->n[j] == 2) && P->n[j+1] != 2) {
int i;
double temp;
temp = P->n[j];
P->n[j] = P->n[j+1];
P->n[j+1] = temp;
temp = P->m[j];
P->m[j] = P->m[j+1];
P->m[j+1] = temp;
for (i = 0; i < P->M; i++) {
if (P->rot[i] == j)
P->rot[i] = j+1;
else if (P->rot[i] == j+1)
P->rot[i] = j;
}
if (P->even != -1) {
if (P->even == j)
P->even = j+1;
else if (P->even == j+1)
P->even = j;
}
J = j;
}
}
}
/*
* Get rid of repeated triangles.
*/
for (J = 0; J < P->N && P->n[J] != 2;J++) {
int k, i;
for (j = J+1; j < P->N && P->n[j]==P->n[J]; j++)
P->m[J] += P->m[j];
k = j - J - 1;
if (k) {
for (i = j; i < P->N; i++) {
P->n[i - k] = P->n[i];
P->m[i - k] = P->m[i];
}
P->N -= k;
for (i = 0; i < P->M; i++) {
if (P->rot[i] >= j)
P->rot[i] -= k;
else if (P->rot[i] > J)
P->rot[i] = J;
}
if (P->even >= j)
P->even -= k;
}
}
/*
* Get rid of trivial triangles.
*/
if (!J)
J = 1; /* hosohedron */
if (J < P->N) {
int i;
P->N = J;
for (i = 0; i < P->M; i++) {
if (P->rot[i] >= P->N) {
for (j = i + 1; j < P->M; j++) {
P->rot[j-1] = P->rot[j];
if (P->snub)
P->snub[j-1] = P->snub[j];
}
P->M--;
}
}
}
/*
* Truncate arrays
*/
Realloc(P->n, P->N, double);
Realloc(P->m, P->N, double);
Realloc(P->rot, P->M, int);
if (P->snub)
Realloc(P->snub, P->M, int);
return 1;
}
static int dihedral(Polyhedron *P, const char *name, const char *dual_name);
/*
* Get the polyhedron name, using standard list or guesswork. Ideally, we
* should try to locate the Wythoff symbol in the standard list (unless, of
* course, it is dihedral), after doing few normalizations, such as sorting
* angles and splitting isoceles triangles.
*/
static int
guessname(Polyhedron *P)
{
if (P->index != -1) {/* tabulated */
P->name = uniform[P->index].name;
P->dual_name = uniform[P->index].dual;
P->group = uniform[P->index].group;
P->class = uniform[P->index].class;
P->dual_class = uniform[P->index].dual_class;
return 1;
} else if (P->K == 2) {/* dihedral nontabulated */
if (!P->p[0]) {
if (P->N == 1) {
Malloc(P->name, sizeof ("Octahedron"), char);
Malloc(P->dual_name, sizeof ("Cube"), char);
strcpy(P->name, "Octahedron");
strcpy(P->dual_name, "Cube");
return 1;
}
P->gon = P->n[0] == 3 ? P->n[1] : P->n[0];
if (P->gon >= 2)
return dihedral(P, "Antiprism", "Deltohedron");
else
return dihedral(P, "Crossed Antiprism", "Concave Deltohedron");
} else if (!P->p[3] ||
(!P->p[2] &&
P->p[3] == 2)) {
if (P->N == 1) {
Malloc(P->name, sizeof("Cube"), char);
Malloc(P->dual_name, sizeof("Octahedron"), char);
strcpy(P->name, "Cube");
strcpy(P->dual_name, "Octahedron");
return 1;
}
P->gon = P->n[0] == 4 ? P->n[1] : P->n[0];
return dihedral(P, "Prism", "Dipyramid");
} else if (!P->p[1] && P->p[0] != 2) {
P->gon = P->m[0];
return dihedral(P, "Hosohedron", "Dihedron");
} else {
P->gon = P->n[0];
return dihedral(P, "Dihedron", "Hosohedron");
}
} else {/* other nontabulated */
static const char *pre[] = {"Tetr", "Oct", "Icos"};
Malloc(P->name, 50, char);
Malloc(P->dual_name, 50, char);
sprintf(P->name, "%sahedral ", pre[P->K - 3]);
if (P->onesided)
strcat (P->name, "One-Sided ");
else if (P->D == 1)
strcat(P->name, "Convex ");
else
strcat(P->name, "Nonconvex ");
strcpy(P->dual_name, P->name);
strcat(P->name, "Isogonal Polyhedron");
strcat(P->dual_name, "Isohedral Polyhedron");
Realloc(P->name, strlen (P->name) + 1, char);
Realloc(P->dual_name, strlen (P->dual_name) + 1, char);
return 1;
}
}
static int
dihedral(Polyhedron *P, const char *name, const char *dual_name)
{
char *s;
int i;
Sprintfrac(s, P->gon < 2 ? compl (P->gon) : P->gon);
i = strlen(s) + sizeof ("-gonal ");
Malloc(P->name, i + strlen (name), char);
Malloc(P->dual_name, i + strlen (dual_name), char);
sprintf(P->name, "%s-gonal %s", s, name);
sprintf(P->dual_name, "%s-gonal %s", s, dual_name);
free(s);
return 1;
}
/*
* Solve the fundamental right spherical triangles.
* If need_approx is set, print iterations on standard error.
*/
static int
newton(Polyhedron *P, int need_approx)
{
/*
* First, we find initial approximations.
*/
int j;
double cosa;
Malloc(P->gamma, P->N, double);
if (P->N == 1) {
P->gamma[0] = M_PI / P->m[0];
return 1;
}
for (j = 0; j < P->N; j++)
P->gamma[j] = M_PI / 2 - M_PI / P->n[j];
errno = 0; /* may be non-zero from some reason */
/*
* Next, iteratively find closer approximations for gamma[0] and compute
* other gamma[j]'s from Napier's equations.
*/
if (need_approx)
fprintf(stderr, "Solving %s\n", P->polyform);
for (;;) {
double delta = M_PI, sigma = 0;
for (j = 0; j < P->N; j++) {
if (need_approx)
fprintf(stderr, "%-20.15f", P->gamma[j]);
delta -= P->m[j] * P->gamma[j];
}
if (need_approx)
printf("(%g)\n", delta);
if (fabs(delta) < 11 * DBL_EPSILON)
return 1;
/* On a RS/6000, fabs(delta)/DBL_EPSILON may occilate between 8 and
* 10. Reported by David W. Sanderson <dws@ssec.wisc.edu> */
for (j = 0; j < P->N; j++)
sigma += P->m[j] * tan(P->gamma[j]);
P->gamma[0] += delta * tan(P->gamma[0]) / sigma;
if (P->gamma[0] < 0 || P->gamma[0] > M_PI)
Err("gamma out of bounds");
cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
for (j = 1; j < P->N; j++)
P->gamma[j] = asin(cos(M_PI / P->n[j]) / cosa);
if (errno)
Err(strerror(errno));
}
}
/*
* Postprocess pqr| where r has an even denominator (cf. Coxeter &al. Sec.9).
* Remove the {2r} and add a retrograde {2p} and retrograde {2q}.
*/
static int
exceptions(Polyhedron *P)
{
int j;
if (P->even != -1) {
P->M = P->N = 4;
Realloc(P->n, P->N, double);
Realloc(P->m, P->N, double);
Realloc(P->gamma, P->N, double);
Realloc(P->rot, P->M, int);
for (j = P->even + 1; j < 3; j++) {
P->n[j-1] = P->n[j];
P->gamma[j-1] = P->gamma[j];
}
P->n[2] = compl(P->n[1]);
P->gamma[2] = - P->gamma[1];
P->n[3] = compl(P->n[0]);
P->m[3] = 1;
P->gamma[3] = - P->gamma[0];
P->rot[0] = 0;
P->rot[1] = 1;
P->rot[2] = 3;
P->rot[3] = 2;
}
/*
* Postprocess the last polyhedron |3/2 5/3 3 5/2 by taking a |5/3 3 5/2,
* replacing the three snub triangles by four equatorial squares and adding
* the missing {3/2} (retrograde triangle, cf. Coxeter &al. Sec. 11).
*/
if (P->index == last_uniform - 1) {
P->N = 5;
P->M = 8;
Realloc(P->n, P->N, double);
Realloc(P->m, P->N, double);
Realloc(P->gamma, P->N, double);
Realloc(P->rot, P->M, int);
Realloc(P->snub, P->M, int);
P->hemi = 1;
P->D = 0;
for (j = 3; j; j--) {
P->m[j] = 1;
P->n[j] = P->n[j-1];
P->gamma[j] = P->gamma[j-1];
}
P->m[0] = P->n[0] = 4;
P->gamma[0] = M_PI / 2;
P->m[4] = 1;
P->n[4] = compl(P->n[1]);
P->gamma[4] = - P->gamma[1];
for (j = 1; j < 6; j += 2) P->rot[j]++;
P->rot[6] = 0;
P->rot[7] = 4;
P->snub[6] = 1;
P->snub[7] = 0;
}
return 1;
}
/*
* Compute edge and face counts, and update D and chi. Update D in the few
* cases the density of the polyhedron is meaningful but different than the
* density of the corresponding Schwarz triangle (cf. Coxeter &al., p. 418 and
* p. 425).
* In these cases, spherical faces of one type are concave (bigger than a
* hemisphere), and the actual density is the number of these faces less the
* computed density. Note that if j != 0, the assignment gamma[j] = asin(...)
* implies gamma[j] cannot be obtuse. Also, compute chi for the only
* non-Wythoffian polyhedron.
*/
static int
count(Polyhedron *P)
{
int j, temp;
Malloc(P->Fi, P->N, int);
for (j = 0; j < P->N; j++) {
P->E += temp = P->V * numerator(P->m[j]);
P->F += P->Fi[j] = temp / numerator(P->n[j]);
}
P->E /= 2;
if (P->D && P->gamma[0] > M_PI / 2)
P->D = P->Fi[0] - P->D;
if (P->index == last_uniform - 1)
P->chi = P->V - P->E + P->F;
return 1;
}
/*
* Generate a printable vertex configuration symbol.
*/
static int
configuration(Polyhedron *P)
{
int j, len = 2;
for (j = 0; j < P->M; j++) {
char *s;
Sprintfrac(s, P->n[P->rot[j]]);
len += strlen (s) + 2;
if (!j) {
Malloc(P->config, len, char);
/* strcpy(P->config, "(");*/
strcpy(P->config, "");
} else {
Realloc(P->config, len, char);
strcat(P->config, ", ");
}
strcat(P->config, s);
free(s);
}
/* strcat (P->config, ")");*/
if ((j = denominator (P->m[0])) != 1) {
char s[MAXDIGITS + 2];
sprintf(s, "/%d", j);
Realloc(P->config, len + strlen (s), char);
strcat(P->config, s);
}
return 1;
}
/*
* Compute polyhedron vertices and vertex adjecency lists.
* The vertices adjacent to v[i] are v[adj[0][i], v[adj[1][i], ...
* v[adj[M-1][i], ordered counterclockwise. The algorith is a BFS on the
* vertices, in such a way that the vetices adjacent to a givem vertex are
* obtained from its BFS parent by a cyclic sequence of rotations. firstrot[i]
* points to the first rotaion in the sequence when applied to v[i]. Note that
* for non-snub polyhedra, the rotations at a child are opposite in sense when
* compared to the rotations at the parent. Thus, we fill adj[*][i] from the
* end to signify clockwise rotations. The firstrot[] array is not needed for
* display thus it is freed after being used for face computations below.
*/
static int
vertices(Polyhedron *P)
{
int i, newV = 2;
double cosa;
Malloc(P->v, P->V, Vector);
Matalloc(P->adj, P->M, P->V, int);
Malloc(P->firstrot, P->V, int); /* temporary , put in Polyhedron
structure so that may be freed on
error */
cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
P->v[0].x = 0;
P->v[0].y = 0;
P->v[0].z = 1;
P->firstrot[0] = 0;
P->adj[0][0] = 1;
P->v[1].x = 2 * cosa * sqrt(1 - cosa * cosa);
P->v[1].y = 0;
P->v[1].z = 2 * cosa * cosa - 1;
if (!P->snub) {
P->firstrot[1] = 0;
P->adj[0][1] = -1;/* start the other side */
P->adj[P->M-1][1] = 0;
} else {
P->firstrot[1] = P->snub[P->M-1] ? 0 : P->M-1 ;
P->adj[0][1] = 0;
}
for (i = 0; i < newV; i++) {
int j, k;
int last, one, start, limit;
if (P->adj[0][i] == -1) {
one = -1; start = P->M-2; limit = -1;
} else {
one = 1; start = 1; limit = P->M;
}
k = P->firstrot[i];
for (j = start; j != limit; j += one) {
Vector temp;
int J;
temp = rotate (P->v[P->adj[j-one][i]], P->v[i],
one * 2 * P->gamma[P->rot[k]]);
for (J=0; J<newV && !same(P->v[J],temp,BIG_EPSILON); J++)
;/* noop */
P->adj[j][i] = J;
last = k;
if (++k == P->M)
k = 0;
if (J == newV) { /* new vertex */
if (newV == P->V) Err ("too many vertices");
P->v[newV++] = temp;
if (!P->snub) {
P->firstrot[J] = k;
if (one > 0) {
P->adj[0][J] = -1;
P->adj[P->M-1][J] = i;
} else {
P->adj[0][J] = i;
}
} else {
P->firstrot[J] = !P->snub[last] ? last :
!P->snub[k] ? (k+1)%P->M : k ;
P->adj[0][J] = i;
}
}
}
}
return 1;
}
/*
* Compute polyhedron faces (dual vertices) and incidence matrices.
* For orientable polyhedra, we can distinguish between the two faces meeting
* at a given directed edge and identify the face on the left and the face on
* the right, as seen from the outside. For one-sided polyhedra, the vertex
* figure is a papillon (in Coxeter &al. terminology, a crossed parallelogram)
* and the two faces meeting at an edge can be identified as the side face
* (n[1] or n[2]) and the diagonal face (n[0] or n[3]).
*/
static int
faces(Polyhedron *P)
{
int i, newF = 0;
Malloc (P->f, P->F, Vector);
Malloc (P->ftype, P->F, int);
Matalloc (P->incid, P->M, P->V, int);
P->minr = 1 / fabs (tan (M_PI / P->n[P->hemi]) * tan (P->gamma[P->hemi]));
for (i = P->M; --i>=0;) {
int j;
for (j = P->V; --j>=0;)
P->incid[i][j] = -1;
}
for (i = 0; i < P->V; i++) {
int j;
for (j = 0; j < P->M; j++) {
int i0, J;
int pap=0;/* papillon edge type */
if (P->incid[j][i] != -1)
continue;
P->incid[j][i] = newF;
if (newF == P->F)
Err("too many faces");
P->f[newF] = pole(P->minr, P->v[i], P->v[P->adj[j][i]],
P->v[P->adj[mod(j + 1, P->M)][i]]);
P->ftype[newF] = P->rot[mod(P->firstrot[i] + ((P->adj[0][i] <
P->adj[P->M - 1][i])
? j
: -j - 2),
P->M)];
if (P->onesided)
pap = (P->firstrot[i] + j) % 2;
i0 = i;
J = j;
for (;;) {
int k;
k = i0;
if ((i0 = P->adj[J][k]) == i) break;
for (J = 0; J < P->M && P->adj[J][i0] != k; J++)
;/* noop */
if (J == P->M)
Err("too many faces");
if (P->onesided && (J + P->firstrot[i0]) % 2 == pap) {
P->incid [J][i0] = newF;
if (++J >= P->M)
J = 0;
} else {
if (--J < 0)
J = P->M - 1;
P->incid [J][i0] = newF;
}
}
newF++;
}
}
Free(P->firstrot);
Free(P->rot);
Free(P->snub);
return 1;
}
/*
* Compute edge list and graph polyhedron and dual.
* If the polyhedron is of the "hemi" type, each edge has one finite vertex and
* one ideal vertex. We make sure the latter is always the out-vertex, so that
* the edge becomes a ray (half-line). Each ideal vertex is represented by a
* unit Vector, and the direction of the ray is either parallel or
* anti-parallel this Vector. We flag this in the array P->anti[E].
*/
static int
edgelist(Polyhedron *P)
{
int i, j, *s, *t, *u;
Matalloc(P->e, 2, P->E, int);
Matalloc(P->dual_e, 2, P->E, int);
s = P->e[0];
t = P->e[1];
for (i = 0; i < P->V; i++)
for (j = 0; j < P->M; j++)
if (i < P->adj[j][i]) {
*s++ = i;
*t++ = P->adj[j][i];
}
s = P->dual_e[0];
t = P->dual_e[1];
if (!P->hemi)
P->anti = 0;
else
Malloc(P->anti, P->E, int);
u = P->anti;
for (i = 0; i < P->V; i++)
for (j = 0; j < P->M; j++)
if (i < P->adj[j][i])
{
if (!u) {
*s++ = P->incid[mod(j-1,P->M)][i];
*t++ = P->incid[j][i];
} else {
if (P->ftype[P->incid[j][i]]) {
*s = P->incid[j][i];
*t = P->incid[mod(j-1,P->M)][i];
} else {
*s = P->incid[mod(j-1,P->M)][i];
*t = P->incid[j][i];
}
*u++ = dot(P->f[*s++], P->f[*t++]) > 0;
}
}
return 1;
}
static char *
sprintfrac(double x)
{
char *s;
frac (x);
if (!frax.d) {
Malloc(s, sizeof ("infinity"), char);
strcpy(s, "infinity");
} else if (frax.d == 1) {
char n[MAXDIGITS + 1];
sprintf(n, "%ld", frax.n);
Malloc(s, strlen (n) + 1, char);
strcpy(s, n);
} else {
char n[MAXDIGITS + 1], d[MAXDIGITS + 1];
sprintf(n, "%ld", frax.n);
sprintf(d, "%ld", frax.d);
Malloc(s, strlen (n) + strlen (d) + 2, char);
sprintf(s, "%s/%s", n, d);
}
return s;
}
static double
dot(Vector a, Vector b)
{
return a.x * b.x + a.y * b.y + a.z * b.z;
}
static Vector
scale(double k, Vector a)
{
a.x *= k;
a.y *= k;
a.z *= k;
return a;
}
static Vector
diff(Vector a, Vector b)
{
a.x -= b.x;
a.y -= b.y;
a.z -= b.z;
return a;
}
static Vector
cross(Vector a, Vector b)
{
Vector p;
p.x = a.y * b.z - a.z * b.y;
p.y = a.z * b.x - a.x * b.z;
p.z = a.x * b.y - a.y * b.x;
return p;
}
static Vector
sum(Vector a, Vector b)
{
a.x += b.x;
a.y += b.y;
a.z += b.z;
return a;
}
static Vector
sum3(Vector a, Vector b, Vector c)
{
a.x += b.x + c.x;
a.y += b.y + c.y;
a.z += b.z + c.z;
return a;
}
static Vector
rotate(Vector vertex, Vector axis, double angle)
{
Vector p;
p = scale(dot (axis, vertex), axis);
return sum3(p, scale(cos(angle), diff(vertex, p)),
scale(sin(angle), cross(axis, vertex)));
}
static Vector x, y, z;
/*
* rotate the standard frame
*/
static void
rotframe(double azimuth, double elevation, double angle)
{
static const Vector X = {1,0,0}, Y = {0,1,0}, Z = {0,0,1};
Vector axis;
axis = rotate(rotate (X, Y, elevation), Z, azimuth);
x = rotate(X, axis, angle);
y = rotate(Y, axis, angle);
z = rotate(Z, axis, angle);
}
/*
* rotate an array of n Vectors
*/
static void
rotarray(Vector *new, Vector *old, int n)
{
while (n--) {
*new++ = sum3(scale(old->x, x), scale(old->y, y), scale(old->z, z));
old++;
}
}
static int
same(Vector a, Vector b, double epsilon)
{
return fabs(a.x - b.x) < epsilon && fabs(a.y - b.y) < epsilon
&& fabs(a.z - b.z) < epsilon;
}
/*
* Compute the polar reciprocal of the plane containing a, b and c:
*
* If this plane does not contain the origin, return p such that
* dot(p,a) = dot(p,b) = dot(p,b) = r.
*
* Otherwise, return p such that
* dot(p,a) = dot(p,b) = dot(p,c) = 0
* and
* dot(p,p) = 1.
*/
static Vector
pole(double r, Vector a, Vector b, Vector c)
{
Vector p;
double k;
p = cross(diff(b, a), diff(c, a));
k = dot(p, a);
if (fabs(k) < 1e-6)
return scale(1 / sqrt(dot(p, p)), p);
else
return scale(r/ k , p);
}
/* output */
static void rotframe(double azimuth, double elevation, double angle);
static void rotarray(Vector *new, Vector *old, int n);
static int mod (int i, int j);
static void
push_point (polyhedron *p, Vector v)
{
p->points[p->npoints].x = v.x;
p->points[p->npoints].y = v.y;
p->points[p->npoints].z = v.z;
p->npoints++;
}
static void
push_face3 (polyhedron *p, int x, int y, int z)
{
p->faces[p->nfaces].npoints = 3;
Malloc (p->faces[p->nfaces].points, 3, int);
p->faces[p->nfaces].points[0] = x;
p->faces[p->nfaces].points[1] = y;
p->faces[p->nfaces].points[2] = z;
p->nfaces++;
}
static void
push_face4 (polyhedron *p, int x, int y, int z, int w)
{
p->faces[p->nfaces].npoints = 4;
Malloc (p->faces[p->nfaces].points, 4, int);
p->faces[p->nfaces].points[0] = x;
p->faces[p->nfaces].points[1] = y;
p->faces[p->nfaces].points[2] = z;
p->faces[p->nfaces].points[3] = w;
p->nfaces++;
}
static polyhedron *
construct_polyhedron (Polyhedron *P, Vector *v, int V, Vector *f, int F,
const char *name, const char *dual,
const char *class, const char *star,
double azimuth, double elevation, double freeze)
{
int i, j, k=0, l, ll, ii, *hit=0, facelets;
polyhedron *result;
Vector *temp;
Malloc (result, 1, polyhedron);
memset (result, 0, sizeof(*result));
/*
* Rotate polyhedron
*/
rotframe(azimuth, elevation, freeze);
Malloc(temp, V, Vector);
rotarray(temp, v, V);
v = temp;
Malloc(temp, F, Vector);
rotarray(temp, f, F);
f = temp;
result->number = P->index + 1;
result->name = strdup (name);
result->dual = strdup (dual);
result->wythoff = strdup (P->polyform);
result->config = strdup (P->config);
result->group = strdup (P->group);
result->class = strdup (class);
/*
* Vertex list
*/
Malloc (result->points, V + F * 13, point);
result->npoints = 0;
result->nedges = P->E;
result->logical_faces = F;
result->logical_vertices = V;
result->density = P->D;
result->chi = P->chi;
for (i = 0; i < V; i++)
push_point (result, v[i]);
/*
* Auxiliary vertices (needed because current VRML browsers cannot handle
* non-simple polygons, i.e., ploygons with self intersections): Each
* non-simple face is assigned an auxiliary vertex. By connecting it to the
* rest of the vertices the face is triangulated. The circum-center is used
* for the regular star faces of uniform polyhedra. The in-center is used for
* the pentagram (#79) and hexagram (#77) of the high-density snub duals, and
* for the pentagrams (#40, #58) and hexagram (#52) of the stellated duals
* with configuration (....)/2. Finally, the self-intersection of the crossed
* parallelogram is used for duals with form p q r| with an even denominator.
*
* This method do not work for the hemi-duals, whose faces are not
* star-shaped and have two self-intersections each.
*
* Thus, for each face we need six auxiliary vertices: The self intersections
* and the terminal points of the truncations of the infinite edges. The
* ideal vertices are listed, but are not used by the face-list.
*
* Note that the face of the last dual (#80) is octagonal, and constists of
* two quadrilaterals of the infinite type.
*/
if (*star && P->even != -1)
Malloc(hit, F, int);
for (i = 0; i < F; i++)
if ((!*star &&
(frac(P->n[P->ftype[i]]), frax.d != 1 && frax.d != frax.n - 1)) ||
(*star &&
P->K == 5 &&
(P->D > 30 ||
denominator (P->m[0]) != 1))) {
/* find the center of the face */
double h;
if (!*star && P->hemi && !P->ftype[i])
h = 0;
else
h = P->minr / dot(f[i],f[i]);
push_point(result, scale (h, f[i]));
} else if (*star && P->even != -1) {
/* find the self-intersection of a crossed parallelogram.
* hit is set if v0v1 intersects v2v3*/
Vector v0, v1, v2, v3, c0, c1, p;
double d0, d1;
v0 = v[P->incid[0][i]];
v1 = v[P->incid[1][i]];
v2 = v[P->incid[2][i]];
v3 = v[P->incid[3][i]];
d0 = sqrt(dot(diff(v0, v2), diff(v0, v2)));
d1 = sqrt(dot (diff(v1, v3), diff(v1, v3)));
c0 = scale(d1, sum(v0, v2));
c1 = scale(d0, sum(v1, v3));
p = scale(0.5 / (d0 + d1), sum(c0, c1));
push_point (result, p);
p = cross(diff(p, v2), diff(p, v3));
hit[i] = (dot(p, p) < 1e-6);
} else if (*star && P->hemi && P->index != last_uniform - 1) {
/* find the terminal points of the truncation and the
* self-intersections.
* v23 v0 v21
* | \ / \ / |
* | v0123 v0321 |
* | / \ / \ |
* v01 v2 v03
*/
Vector v0, v1, v2, v3, v01, v03, v21, v23, v0123, v0321 ;
Vector u;
double t = 1.5;/* truncation adjustment factor */
j = !P->ftype[P->incid[0][i]];
v0 = v[P->incid[j][i]];/* real vertex */
v1 = v[P->incid[j+1][i]];/* ideal vertex (unit vector) */
v2 = v[P->incid[j+2][i]];/* real */
v3 = v[P->incid[(j+3)%4][i]];/* ideal */
/* compute intersections
* this uses the following linear algebra:
* v0123 = v0 + a v1 = v2 + b v3
* v0 x v3 + a (v1 x v3) = v2 x v3
* a (v1 x v3) = (v2 - v0) x v3
* a (v1 x v3) . (v1 x v3) = (v2 - v0) x v3 . (v1 x v3)
*/
u = cross(v1, v3);
v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
v1));
v0321 = sum(v0, scale(dot(cross(diff(v0, v2), v1), u) / dot(u,u),
v3));
/* compute truncations */
v01 = sum(v0 , scale(t, diff(v0123, v0)));
v23 = sum(v2 , scale(t, diff(v0123, v2)));
v03 = sum(v0 , scale(t, diff(v0321, v0)));
v21 = sum(v2 , scale(t, diff(v0321, v2)));
push_point(result, v01);
push_point(result, v23);
push_point(result, v0123);
push_point(result, v03);
push_point(result, v21);
push_point(result, v0321);
} else if (*star && P->index == last_uniform - 1) {
/* find the terminal points of the truncation and the
* self-intersections.
* v23 v0 v21
* | \ / \ / |
* | v0123 v0721 |
* | / \ / \ |
* v01 v2 v07
*
* v65 v4 v67
* | \ / \ / |
* | v4365 v4567 |
* | / \ / \ |
* v43 v6 v45
*/
Vector v0, v1, v2, v3, v4, v5, v6, v7, v01, v07, v21, v23;
Vector v43, v45, v65, v67, v0123, v0721, v4365, v4567;
double t = 1.5;/* truncation adjustment factor */
Vector u;
for (j = 0; j < 8; j++)
if (P->ftype[P->incid[j][i]] == 3)
break;
v0 = v[P->incid[j][i]];/* real {5/3} */
v1 = v[P->incid[(j+1)%8][i]];/* ideal */
v2 = v[P->incid[(j+2)%8][i]];/* real {3} */
v3 = v[P->incid[(j+3)%8][i]];/* ideal */
v4 = v[P->incid[(j+4)%8][i]];/* real {5/2} */
v5 = v[P->incid[(j+5)%8][i]];/* ideal */
v6 = v[P->incid[(j+6)%8][i]];/* real {3/2} */
v7 = v[P->incid[(j+7)%8][i]];/* ideal */
/* compute intersections */
u = cross(v1, v3);
v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
v1));
u = cross(v7, v1);
v0721 = sum(v0, scale(dot(cross(diff(v2, v0), v1), u) / dot(u,u),
v7));
u = cross(v5, v7);
v4567 = sum(v4, scale(dot(cross(diff(v6, v4), v7), u) / dot(u,u),
v5));
u = cross(v3, v5);
v4365 = sum(v4, scale(dot(cross(diff(v6, v4), v5), u) / dot(u,u),
v3));
/* compute truncations */
v01 = sum(v0 , scale(t, diff(v0123, v0)));
v23 = sum(v2 , scale(t, diff(v0123, v2)));
v07 = sum(v0 , scale(t, diff(v0721, v0)));
v21 = sum(v2 , scale(t, diff(v0721, v2)));
v45 = sum(v4 , scale(t, diff(v4567, v4)));
v67 = sum(v6 , scale(t, diff(v4567, v6)));
v43 = sum(v4 , scale(t, diff(v4365, v4)));
v65 = sum(v6 , scale(t, diff(v4365, v6)));
push_point(result, v01);
push_point(result, v23);
push_point(result, v0123);
push_point(result, v07);
push_point(result, v21);
push_point(result, v0721);
push_point(result, v45);
push_point(result, v67);
push_point(result, v4567);
push_point(result, v43);
push_point(result, v65);
push_point(result, v4365);
}
/*
* Face list:
* Each face is printed in a separate line, by listing the indices of its
* vertices. In the non-simple case, the polygon is represented by the
* triangulation, each triangle consists of two polyhedron vertices and one
* auxiliary vertex.
*/
Malloc (result->faces, F * 10, face);
result->nfaces = 0;
ii = V;
facelets = 0;
for (i = 0; i < F; i++) {
if (*star) {
if (P->K == 5 &&
(P->D > 30 ||
denominator (P->m[0]) != 1)) {
for (j = 0; j < P->M - 1; j++) {
push_face3 (result, P->incid[j][i], P->incid[j+1][i], ii);
facelets++;
}
push_face3 (result, P->incid[j][i], P->incid[0][i], ii++);
facelets++;
} else if (P->even != -1) {
if (hit && hit[i]) {
push_face3 (result, P->incid[3][i], P->incid[0][i], ii);
push_face3 (result, P->incid[1][i], P->incid[2][i], ii);
} else {
push_face3 (result, P->incid[0][i], P->incid[1][i], ii);
push_face3 (result, P->incid[2][i], P->incid[3][i], ii);
}
ii++;
facelets += 2;
} else if (P->hemi && P->index != last_uniform - 1) {
j = !P->ftype[P->incid[0][i]];
push_face3 (result, ii, ii + 1, ii + 2);
push_face4 (result, P->incid[j][i], ii + 2, P->incid[j+2][i], ii + 5);
push_face3 (result, ii + 3, ii + 4, ii + 5);
ii += 6;
facelets += 3;
} else if (P->index == last_uniform - 1) {
for (j = 0; j < 8; j++)
if (P->ftype[P->incid[j][i]] == 3)
break;
push_face3 (result, ii, ii + 1, ii + 2);
push_face4 (result,
P->incid[j][i], ii + 2, P->incid[(j+2)%8][i], ii + 5);
push_face3 (result, ii + 3, ii + 4, ii + 5);
push_face3 (result, ii + 6, ii + 7, ii + 8);
push_face4 (result,
P->incid[(j+4)%8][i], ii + 8, P->incid[(j+6)%8][i],
ii + 11);
push_face3 (result, ii + 9, ii + 10, ii + 11);
ii += 12;
facelets += 6;
} else {
result->faces[result->nfaces].npoints = P->M;
Malloc (result->faces[result->nfaces].points, P->M, int);
for (j = 0; j < P->M; j++)
result->faces[result->nfaces].points[j] = P->incid[j][i];
result->nfaces++;
facelets++;
}
} else {
int split = (frac(P->n[P->ftype[i]]),
frax.d != 1 && frax.d != frax.n - 1);
for (j = 0; j < V; j++) {
for (k = 0; k < P->M; k++)
if (P->incid[k][j] == i)
break;
if (k != P->M)
break;
}
if (split) {
ll = j;
for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
for (k = 0; k < P->M; k++)
if (P->incid[k][l] == i)
break;
if (P->adj[k][l] == ll)
k = mod(k + 1 , P->M);
push_face3 (result, ll, l, ii);
facelets++;
ll = l;
}
push_face3 (result, ll, j, ii++);
facelets++;
} else {
int *pp;
int pi = 0;
Malloc (pp, 100, int);
pp[pi++] = j;
ll = j;
for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
for (k = 0; k < P->M; k++)
if (P->incid[k][l] == i)
break;
if (P->adj[k][l] == ll)
k = mod(k + 1 , P->M);
pp[pi++] = l;
ll = l;
}
result->faces[result->nfaces].npoints = pi;
result->faces[result->nfaces].points = pp;
result->nfaces++;
facelets++;
}
}
}
/*
* Face color indices - for polyhedra with multiple face types
* For non-simple faces, the index is repeated as many times as needed by the
* triangulation.
*/
{
int ff = 0;
if (!*star && P->N != 1) {
for (i = 0; i < F; i++)
if (frac(P->n[P->ftype[i]]), frax.d == 1 || frax.d == frax.n - 1)
result->faces[ff++].color = P->ftype[i];
else
for (j = 0; j < frax.n; j++)
result->faces[ff++].color = P->ftype[i];
} else {
for (i = 0; i < facelets; i++)
result->faces[ff++].color = 0;
}
}
if (*star && P->even != -1)
free(hit);
free(v);
free(f);
return result;
}
/* External interface (jwz)
*/
void
free_polyhedron (polyhedron *p)
{
if (!p) return;
Free (p->wythoff);
Free (p->name);
Free (p->group);
Free (p->class);
if (p->faces)
{
int i;
for (i = 0; i < p->nfaces; i++)
Free (p->faces[i].points);
Free (p->faces);
}
Free (p);
}
int
construct_polyhedra (polyhedron ***polyhedra_ret)
{
double freeze = 0;
double azimuth = AZ;
double elevation = EL;
int index = 0;
int count = 0;
polyhedron **result;
Malloc (result, last_uniform * 2 + 3, polyhedron*);
while (index < last_uniform) {
char sym[4];
Polyhedron *P;
sprintf(sym, "#%d", index + 1);
if (!(P = kaleido(sym, 1, 0, 0, 0))) {
Err (strerror(errno));
}
result[count++] = construct_polyhedron (P, P->v, P->V, P->f, P->F,
P->name, P->dual_name,
P->class, "",
azimuth, elevation, freeze);
result[count++] = construct_polyhedron (P, P->f, P->F, P->v, P->V,
P->dual_name, P->name,
P->dual_class, "*",
azimuth, elevation, freeze);
polyfree(P);
index++;
}
*polyhedra_ret = result;
count++; /* leave room for teapot */
return count;
}