/*****************************************************************************
 * #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(void)
{
  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);
  Free (p->points);
  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;
}
