git.hungrycats.org Git - xscreensaver/blob - hacks/glx/polyhedra.c

   1 /*****************************************************************************
   2  * #ident "Id: main.c,v 3.27 2002-01-06 16:23:01+02 rl Exp "
   3  * kaleido
   4  *
   5  *      Kaleidoscopic construction of uniform polyhedra
   6  *      Copyright (c) 1991-2002 Dr. Zvi Har'El <rl@math.technion.ac.il>
   7  *
   8  *      Redistribution and use in source and binary forms, with or without
   9  *      modification, are permitted provided that the following conditions
  10  *      are met:
  11  *
  12  *      1. Redistributions of source code must retain the above copyright
  13  *         notice, this list of conditions and the following disclaimer.
  14  *
  15  *      2. Redistributions in binary form must reproduce the above copyright
  16  *         notice, this list of conditions and the following disclaimer in
  17  *         the documentation and/or other materials provided with the
  18  *         distribution.
  19  *
  20  *      3. The end-user documentation included with the redistribution,
  21  *         if any, must include the following acknowledgment:
  22  *              "This product includes software developed by
  23  *               Dr. Zvi Har'El (http://www.math.technion.ac.il/~rl/)."
  24  *         Alternately, this acknowledgment may appear in the software itself,
  25  *         if and wherever such third-party acknowledgments normally appear.
  26  *
  27  *      This software is provided 'as-is', without any express or implied
  28  *      warranty.  In no event will the author be held liable for any
  29  *      damages arising from the use of this software.
  30  *
  31  *      Author:
  32  *              Dr. Zvi Har'El,
  33  *              Deptartment of Mathematics,
  34  *              Technion, Israel Institue of Technology,
  35  *              Haifa 32000, Israel.
  36  *              E-Mail: rl@math.technion.ac.il
  37  *
  38  *      ftp://ftp.math.technion.ac.il/kaleido/
  39  *      http://www.mathconsult.ch/showroom/unipoly/
  40  *
  41  *      Adapted for xscreensaver by Jamie Zawinski <jwz@jwz.org> 25-Apr-2004
  42  *
  43  *****************************************************************************
  44  */
  45
  46 #include <stdio.h>
  47 #include <ctype.h>
  48 #include <string.h>
  49 #include <math.h>
  50 #include <stdlib.h>
  51 #include <errno.h>
  52
  53 #include "config.h"
  54 #include "polyhedra.h"
  55
  56 extern const char *progname;
  57
  58 #ifndef MAXLONG
  59 #define MAXLONG 0x7FFFFFFF
  60 #endif
  61 #ifndef MAXDIGITS
  62 #define MAXDIGITS 10 /* (int)log10((double)MAXLONG) + 1 */
  63 #endif
  64
  65 #ifndef DBL_EPSILON
  66 #define DBL_EPSILON 2.2204460492503131e-16
  67 #endif
  68 #define BIG_EPSILON 3e-2
  69 #define AZ M_PI/7  /* axis azimuth */
  70 #define EL M_PI/17 /* axis elevation */
  71
  72 #define Err(x) do {\
  73         fprintf (stderr, "%s: %s\n", progname, (x)); \
  74         exit (1); \
  75     } while(0)
  76
  77 #define Free(lvalue) do {\
  78         if (lvalue) {\
  79             free((char*) lvalue);\
  80             lvalue=0;\
  81         }\
  82     } while(0)
  83
  84 #define Matfree(lvalue,n) do {\
  85         if (lvalue) \
  86             matfree((char*) lvalue, n);\
  87         lvalue=0;\
  88     } while(0)
  89
  90 #define Malloc(lvalue,n,type) do {\
  91         if (!(lvalue = (type*) malloc((n) * sizeof(type)))) \
  92             Err("out of memory");\
  93     } while(0)
  94
  95 #define Realloc(lvalue,n,type) do {\
  96         if (!(lvalue = (type*) realloc(lvalue, (n) * sizeof(type)))) \
  97             Err("out of memory");\
  98     } while(0)
  99
 100 #define Calloc(lvalue,n,type) do {\
 101         if (!(lvalue = (type*) calloc(n, sizeof(type))))\
 102             Err("out of memory");\
 103     } while(0)
 104
 105 #define Matalloc(lvalue,n,m,type) do {\
 106         if (!(lvalue = (type**) matalloc(n, (m) * sizeof(type))))\
 107             Err("out of memory");\
 108     } while(0)
 109
 110 #define Sprintfrac(lvalue,x) do {\
 111         if (!(lvalue=sprintfrac(x)))\
 112             return 0;\
 113     } while(0)
 114
 115 #define numerator(x) (frac(x), frax.n)
 116 #define denominator(x) (frac(x), frax.d)
 117 #define compl(x) (frac(x), (double) frax.n / (frax.n-frax.d))
 118
 119 typedef struct {
 120     double x, y, z;
 121 } Vector;
 122
 123 typedef struct {
 124   /* NOTE: some of the int's can be replaced by short's, char's,
 125      or even bit fields, at the expense of readability!!!*/
 126   int index; /* index to the standard list, the array uniform[] */
 127   int N; /* number of faces types (atmost 5)*/
 128   int M; /* vertex valency  (may be big for dihedral polyhedra) */
 129   int V; /* vertex count */
 130   int E; /* edge count */
 131   int F; /* face count */
 132   int D; /* density */
 133   int chi; /* Euler characteristic */
 134   int g; /* order of symmetry group */
 135   int K; /* symmetry type: D=2, T=3, O=4, I=5 */
 136   int hemi;/* flag hemi polyhedron */
 137   int onesided;/* flag onesided polyhedron */
 138   int even; /* removed face in pqr| */
 139   int *Fi; /* face counts by type (array N)*/
 140   int *rot; /* vertex configuration (array M of 0..N-1) */
 141   int *snub; /* snub triangle configuration (array M of 0..1) */
 142   int *firstrot; /* temporary for vertex generation (array V) */
 143   int *anti; /* temporary for direction of ideal vertices (array E) */
 144   int *ftype; /* face types (array F) */
 145   int **e; /* edges (matrix 2 x E of 0..V-1)*/
 146   int **dual_e; /* dual edges (matrix 2 x E of 0..F-1)*/
 147   int **incid; /* vertex-face incidence (matrix M x V of 0..F-1)*/
 148   int **adj; /* vertex-vertex adjacency (matrix M x V of 0..V-1)*/
 149   double p[4]; /* p, q and r; |=0 */
 150   double minr; /* smallest nonzero inradius */
 151   double gon; /* basis type for dihedral polyhedra */
 152   double *n; /* number of side of a face of each type (array N) */
 153   double *m; /* number of faces at a vertex of each type (array N) */
 154   double *gamma; /* fundamental angles in radians (array N) */
 155   char *polyform; /* printable Wythoff symbol */
 156   char *config; /* printable vertex configuration */
 157   char *group; /* printable group name */
 158   char *name; /* name, standard or manifuctured */
 159   char *dual_name; /* dual name, standard or manifuctured */
 160   char *class;
 161   char *dual_class;
 162   Vector *v; /* vertex coordinates (array V) */
 163   Vector *f; /* face coordinates (array F)*/
 164 } Polyhedron;
 165
 166 typedef struct {
 167   long n,d;
 168 } Fraction;
 169
 170 static Polyhedron *polyalloc(void);
 171 static Vector rotate(Vector vertex, Vector axis, double angle);
 172
 173 static Vector sum3(Vector a, Vector b, Vector c);
 174 static Vector scale(double k, Vector a);
 175 static Vector sum(Vector a, Vector b);
 176 static Vector diff(Vector a, Vector b);
 177 static Vector pole (double r, Vector a, Vector b, Vector c);
 178 static Vector cross(Vector a, Vector b);
 179 static double dot(Vector a, Vector b);
 180 static int same(Vector a, Vector b, double epsilon);
 181
 182 static char *sprintfrac(double x);
 183
 184 static void frac(double x);
 185 static void matfree(void *mat, int rows);
 186 static void *matalloc(int rows, int row_size);
 187
 188 static Fraction frax;
 189
 190
 191 static struct {
 192   char *Wythoff, *name, *dual, *group, *class, *dual_class;
 193   short Coxeter, Wenninger;
 194 } uniform[] = {
 195
 196   /****************************************************************************
 197    *            Dihedral Schwarz Triangles (D5 only)
 198    ***************************************************************************/
 199
 200                                                            /* (2 2 5) (D1/5) */
 201   /*  1 */      {"2 5|2",       "Pentagonal Prism",
 202                                 "Pentagonal Dipyramid",
 203                                 "Dihedral (D[1/5])",
 204                                 "",
 205                                 "",
 206                                 0, 0},
 207
 208   /*  2 */      {"|2 2 5",      "Pentagonal Antiprism",
 209                                 "Pentagonal Deltohedron",
 210                                 "Dihedral (D[1/5])",
 211                                 "",
 212                                 "",
 213                                 0, 0},
 214                                                          /* (2 2 5/2) (D2/5) */
 215   /*  3 */      {"2 5/2|2",     "Pentagrammic Prism",
 216                                 "Pentagrammic Dipyramid",
 217                                 "Dihedral (D[2/5])",
 218                                 "",
 219                                 "",
 220                                 0, 0},
 221
 222   /*  4 */      {"|2 2 5/2",    "Pentagrammic Antiprism",
 223                                 "Pentagrammic Deltohedron",
 224                                 "Dihedral (D[2/5])",
 225                                 "",
 226                                 "",
 227                                 0, 0},
 228                                                          /* (5/3 2 2) (D3/5) */
 229
 230   /*  5 */      {"|2 2 5/3",    "Pentagrammic Crossed Antiprism",
 231                                 "Pentagrammic Concave Deltohedron",
 232                                 "Dihedral (D[3/5])",
 233                                 "",
 234                                 "",
 235                                 0, 0},
 236
 237   /****************************************************************************
 238    *            Tetrahedral
 239    ***************************************************************************/
 240
 241                                                              /* (2 3 3) (T1) */
 242   /*  6 */      {"3|2 3",       "Tetrahedron",
 243                                 "Tetrahedron",
 244                                 "Tetrahedral (T[1])",
 245                                 "Platonic Solid",
 246                                 "Platonic Solid",
 247                                 15, 1},
 248
 249   /*  7 */      {"2 3|3",       "Truncated Tetrahedron",
 250                                 "Triakistetrahedron",
 251                                 "Tetrahedral (T[1])",
 252                                 "Archimedian Solid",
 253                                 "Catalan Solid",
 254                                 16, 6},
 255                                                            /* (3/2 3 3) (T2) */
 256   /*  8 */      {"3/2 3|3",     "Octahemioctahedron",
 257                                 "Octahemioctacron",
 258                                 "Tetrahedral (T[2])",
 259                                 "",
 260                                 "",
 261                                 37, 68},
 262
 263                                                            /* (3/2 2 3) (T3) */
 264   /*  9 */      {"3/2 3|2",     "Tetrahemihexahedron",
 265                                 "Tetrahemihexacron",
 266                                 "Tetrahedral (T[3])",
 267                                 "",
 268                                 "",
 269                                 36, 67},
 270
 271   /****************************************************************************
 272    *            Octahedral
 273    ***************************************************************************/
 274
 275                                                              /* (2 3 4) (O1) */
 276   /* 10 */      {"4|2 3",       "Octahedron",
 277                                 "Cube",
 278                                 "Octahedral (O[1])",
 279                                 "Platonic Solid",
 280                                 "Platonic Solid",
 281                                 17, 2},
 282
 283   /* 11 */      {"3|2 4",       "Cube",
 284                                 "Octahedron",
 285                                 "Octahedral (O[1])",
 286                                 "Platonic Solid",
 287                                 "Platonic Solid",
 288                                 18, 3},
 289
 290   /* 12 */      {"2|3 4",       "Cuboctahedron",
 291                                 "Rhombic Dodecahedron",
 292                                 "Octahedral (O[1])",
 293                                 "Archimedian Solid",
 294                                 "Catalan Solid",
 295                                 19, 11},
 296
 297   /* 13 */      {"2 4|3",       "Truncated Octahedron",
 298                                 "Tetrakishexahedron",
 299                                 "Octahedral (O[1])",
 300                                 "Archimedian Solid",
 301                                 "Catalan Solid",
 302                                 20, 7},
 303
 304   /* 14 */      {"2 3|4",       "Truncated Cube",
 305                                 "Triakisoctahedron",
 306                                 "Octahedral (O[1])",
 307                                 "Archimedian Solid",
 308                                 "Catalan Solid",
 309                                 21, 8},
 310
 311   /* 15 */      {"3 4|2",       "Rhombicuboctahedron",
 312                                 "Deltoidal Icositetrahedron",
 313                                 "Octahedral (O[1])",
 314                                 "Archimedian Solid",
 315                                 "Catalan Solid",
 316                                 22, 13},
 317
 318   /* 16 */      {"2 3 4|",      "Truncated Cuboctahedron",
 319                                 "Disdyakisdodecahedron",
 320                                 "Octahedral (O[1])",
 321                                 "Archimedian Solid",
 322                                 "Catalan Solid",
 323                                 23, 15},
 324
 325   /* 17 */      {"|2 3 4",      "Snub Cube",
 326                                 "Pentagonal Icositetrahedron",
 327                                 "Octahedral (O[1])",
 328                                 "Archimedian Solid",
 329                                 "Catalan Solid",
 330                                 24, 17},
 331                                                           /* (3/2 4 4) (O2b) */
 332
 333   /* 18 */      {"3/2 4|4",     "Small Cubicuboctahedron",
 334                                 "Small Hexacronic Icositetrahedron",
 335                                 "Octahedral (O[2b])",
 336                                 "",
 337                                 "",
 338                                 38, 69},
 339                                                            /* (4/3 3 4) (O4) */
 340
 341   /* 19 */      {"3 4|4/3",     "Great Cubicuboctahedron",
 342                                 "Great Hexacronic Icositetrahedron",
 343                                 "Octahedral (O[4])",
 344                                 "",
 345                                 "",
 346                                 50, 77},
 347
 348   /* 20 */      {"4/3 4|3",     "Cubohemioctahedron",
 349                                 "Hexahemioctacron",
 350                                 "Octahedral (O[4])",
 351                                 "",
 352                                 "",
 353                                 51, 78},
 354
 355   /* 21 */      {"4/3 3 4|",    "Cubitruncated Cuboctahedron",
 356                                 "Tetradyakishexahedron",
 357                                 "Octahedral (O[4])",
 358                                 "",
 359                                 "",
 360                                 52, 79},
 361                                                            /* (3/2 2 4) (O5) */
 362
 363   /* 22 */      {"3/2 4|2",     "Great Rhombicuboctahedron",
 364                                 "Great Deltoidal Icositetrahedron",
 365                                 "Octahedral (O[5])",
 366                                 "",
 367                                 "",
 368                                 59, 85},
 369
 370   /* 23 */      {"3/2 2 4|",    "Small Rhombihexahedron",
 371                                 "Small Rhombihexacron",
 372                                 "Octahedral (O[5])",
 373                                 "",
 374                                 "",
 375                                 60, 86},
 376                                                            /* (4/3 2 3) (O7) */
 377
 378   /* 24 */      {"2 3|4/3",     "Stellated Truncated Hexahedron",
 379                                 "Great Triakisoctahedron",
 380                                 "Octahedral (O[7])",
 381                                 "",
 382                                 "",
 383                                 66, 92},
 384
 385   /* 25 */      {"4/3 2 3|",    "Great Truncated Cuboctahedron",
 386                                 "Great Disdyakisdodecahedron",
 387                                 "Octahedral (O[7])",
 388                                 "",
 389                                 "",
 390                                 67, 93},
 391                                                         /* (4/3 3/2 2) (O11) */
 392
 393   /* 26 */      {"4/3 3/2 2|",  "Great Rhombihexahedron",
 394                                 "Great Rhombihexacron",
 395                                 "Octahedral (O[11])",
 396                                 "",
 397                                 "",
 398                                 82, 103},
 399
 400   /****************************************************************************
 401    *            Icosahedral
 402    ***************************************************************************/
 403
 404                                                              /* (2 3 5) (I1) */
 405   /* 27 */      {"5|2 3",       "Icosahedron",
 406                                 "Dodecahedron",
 407                                 "Icosahedral (I[1])",
 408                                 "Platonic Solid",
 409                                 "Platonic Solid",
 410                                 25, 4},
 411
 412   /* 28 */      {"3|2 5",       "Dodecahedron",
 413                                 "Icosahedron",
 414                                 "Icosahedral (I[1])",
 415                                 "Platonic Solid",
 416                                 "Platonic Solid",
 417                                 26, 5},
 418
 419   /* 29 */      {"2|3 5",       "Icosidodecahedron",
 420                                 "Rhombic Triacontahedron",
 421                                 "Icosahedral (I[1])",
 422                                 "Archimedian Solid",
 423                                 "Catalan Solid",
 424                                 28, 12},
 425
 426   /* 30 */      {"2 5|3",       "Truncated Icosahedron",
 427                                 "Pentakisdodecahedron",
 428                                 "Icosahedral (I[1])",
 429                                 "Archimedian Solid",
 430                                 "Catalan Solid",
 431                                 27, 9},
 432
 433   /* 31 */      {"2 3|5",       "Truncated Dodecahedron",
 434                                 "Triakisicosahedron",
 435                                 "Icosahedral (I[1])",
 436                                 "Archimedian Solid",
 437                                 "Catalan Solid",
 438                                 29, 10},
 439
 440   /* 32 */      {"3 5|2",       "Rhombicosidodecahedron",
 441                                 "Deltoidal Hexecontahedron",
 442                                 "Icosahedral (I[1])",
 443                                 "Archimedian Solid",
 444                                 "Catalan Solid",
 445                                 30, 14},
 446
 447   /* 33 */      {"2 3 5|",      "Truncated Icosidodechedon",
 448                                 "Disdyakistriacontahedron",
 449                                 "Icosahedral (I[1])",
 450                                 "Archimedian Solid",
 451                                 "Catalan Solid",
 452                                 31, 16},
 453
 454   /* 34 */      {"|2 3 5",      "Snub Dodecahedron",
 455                                 "Pentagonal Hexecontahedron",
 456                                 "Icosahedral (I[1])",
 457                                 "Archimedian Solid",
 458                                 "Catalan Solid",
 459                                 32, 18},
 460                                                           /* (5/2 3 3) (I2a) */
 461
 462   /* 35 */      {"3|5/2 3",     "Small Ditrigonal Icosidodecahedron",
 463                                 "Small Triambic Icosahedron",
 464                                 "Icosahedral (I[2a])",
 465                                 "",
 466                                 "",
 467                                 39, 70},
 468
 469   /* 36 */      {"5/2 3|3",     "Small Icosicosidodecahedron",
 470                                 "Small Icosacronic Hexecontahedron",
 471                                 "Icosahedral (I[2a])",
 472                                 "",
 473                                 "",
 474                                 40, 71},
 475
 476   /* 37 */      {"|5/2 3 3",    "Small Snub Icosicosidodecahedron",
 477                                 "Small Hexagonal Hexecontahedron",
 478                                 "Icosahedral (I[2a])",
 479                                 "",
 480                                 "",
 481                                 41, 110},
 482                                                           /* (3/2 5 5) (I2b) */
 483
 484   /* 38 */      {"3/2 5|5",     "Small Dodecicosidodecahedron",
 485                                 "Small Dodecacronic Hexecontahedron",
 486                                 "Icosahedral (I[2b])",
 487                                 "",
 488                                 "",
 489                                 42, 72},
 490                                                            /* (2 5/2 5) (I3) */
 491
 492   /* 39 */      {"5|2 5/2",     "Small Stellated Dodecahedron",
 493                                 "Great Dodecahedron",
 494                                 "Icosahedral (I[3])",
 495                                 "Truncated Kepler-Poinsot Solid",
 496                                 "",
 497                                 43, 20},
 498
 499   /* 40 */      {"5/2|2 5",     "Great Dodecahedron",
 500                                 "Small Stellated Dodecahedron",
 501                                 "Icosahedral (I[3])",
 502                                 "",
 503                                 "",
 504                                 44, 21},
 505
 506   /* 41 */      {"2|5/2 5",     "Great Dodecadodecahedron",
 507                                 "Medial Rhombic Triacontahedron",
 508                                 "Icosahedral (I[3])",
 509                                 "",
 510                                 "",
 511                                 45, 73},
 512
 513   /* 42 */      {"2 5/2|5",     "Truncated Great Dodecahedron",
 514                                 "Small Stellapentakisdodecahedron",
 515                                 "Icosahedral (I[3])",
 516                                 "Truncated Kepler-Poinsot Solid",
 517                                 "",
 518                                 47, 75},
 519
 520   /* 43 */      {"5/2 5|2",     "Rhombidodecadodecahedron",
 521                                 "Medial Deltoidal Hexecontahedron",
 522                                 "Icosahedral (I[3])",
 523                                 "",
 524                                 "",
 525                                 48, 76},
 526
 527   /* 44 */      {"2 5/2 5|",    "Small Rhombidodecahedron",
 528                                 "Small Rhombidodecacron",
 529                                 "Icosahedral (I[3])",
 530                                 "",
 531                                 "",
 532                                 46, 74},
 533
 534   /* 45 */      {"|2 5/2 5",    "Snub Dodecadodecahedron",
 535                                 "Medial Pentagonal Hexecontahedron",
 536                                 "Icosahedral (I[3])",
 537                                 "",
 538                                 "",
 539                                 49, 111},
 540                                                            /* (5/3 3 5) (I4) */
 541
 542   /* 46 */      {"3|5/3 5",     "Ditrigonal Dodecadodecahedron",
 543                                 "Medial Triambic Icosahedron",
 544                                 "Icosahedral (I[4])",
 545                                 "",
 546                                 "",
 547                                 53, 80},
 548
 549   /* 47 */      {"3 5|5/3",     "Great Ditrigonal Dodecicosidodecahedron",
 550                                "Great Ditrigonal Dodecacronic Hexecontahedron",
 551                                 "Icosahedral (I[4])",
 552                                 "",
 553                                 "",
 554                                 54, 81},
 555
 556   /* 48 */      {"5/3 3|5",     "Small Ditrigonal Dodecicosidodecahedron",
 557                                "Small Ditrigonal Dodecacronic Hexecontahedron",
 558                                 "Icosahedral (I[4])",
 559                                 "",
 560                                 "",
 561                                 55, 82},
 562
 563   /* 49 */      {"5/3 5|3",     "Icosidodecadodecahedron",
 564                                 "Medial Icosacronic Hexecontahedron",
 565                                 "Icosahedral (I[4])",
 566                                 "",
 567                                 "",
 568                                 56, 83},
 569
 570   /* 50 */      {"5/3 3 5|",    "Icositruncated Dodecadodecahedron",
 571                                 "Tridyakisicosahedron",
 572                                 "Icosahedral (I[4])",
 573                                 "",
 574                                 "",
 575                                 57, 84},
 576
 577   /* 51 */      {"|5/3 3 5",    "Snub Icosidodecadodecahedron",
 578                                 "Medial Hexagonal Hexecontahedron",
 579                                 "Icosahedral (I[4])",
 580                                 "",
 581                                 "",
 582                                 58, 112},
 583                                                           /* (3/2 3 5) (I6b) */
 584
 585   /* 52 */      {"3/2|3 5",     "Great Ditrigonal Icosidodecahedron",
 586                                 "Great Triambic Icosahedron",
 587                                 "Icosahedral (I[6b])",
 588                                 "",
 589                                 "",
 590                                 61, 87},
 591
 592   /* 53 */      {"3/2 5|3",     "Great Icosicosidodecahedron",
 593                                 "Great Icosacronic Hexecontahedron",
 594                                 "Icosahedral (I[6b])",
 595                                 "",
 596                                 "",
 597                                 62, 88},
 598
 599   /* 54 */      {"3/2 3|5",     "Small Icosihemidodecahedron",
 600                                 "Small Icosihemidodecacron",
 601                                 "Icosahedral (I[6b])",
 602                                 "",
 603                                 "",
 604                                 63, 89},
 605
 606   /* 55 */      {"3/2 3 5|",    "Small Dodecicosahedron",
 607                                 "Small Dodecicosacron",
 608                                 "Icosahedral (I[6b])",
 609                                 "",
 610                                 "",
 611                                 64, 90},
 612                                                           /* (5/4 5 5) (I6c) */
 613
 614   /* 56 */      {"5/4 5|5",     "Small Dodecahemidodecahedron",
 615                                 "Small Dodecahemidodecacron",
 616                                 "Icosahedral (I[6c])",
 617                                 "",
 618                                 "",
 619                                 65, 91},
 620                                                            /* (2 5/2 3) (I7) */
 621
 622   /* 57 */      {"3|2 5/2",     "Great Stellated Dodecahedron",
 623                                 "Great Icosahedron",
 624                                 "Icosahedral (I[7])",
 625                                 "",
 626                                 "",
 627                                 68, 22},
 628
 629   /* 58 */      {"5/2|2 3",     "Great Icosahedron",
 630                                 "Great Stellated Dodecahedron",
 631                                 "Icosahedral (I[7])",
 632                                 "",
 633                                 "",
 634                                 69, 41},
 635
 636   /* 59 */      {"2|5/2 3",     "Great Icosidodecahedron",
 637                                 "Great Rhombic Triacontahedron",
 638                                 "Icosahedral (I[7])",
 639                                 "Truncated Kepler-Poinsot Solid",
 640                                 "",
 641                                 70, 94},
 642
 643   /* 60 */      {"2 5/2|3",     "Great Truncated Icosahedron",
 644                                 "Great Stellapentakisdodecahedron",
 645                                 "Icosahedral (I[7])",
 646                                 "Truncated Kepler-Poinsot Solid",
 647                                 "",
 648                                 71, 95},
 649
 650   /* 61 */      {"2 5/2 3|",    "Rhombicosahedron",
 651                                 "Rhombicosacron",
 652                                 "Icosahedral (I[7])",
 653                                 "",
 654                                 "",
 655                                 72, 96},
 656
 657   /* 62 */      {"|2 5/2 3",    "Great Snub Icosidodecahedron",
 658                                 "Great Pentagonal Hexecontahedron",
 659                                 "Icosahedral (I[7])",
 660                                 "",
 661                                 "",
 662                                 73, 113},
 663                                                            /* (5/3 2 5) (I9) */
 664
 665   /* 63 */      {"2 5|5/3",     "Small Stellated Truncated Dodecahedron",
 666                                 "Great Pentakisdodekahedron",
 667                                 "Icosahedral (I[9])",
 668                                 "",
 669                                 "",
 670                                 74, 97},
 671
 672   /* 64 */      {"5/3 2 5|",    "Truncated Dodecadodecahedron",
 673                                 "Medial Disdyakistriacontahedron",
 674                                 "Icosahedral (I[9])",
 675                                 "",
 676                                 "",
 677                                 75, 98},
 678
 679   /* 65 */      {"|5/3 2 5",    "Inverted Snub Dodecadodecahedron",
 680                                 "Medial Inverted Pentagonal Hexecontahedron",
 681                                 "Icosahedral (I[9])",
 682                                 "",
 683                                 "",
 684                                 76, 114},
 685                                                        /* (5/3 5/2 3) (I10a) */
 686
 687   /* 66 */      {"5/2 3|5/3",   "Great Dodecicosidodecahedron",
 688                                 "Great Dodecacronic Hexecontahedron",
 689                                 "Icosahedral (I[10a])",
 690                                 "",
 691                                 "",
 692                                 77, 99},
 693
 694   /* 67 */      {"5/3 5/2|3",   "Small Dodecahemicosahedron",
 695                                 "Small Dodecahemicosacron",
 696                                 "Icosahedral (I[10a])",
 697                                 "",
 698                                 "",
 699                                 78, 100},
 700
 701   /* 68 */      {"5/3 5/2 3|",  "Great Dodecicosahedron",
 702                                 "Great Dodecicosacron",
 703                                 "Icosahedral (I[10a])",
 704                                 "",
 705                                 "",
 706                                 79, 101},
 707
 708   /* 69 */      {"|5/3 5/2 3",  "Great Snub Dodecicosidodecahedron",
 709                                 "Great Hexagonal Hexecontahedron",
 710                                 "Icosahedral (I[10a])",
 711                                 "",
 712                                 "",
 713                                 80, 115},
 714                                                          /* (5/4 3 5) (I10b) */
 715
 716   /* 70 */      {"5/4 5|3",     "Great Dodecahemicosahedron",
 717                                 "Great Dodecahemicosacron",
 718                                 "Icosahedral (I[10b])",
 719                                 "",
 720                                 "",
 721                                 81, 102},
 722                                                           /* (5/3 2 3) (I13) */
 723
 724   /* 71 */      {"2 3|5/3",     "Great Stellated Truncated Dodecahedron",
 725                                 "Great Triakisicosahedron",
 726                                 "Icosahedral (I[13])",
 727                                 "",
 728                                 "",
 729                                 83, 104},
 730
 731   /* 72 */      {"5/3 3|2",     "Great Rhombicosidodecahedron",
 732                                 "Great Deltoidal Hexecontahedron",
 733                                 "Icosahedral (I[13])",
 734                                 "",
 735                                 "",
 736                                 84, 105},
 737
 738   /* 73 */      {"5/3 2 3|",    "Great Truncated Icosidodecahedron",
 739                                 "Great Disdyakistriacontahedron",
 740                                 "Icosahedral (I[13])",
 741                                 "",
 742                                 "",
 743                                 87, 108},
 744
 745   /* 74 */      {"|5/3 2 3",    "Great Inverted Snub Icosidodecahedron",
 746                                 "Great Inverted Pentagonal Hexecontahedron",
 747                                 "Icosahedral (I[13])",
 748                                 "",
 749                                 "",
 750                                 88, 116},
 751                                                      /* (5/3 5/3 5/2) (I18a) */
 752
 753   /* 75 */      {"5/3 5/2|5/3", "Great Dodecahemidodecahedron",
 754                                 "Great Dodecahemidodecacron",
 755                                 "Icosahedral (I[18a])",
 756                                 "",
 757                                 "",
 758                                 86, 107},
 759                                                        /* (3/2 5/3 3) (I18b) */
 760
 761   /* 76 */      {"3/2 3|5/3",   "Great Icosihemidodecahedron",
 762                                 "Great Icosihemidodecacron",
 763                                 "Icosahedral (I[18b])",
 764                                 "",
 765                                 "",
 766                                 85, 106},
 767                                                       /* (3/2 3/2 5/3) (I22) */
 768
 769   /* 77 */      {"|3/2 3/2 5/2","Small Retrosnub Icosicosidodecahedron",
 770                                 "Small Hexagrammic Hexecontahedron",
 771                                 "Icosahedral (I[22])",
 772                                 "",
 773                                 "",
 774                                 91, 118},
 775                                                         /* (3/2 5/3 2) (I23) */
 776
 777   /* 78 */      {"3/2 5/3 2|",  "Great Rhombidodecahedron",
 778                                 "Great Rhombidodecacron",
 779                                 "Icosahedral (I[23])",
 780                                 "",
 781                                 "",
 782                                 89, 109},
 783
 784   /* 79 */      {"|3/2 5/3 2",  "Great Retrosnub Icosidodecahedron",
 785                                 "Great Pentagrammic Hexecontahedron",
 786                                 "Icosahedral (I[23])",
 787                                 "",
 788                                 "",
 789                                 90, 117},
 790
 791   /****************************************************************************
 792    *            Last But Not Least
 793    ***************************************************************************/
 794
 795   /* 80 */    {"3/2 5/3 3 5/2", "Great Dirhombicosidodecahedron",
 796                                 "Great Dirhombicosidodecacron",
 797                                 "Non-Wythoffian",
 798                                 "",
 799                                 "",
 800                                 92, 119}
 801 };
 802
 803 static int last_uniform = sizeof (uniform) / sizeof (uniform[0]);
 804
 805
 806
 807 static int unpacksym(char *sym, Polyhedron *P);
 808 static int moebius(Polyhedron *P);
 809 static int decompose(Polyhedron *P);
 810 static int guessname(Polyhedron *P);
 811 static int newton(Polyhedron *P, int need_approx);
 812 static int exceptions(Polyhedron *P);
 813 static int count(Polyhedron *P);
 814 static int configuration(Polyhedron *P);
 815 static int vertices(Polyhedron *P);
 816 static int faces(Polyhedron *P);
 817 static int edgelist(Polyhedron *P);
 818
 819 static Polyhedron *
 820 kaleido(char *sym, int need_coordinates, int need_edgelist, int need_approx,
 821         int just_list)
 822 {
 823   Polyhedron *P;
 824   /*
 825    * Allocate a Polyhedron structure P.
 826    */
 827   if (!(P = polyalloc()))
 828     return 0;
 829   /*
 830    * Unpack input symbol into P.
 831    */
 832   if (!unpacksym(sym, P))
 833     return 0;
 834   /*
 835    * Find Mebius triangle, its density and Euler characteristic.
 836    */
 837   if (!moebius(P))
 838     return 0;
 839   /*
 840    * Decompose Schwarz triangle.
 841    */
 842   if (!decompose(P))
 843     return 0;
 844   /*
 845    * Find the names of the polyhedron and its dual.
 846    */
 847   if (!guessname(P))
 848     return 0;
 849   if (just_list)
 850     return P;
 851   /*
 852    * Solve Fundamental triangles, optionally printing approximations.
 853    */
 854   if (!newton(P,need_approx))
 855     return 0;
 856   /*
 857    * Deal with exceptional polyhedra.
 858    */
 859   if (!exceptions(P))
 860     return 0;
 861   /*
 862    * Count edges and faces, update density and characteristic if needed.
 863    */
 864   if (!count(P))
 865     return 0;
 866   /*
 867    * Generate printable vertex configuration.
 868    */
 869   if (!configuration(P))
 870     return 0;
 871   /*
 872    * Compute coordinates.
 873    */
 874   if (!need_coordinates && !need_edgelist)
 875     return P;
 876   if (!vertices(P))
 877     return 0;
 878   if (!faces (P))
 879     return 0;
 880   /*
 881    * Compute edgelist.
 882    */
 883   if (!need_edgelist)
 884     return P;
 885   if (!edgelist(P))
 886     return 0;
 887   return P;
 888 }
 889
 890 /*
 891  * Allocate a blank Polyhedron structure and initialize some of its nonblank
 892  * fields.
 893  *
 894  * Array and matrix field are allocated when needed.
 895  */
 896 static Polyhedron *
 897 polyalloc()
 898 {
 899   Polyhedron *P;
 900   Calloc(P, 1, Polyhedron);
 901   P->index = -1;
 902   P->even = -1;
 903   P->K = 2;
 904   return P;
 905 }
 906
 907 /*
 908  * Free the struture allocated by polyalloc(), as well as all the array and
 909  * matrix fields.
 910  */
 911 static void
 912 polyfree(Polyhedron *P)
 913 {
 914   Free(P->Fi);
 915   Free(P->n);
 916   Free(P->m);
 917   Free(P->gamma);
 918   Free(P->rot);
 919   Free(P->snub);
 920   Free(P->firstrot);
 921   Free(P->anti);
 922   Free(P->ftype);
 923   Free(P->polyform);
 924   Free(P->config);
 925   if (P->index < 0) {
 926     Free(P->name);
 927     Free(P->dual_name);
 928   }
 929   Free(P->v);
 930   Free(P->f);
 931   Matfree(P->e, 2);
 932   Matfree(P->dual_e, 2);
 933   Matfree(P->incid, P->M);
 934   Matfree(P->adj, P->M);
 935   free(P);
 936 }
 937
 938 static void *
 939 matalloc(int rows, int row_size)
 940 {
 941   void **mat;
 942   int i = 0;
 943   if (!(mat = malloc(rows * sizeof (void *))))
 944     return 0;
 945   while ((mat[i] = malloc(row_size)) && ++i < rows)
 946     ;
 947   if (i == rows)
 948     return (void *)mat;
 949   while (--i >= 0)
 950     free(mat[i]);
 951   free(mat);
 952   return 0;
 953 }
 954
 955 static void
 956 matfree(void *mat, int rows)
 957 {
 958   while (--rows >= 0)
 959     free(((void **)mat)[rows]);
 960   free(mat);
 961 }
 962
 963 /*
 964  * compute the mathematical modulus function.
 965  */
 966 static int
 967 mod (int i, int j)
 968 {
 969   return (i%=j)>=0?i:j<0?i-j:i+j;
 970 }
 971
 972
 973 /*
 974  * Find the numerator and the denominator using the Euclidean algorithm.
 975  */
 976 static void
 977 frac(double x)
 978 {
 979   static Fraction zero = {0,1}, inf = {1,0};
 980   Fraction r0, r;
 981   long f;
 982   double s = x;
 983   r = zero;
 984   frax = inf;
 985   for (;;) {
 986     if (fabs(s) > (double) MAXLONG)
 987       return;
 988     f = (long) floor (s);
 989     r0 = r;
 990     r = frax;
 991     frax.n = frax.n * f + r0.n;
 992     frax.d = frax.d * f + r0.d;
 993     if (x == (double)frax.n/(double)frax.d)
 994       return;
 995     s = 1 / (s - f);
 996   }
 997 }
 998
 999
1000 /*
1001  * Unpack input symbol: Wythoff symbol or an index to uniform[]. The symbol is
1002  * a # followed by a number, or a three fractions and a bar in some order. We
1003  * allow no bars only if it result from the input symbol #80.
1004  */
1005 static int
1006 unpacksym(char *sym, Polyhedron *P)
1007 {
1008   int i = 0, n, d, bars = 0;
1009   char c;
1010   while ((c = *sym++) && isspace(c))
1011     ;
1012   if (!c) Err("no data");
1013   if (c == '#') {
1014     while ((c = *sym++) && isspace(c))
1015       ;
1016     if (!c)
1017       Err("no digit after #");
1018     if (!isdigit(c))
1019       Err("not a digit");
1020     n = c - '0';
1021     while ((c = *sym++) && isdigit(c))
1022       n = n * 10 + c - '0';
1023     if (!n)
1024       Err("zero index");
1025     if (n > last_uniform)
1026       Err("index too big");
1027     sym--;
1028     while ((c = *sym++) && isspace(c))
1029       ;
1030     if (c)
1031       Err("data exceeded");
1032     sym = uniform[P->index = n - 1].Wythoff;
1033   } else
1034     sym--;
1035
1036   for (;;) {
1037     while ((c = *sym++) && isspace(c))
1038       ;
1039     if (!c) {
1040       if (i == 4 && (bars || P->index == last_uniform - 1))
1041         return 1;
1042       if (!bars)
1043         Err("no bars");
1044       Err("not enough fractions");
1045     }
1046     if (i == 4)
1047       Err("data exceeded");
1048     if (c == '|'){
1049       if (++bars > 1)
1050         Err("too many bars");
1051       P->p[i++] = 0;
1052       continue;
1053     }
1054     if (!isdigit(c))
1055       Err("not a digit");
1056     n = c - '0';
1057     while ((c = *sym++) && isdigit(c))
1058       n = n * 10 + c - '0';
1059     if (c && isspace (c))
1060       while ((c = *sym++) && isspace(c))
1061         ;
1062     if (c != '/') {
1063       sym--;
1064       if ((P->p[i++] = n) <= 1)
1065         Err("fraction<=1");
1066       continue;
1067     }
1068     while ((c = *sym++) && isspace(c))
1069       ;
1070     if (!c || !isdigit(c))
1071       return 0;
1072     d = c - '0';
1073     while ((c = *sym++) && isdigit(c))
1074       d = d * 10 + c - '0';
1075     if (!d)
1076       Err("zero denominator");
1077     sym--;
1078     if ((P->p[i++] = (double) n / d) <= 1)
1079       Err("fraction<=1");
1080   }
1081 }
1082
1083 /*
1084  * Using Wythoff symbol (p|qr, pq|r, pqr| or |pqr), find the Moebius triangle
1085  * (2 3 K) (or (2 2 n)) of the Schwarz triangle (pqr), the order g of its
1086  * symmetry group, its Euler characteristic chi, and its covering density D.
1087  * g is the number of copies of (2 3 K) covering the sphere, i.e.,
1088  *
1089  *              g * pi * (1/2 + 1/3 + 1/K - 1) = 4 * pi
1090  *
1091  * D is the number of times g copies of (pqr) cover the sphere, i.e.
1092  *
1093  *              D * 4 * pi = g * pi * (1/p + 1/q + 1/r - 1)
1094  *
1095  * chi is V - E + F, where F = g is the number of triangles, E = 3*g/2 is the
1096  * number of triangle edges, and V = Vp+ Vq+ Vr, with Vp = g/(2*np) being the
1097  * number of vertices with angle pi/p (np is the numerator of p).
1098  */
1099 static int
1100 moebius(Polyhedron *P)
1101 {
1102   int twos = 0, j, len = 1;
1103   /*
1104    * Arrange Wythoff symbol in a presentable form. In the same time check the
1105    * restrictions on the three fractions: They all have to be greater then one,
1106    * and the numerators 4 or 5 cannot occur together.  We count the ocurrences
1107    * of 2 in `two', and save the largest numerator in `P->K', since they
1108    * reflect on the symmetry group.
1109    */
1110   P->K = 2;
1111   if (P->index == last_uniform - 1) {
1112     Malloc(P->polyform, ++len, char);
1113     strcpy(P->polyform, "|");
1114   } else
1115     Calloc(P->polyform, len, char);
1116   for (j = 0; j < 4; j++) {
1117     if (P->p[j]) {
1118       char *s;
1119       Sprintfrac(s, P->p[j]);
1120       if (j && P->p[j-1]) {
1121         Realloc(P->polyform, len += strlen (s) + 1, char);
1122         strcat(P->polyform, " ");
1123       } else
1124         Realloc (P->polyform, len += strlen (s), char);
1125       strcat(P->polyform, s);
1126       free(s);
1127       if (P->p[j] != 2) {
1128         int k;
1129         if ((k = numerator (P->p[j])) > P->K) {
1130           if (P->K == 4)
1131             break;
1132           P->K = k;
1133         } else if (k < P->K && k == 4)
1134           break;
1135       } else
1136         twos++;
1137     } else  {
1138       Realloc(P->polyform, ++len, char);
1139       strcat(P->polyform, "|");
1140     }
1141   }
1142   /*
1143    * Find the symmetry group P->K (where 2, 3, 4, 5 represent the dihedral,
1144    * tetrahedral, octahedral and icosahedral groups, respectively), and its
1145    * order P->g.
1146    */
1147   if (twos >= 2) {/* dihedral */
1148     P->g = 4 * P->K;
1149     P->K = 2;
1150   } else {
1151     if (P->K > 5)
1152       Err("numerator too large");
1153     P->g = 24 * P->K / (6 - P->K);
1154   }
1155   /*
1156    * Compute the nominal density P->D and Euler characteristic P->chi.
1157    * In few exceptional cases, these values will be modified later.
1158    */
1159   if (P->index != last_uniform - 1) {
1160     int i;
1161     P->D = P->chi = - P->g;
1162     for (j = 0; j < 4; j++) if (P->p[j]) {
1163       P->chi += i = P->g / numerator(P->p[j]);
1164       P->D += i * denominator(P->p[j]);
1165     }
1166     P->chi /= 2;
1167     P->D /= 4;
1168     if (P->D <= 0)
1169       Err("nonpositive density");
1170   }
1171   return 1;
1172 }
1173
1174 /*
1175  * Decompose Schwarz triangle into N right triangles and compute the vertex
1176  * count V and the vertex valency M.  V is computed from the number g of
1177  * Schwarz triangles in the cover, divided by the number of triangles which
1178  * share a vertex. It is halved for one-sided polyhedra, because the
1179  * kaleidoscopic construction really produces a double orientable covering of
1180  * such polyhedra. All q' q|r are of the "hemi" type, i.e. have equatorial {2r}
1181  * faces, and therefore are (except 3/2 3|3 and the dihedra 2 2|r) one-sided. A
1182  * well known example is 3/2 3|4, the "one-sided heptahedron". Also, all p q r|
1183  * with one even denominator have a crossed parallelogram as a vertex figure,
1184  * and thus are one-sided as well.
1185  */
1186 static int
1187 decompose(Polyhedron *P)
1188 {
1189   int j, J, *s, *t;
1190   if (!P->p[1]) { /* p|q r */
1191     P->N = 2;
1192     P->M = 2 * numerator(P->p[0]);
1193     P->V = P->g / P->M;
1194     Malloc(P->n, P->N, double);
1195     Malloc(P->m, P->N, double);
1196     Malloc(P->rot, P->M, int);
1197     s = P->rot;
1198     for (j = 0; j < 2; j++) {
1199       P->n[j] = P->p[j+2];
1200       P->m[j] = P->p[0];
1201     }
1202     for (j = P->M / 2; j--;) {
1203       *s++ = 0;
1204       *s++ = 1;
1205     }
1206   } else if (!P->p[2]) { /* p q|r */
1207     P->N = 3;
1208     P->M = 4;
1209     P->V = P->g / 2;
1210     Malloc(P->n, P->N, double);
1211     Malloc(P->m, P->N, double);
1212     Malloc(P->rot, P->M, int);
1213     s = P->rot;
1214     P->n[0] = 2 * P->p[3];
1215     P->m[0] = 2;
1216     for (j = 1; j < 3; j++) {
1217       P->n[j] = P->p[j-1];
1218       P->m[j] = 1;
1219       *s++ = 0;
1220       *s++ = j;
1221     }
1222     if (fabs(P->p[0] - compl (P->p[1])) < DBL_EPSILON) {/* p = q' */
1223       /* P->p[0]==compl(P->p[1]) should work.  However, MSDOS
1224        * yeilds a 7e-17 difference! Reported by Jim Buddenhagen
1225        * <jb1556@daditz.sbc.com> */
1226       P->hemi = 1;
1227       P->D = 0;
1228       if (P->p[0] != 2 && !(P->p[3] == 3 && (P->p[0] == 3 ||
1229                                              P->p[1] == 3))) {
1230         P->onesided = 1;
1231         P->V /= 2;
1232         P->chi /= 2;
1233       }
1234     }
1235   } else if (!P->p[3]) { /* p q r| */
1236     P->M = P->N = 3;
1237     P->V = P->g;
1238     Malloc(P->n, P->N, double);
1239     Malloc(P->m, P->N, double);
1240     Malloc(P->rot, P->M, int);
1241     s = P->rot;
1242     for (j = 0; j < 3; j++) {
1243       if (!(denominator(P->p[j]) % 2)) {
1244         /* what happens if there is more then one even denominator? */
1245         if (P->p[(j+1)%3] != P->p[(j+2)%3]) { /* needs postprocessing */
1246           P->even = j;/* memorize the removed face */
1247           P->chi -= P->g / numerator(P->p[j]) / 2;
1248           P->onesided = 1;
1249           P->D = 0;
1250         } else {/* for p = q we get a double 2 2r|p */
1251                 /* noted by Roman Maeder <maeder@inf.ethz.ch> for 4 4 3/2| */
1252                 /* Euler characteristic is still wrong */
1253           P->D /= 2;
1254         }
1255         P->V /= 2;
1256       }
1257       P->n[j] = 2 * P->p[j];
1258       P->m[j] = 1;
1259       *s++ = j;
1260     }
1261   } else { /* |p q r - snub polyhedron */
1262     P->N = 4;
1263     P->M = 6;
1264     P->V = P->g / 2;/* Only "white" triangles carry a vertex */
1265     Malloc(P->n, P->N, double);
1266     Malloc(P->m, P->N, double);
1267     Malloc(P->rot, P->M, int);
1268     Malloc(P->snub, P->M, int);
1269     s = P->rot;
1270     t = P->snub;
1271     P->m[0] = P->n[0] = 3;
1272     for (j = 1; j < 4; j++) {
1273       P->n[j] = P->p[j];
1274       P->m[j] = 1;
1275       *s++ = 0;
1276       *s++ = j;
1277       *t++ = 1;
1278       *t++ = 0;
1279     }
1280   }
1281   /*
1282    * Sort the fundamental triangles (using bubble sort) according to decreasing
1283    * n[i], while pushing the trivial triangles (n[i] = 2) to the end.
1284    */
1285   J = P->N - 1;
1286   while (J) {
1287     int last;
1288     last = J;
1289     J = 0;
1290     for (j = 0; j < last; j++) {
1291       if ((P->n[j] < P->n[j+1] || P->n[j] == 2) && P->n[j+1] != 2) {
1292         int i;
1293         double temp;
1294         temp = P->n[j];
1295         P->n[j] = P->n[j+1];
1296         P->n[j+1] = temp;
1297         temp = P->m[j];
1298         P->m[j] = P->m[j+1];
1299         P->m[j+1] = temp;
1300         for (i = 0; i < P->M; i++) {
1301           if (P->rot[i] == j)
1302             P->rot[i] = j+1;
1303           else if (P->rot[i] == j+1)
1304             P->rot[i] = j;
1305         }
1306         if (P->even != -1) {
1307           if (P->even == j)
1308             P->even = j+1;
1309           else if (P->even == j+1)
1310             P->even = j;
1311         }
1312         J = j;
1313       }
1314     }
1315   }
1316   /*
1317    *  Get rid of repeated triangles.
1318    */
1319   for (J = 0; J < P->N && P->n[J] != 2;J++) {
1320     int k, i;
1321     for (j = J+1; j < P->N && P->n[j]==P->n[J]; j++)
1322       P->m[J] += P->m[j];
1323     k = j - J - 1;
1324     if (k) {
1325       for (i = j; i < P->N; i++) {
1326         P->n[i - k] = P->n[i];
1327         P->m[i - k] = P->m[i];
1328       }
1329       P->N -= k;
1330       for (i = 0; i < P->M; i++) {
1331         if (P->rot[i] >= j)
1332           P->rot[i] -= k;
1333         else if (P->rot[i] > J)
1334           P->rot[i] = J;
1335       }
1336       if (P->even >= j)
1337         P->even -= k;
1338     }
1339   }
1340   /*
1341    * Get rid of trivial triangles.
1342    */
1343   if (!J)
1344     J = 1; /* hosohedron */
1345   if (J < P->N) {
1346     int i;
1347     P->N = J;
1348     for (i = 0; i < P->M; i++) {
1349       if (P->rot[i] >= P->N) {
1350         for (j = i + 1; j < P->M; j++) {
1351           P->rot[j-1] = P->rot[j];
1352           if (P->snub)
1353             P->snub[j-1] = P->snub[j];
1354         }
1355         P->M--;
1356       }
1357     }
1358   }
1359   /*
1360    * Truncate arrays
1361    */
1362   Realloc(P->n, P->N, double);
1363   Realloc(P->m, P->N, double);
1364   Realloc(P->rot, P->M, int);
1365   if (P->snub)
1366     Realloc(P->snub, P->M, int);
1367   return 1;
1368 }
1369
1370
1371 static int dihedral(Polyhedron *P, char *name, char *dual_name);
1372
1373
1374 /*
1375  * Get the polyhedron name, using standard list or guesswork. Ideally, we
1376  * should try to locate the Wythoff symbol in the standard list (unless, of
1377  * course, it is dihedral), after doing few normalizations, such as sorting
1378  * angles and splitting isoceles triangles.
1379  */
1380 static int
1381 guessname(Polyhedron *P)
1382 {
1383   if (P->index != -1) {/* tabulated */
1384     P->name = uniform[P->index].name;
1385     P->dual_name = uniform[P->index].dual;
1386     P->group = uniform[P->index].group;
1387     P->class = uniform[P->index].class;
1388     P->dual_class = uniform[P->index].dual_class;
1389     return 1;
1390   } else if (P->K == 2) {/* dihedral nontabulated */
1391     if (!P->p[0]) {
1392       if (P->N == 1) {
1393         Malloc(P->name, sizeof ("Octahedron"), char);
1394         Malloc(P->dual_name, sizeof ("Cube"), char);
1395         strcpy(P->name, "Octahedron");
1396         strcpy(P->dual_name, "Cube");
1397         return 1;
1398       }
1399       P->gon = P->n[0] == 3 ? P->n[1] : P->n[0];
1400       if (P->gon >= 2)
1401         return dihedral(P, "Antiprism", "Deltohedron");
1402       else
1403         return dihedral(P, "Crossed Antiprism", "Concave Deltohedron");
1404     } else if (!P->p[3] ||
1405                (!P->p[2] &&
1406                 P->p[3] == 2)) {
1407       if (P->N == 1) {
1408         Malloc(P->name, sizeof("Cube"), char);
1409         Malloc(P->dual_name, sizeof("Octahedron"), char);
1410         strcpy(P->name, "Cube");
1411         strcpy(P->dual_name, "Octahedron");
1412         return 1;
1413       }
1414       P->gon = P->n[0] == 4 ? P->n[1] : P->n[0];
1415       return dihedral(P, "Prism", "Dipyramid");
1416     } else if (!P->p[1] && P->p[0] != 2) {
1417       P->gon = P->m[0];
1418       return dihedral(P, "Hosohedron", "Dihedron");
1419     } else {
1420       P->gon = P->n[0];
1421       return dihedral(P, "Dihedron", "Hosohedron");
1422     }
1423   } else {/* other nontabulated */
1424     static char *pre[] = {"Tetr", "Oct", "Icos"};
1425     Malloc(P->name, 50, char);
1426     Malloc(P->dual_name, 50, char);
1427     sprintf(P->name, "%sahedral ", pre[P->K - 3]);
1428     if (P->onesided)
1429       strcat (P->name, "One-Sided ");
1430     else if (P->D == 1)
1431       strcat(P->name, "Convex ");
1432     else
1433       strcat(P->name, "Nonconvex ");
1434     strcpy(P->dual_name, P->name);
1435     strcat(P->name, "Isogonal Polyhedron");
1436     strcat(P->dual_name, "Isohedral Polyhedron");
1437     Realloc(P->name, strlen (P->name) + 1, char);
1438     Realloc(P->dual_name, strlen (P->dual_name) + 1, char);
1439     return 1;
1440   }
1441 }
1442
1443 static int
1444 dihedral(Polyhedron *P, char *name, char *dual_name)
1445 {
1446   char *s;
1447   int i;
1448   Sprintfrac(s, P->gon < 2 ? compl (P->gon) : P->gon);
1449   i = strlen(s) + sizeof ("-gonal ");
1450   Malloc(P->name, i + strlen (name), char);
1451   Malloc(P->dual_name, i + strlen (dual_name), char);
1452   sprintf(P->name, "%s-gonal %s", s, name);
1453   sprintf(P->dual_name, "%s-gonal %s", s, dual_name);
1454   free(s);
1455   return 1;
1456 }
1457
1458 /*
1459  * Solve the fundamental right spherical triangles.
1460  * If need_approx is set, print iterations on standard error.
1461  */
1462 static int
1463 newton(Polyhedron *P, int need_approx)
1464 {
1465   /*
1466    * First, we find initial approximations.
1467    */
1468   int j;
1469   double cosa;
1470   Malloc(P->gamma, P->N, double);
1471   if (P->N == 1) {
1472     P->gamma[0] = M_PI / P->m[0];
1473     return 1;
1474   }
1475   for (j = 0; j < P->N; j++)
1476     P->gamma[j] = M_PI / 2 - M_PI / P->n[j];
1477   errno = 0; /* may be non-zero from some reason */
1478   /*
1479    * Next, iteratively find closer approximations for gamma[0] and compute
1480    * other gamma[j]'s from Napier's equations.
1481    */
1482   if (need_approx)
1483     fprintf(stderr, "Solving %s\n", P->polyform);
1484   for (;;) {
1485     double delta = M_PI, sigma = 0;
1486     for (j = 0; j < P->N; j++) {
1487       if (need_approx)
1488         fprintf(stderr, "%-20.15f", P->gamma[j]);
1489       delta -= P->m[j] * P->gamma[j];
1490     }
1491     if (need_approx)
1492       printf("(%g)\n", delta);
1493     if (fabs(delta) < 11 * DBL_EPSILON)
1494       return 1;
1495     /* On a RS/6000, fabs(delta)/DBL_EPSILON may occilate between 8 and
1496      * 10. Reported by David W. Sanderson <dws@ssec.wisc.edu> */
1497     for (j = 0; j < P->N; j++)
1498       sigma += P->m[j] * tan(P->gamma[j]);
1499     P->gamma[0] += delta * tan(P->gamma[0]) / sigma;
1500     if (P->gamma[0] < 0 || P->gamma[0] > M_PI)
1501       Err("gamma out of bounds");
1502     cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
1503     for (j = 1; j < P->N; j++)
1504       P->gamma[j] = asin(cos(M_PI / P->n[j]) / cosa);
1505     if (errno)
1506       Err(strerror(errno));
1507   }
1508 }
1509
1510 /*
1511  * Postprocess pqr| where r has an even denominator (cf. Coxeter &al. Sec.9).
1512  * Remove the {2r} and add a retrograde {2p} and retrograde {2q}.
1513  */
1514 static int
1515 exceptions(Polyhedron *P)
1516 {
1517   int j;
1518   if (P->even != -1) {
1519     P->M = P->N = 4;
1520     Realloc(P->n, P->N, double);
1521     Realloc(P->m, P->N, double);
1522     Realloc(P->gamma, P->N, double);
1523     Realloc(P->rot, P->M, int);
1524     for (j = P->even + 1; j < 3; j++) {
1525       P->n[j-1] = P->n[j];
1526       P->gamma[j-1] = P->gamma[j];
1527     }
1528     P->n[2] = compl(P->n[1]);
1529     P->gamma[2] = - P->gamma[1];
1530     P->n[3] = compl(P->n[0]);
1531     P->m[3] = 1;
1532     P->gamma[3] = - P->gamma[0];
1533     P->rot[0] = 0;
1534     P->rot[1] = 1;
1535     P->rot[2] = 3;
1536     P->rot[3] = 2;
1537   }
1538
1539   /*
1540    * Postprocess the last polyhedron |3/2 5/3 3 5/2 by taking a |5/3 3 5/2,
1541    * replacing the three snub triangles by four equatorial squares and adding
1542    * the missing {3/2} (retrograde triangle, cf. Coxeter &al. Sec. 11).
1543    */
1544   if (P->index == last_uniform - 1) {
1545     P->N = 5;
1546     P->M = 8;
1547     Realloc(P->n, P->N, double);
1548     Realloc(P->m, P->N, double);
1549     Realloc(P->gamma, P->N, double);
1550     Realloc(P->rot, P->M, int);
1551     Realloc(P->snub, P->M, int);
1552     P->hemi = 1;
1553     P->D = 0;
1554     for (j = 3; j; j--) {
1555       P->m[j] = 1;
1556       P->n[j] = P->n[j-1];
1557       P->gamma[j] = P->gamma[j-1];
1558     }
1559     P->m[0] = P->n[0] = 4;
1560     P->gamma[0] = M_PI / 2;
1561     P->m[4] = 1;
1562     P->n[4] = compl(P->n[1]);
1563     P->gamma[4] = - P->gamma[1];
1564     for (j = 1; j < 6; j += 2) P->rot[j]++;
1565     P->rot[6] = 0;
1566     P->rot[7] = 4;
1567     P->snub[6] = 1;
1568     P->snub[7] = 0;
1569   }
1570   return 1;
1571 }
1572
1573 /*
1574  * Compute edge and face counts, and update D and chi.  Update D in the few
1575  * cases the density of the polyhedron is meaningful but different than the
1576  * density of the corresponding Schwarz triangle (cf. Coxeter &al., p. 418 and
1577  * p. 425).
1578  * In these cases, spherical faces of one type are concave (bigger than a
1579  * hemisphere), and the actual density is the number of these faces less the
1580  * computed density.  Note that if j != 0, the assignment gamma[j] = asin(...)
1581  * implies gamma[j] cannot be obtuse.  Also, compute chi for the only
1582  * non-Wythoffian polyhedron.
1583  */
1584 static int
1585 count(Polyhedron *P)
1586 {
1587   int j, temp;
1588   Malloc(P->Fi, P->N, int);
1589   for (j = 0; j < P->N; j++) {
1590     P->E += temp = P->V * numerator(P->m[j]);
1591     P->F += P->Fi[j] = temp / numerator(P->n[j]);
1592   }
1593   P->E /= 2;
1594   if (P->D && P->gamma[0] > M_PI / 2)
1595     P->D = P->Fi[0] - P->D;
1596   if (P->index == last_uniform - 1)
1597     P->chi = P->V - P->E + P->F;
1598   return 1;
1599 }
1600
1601 /*
1602  * Generate a printable vertex configuration symbol.
1603  */
1604 static int
1605 configuration(Polyhedron *P)
1606 {
1607   int j, len = 2;
1608   for (j = 0; j < P->M; j++) {
1609     char *s;
1610     Sprintfrac(s, P->n[P->rot[j]]);
1611     len += strlen (s) + 2;
1612     if (!j) {
1613       Malloc(P->config, len, char);
1614 /*      strcpy(P->config, "(");*/
1615       strcpy(P->config, "");
1616     } else {
1617       Realloc(P->config, len, char);
1618       strcat(P->config, ", ");
1619     }
1620     strcat(P->config, s);
1621     free(s);
1622   }
1623 /*  strcat (P->config, ")");*/
1624   if ((j = denominator (P->m[0])) != 1) {
1625     char s[MAXDIGITS + 2];
1626     sprintf(s, "/%d", j);
1627     Realloc(P->config, len + strlen (s), char);
1628     strcat(P->config, s);
1629   }
1630   return 1;
1631 }
1632
1633 /*
1634  * Compute polyhedron vertices and vertex adjecency lists.
1635  * The vertices adjacent to v[i] are v[adj[0][i], v[adj[1][i], ...
1636  * v[adj[M-1][i], ordered counterclockwise.  The algorith is a BFS on the
1637  * vertices, in such a way that the vetices adjacent to a givem vertex are
1638  * obtained from its BFS parent by a cyclic sequence of rotations. firstrot[i]
1639  * points to the first  rotaion in the sequence when applied to v[i]. Note that
1640  * for non-snub polyhedra, the rotations at a child are opposite in sense when
1641  * compared to the rotations at the parent. Thus, we fill adj[*][i] from the
1642  * end to signify clockwise rotations. The firstrot[] array is not needed for
1643  * display thus it is freed after being used for face computations below.
1644  */
1645 static int
1646 vertices(Polyhedron *P)
1647 {
1648   int i, newV = 2;
1649   double cosa;
1650   Malloc(P->v, P->V, Vector);
1651   Matalloc(P->adj, P->M, P->V, int);
1652   Malloc(P->firstrot, P->V, int); /* temporary , put in Polyhedron
1653                                      structure so that may be freed on
1654                                      error */
1655   cosa = cos(M_PI / P->n[0]) / sin(P->gamma[0]);
1656   P->v[0].x = 0;
1657   P->v[0].y = 0;
1658   P->v[0].z = 1;
1659   P->firstrot[0] = 0;
1660   P->adj[0][0] = 1;
1661   P->v[1].x = 2 * cosa * sqrt(1 - cosa * cosa);
1662   P->v[1].y = 0;
1663   P->v[1].z = 2 * cosa * cosa - 1;
1664   if (!P->snub) {
1665     P->firstrot[1] = 0;
1666     P->adj[0][1] = -1;/* start the other side */
1667     P->adj[P->M-1][1] = 0;
1668   } else {
1669     P->firstrot[1] = P->snub[P->M-1] ? 0 : P->M-1 ;
1670     P->adj[0][1] = 0;
1671   }
1672   for (i = 0; i < newV; i++) {
1673     int j, k;
1674     int last, one, start, limit;
1675     if (P->adj[0][i] == -1) {
1676       one = -1; start = P->M-2; limit = -1;
1677     } else {
1678       one = 1; start = 1; limit = P->M;
1679     }
1680     k = P->firstrot[i];
1681     for (j = start; j != limit; j += one) {
1682       Vector temp;
1683       int J;
1684       temp = rotate (P->v[P->adj[j-one][i]], P->v[i],
1685                      one * 2 * P->gamma[P->rot[k]]);
1686       for (J=0; J<newV && !same(P->v[J],temp,BIG_EPSILON); J++)
1687         ;/* noop */
1688       P->adj[j][i] = J;
1689       last = k;
1690       if (++k == P->M)
1691         k = 0;
1692       if (J == newV) { /* new vertex */
1693         if (newV == P->V) Err ("too many vertices");
1694         P->v[newV++] = temp;
1695         if (!P->snub) {
1696           P->firstrot[J] = k;
1697           if (one > 0) {
1698             P->adj[0][J] = -1;
1699             P->adj[P->M-1][J] = i;
1700           } else {
1701             P->adj[0][J] = i;
1702           }
1703         } else {
1704           P->firstrot[J] = !P->snub[last] ? last :
1705             !P->snub[k] ? (k+1)%P->M : k ;
1706           P->adj[0][J] = i;
1707         }
1708       }
1709     }
1710   }
1711   return 1;
1712 }
1713
1714 /*
1715  * Compute polyhedron faces (dual vertices) and incidence matrices.
1716  * For orientable polyhedra, we can distinguish between the two faces meeting
1717  * at a given directed edge and identify the face on the left and the face on
1718  * the right, as seen from the outside.  For one-sided polyhedra, the vertex
1719  * figure is a papillon (in Coxeter &al.  terminology, a crossed parallelogram)
1720  * and the two faces meeting at an edge can be identified as the side face
1721  * (n[1] or n[2]) and the diagonal face (n[0] or n[3]).
1722  */
1723 static int
1724 faces(Polyhedron *P)
1725 {
1726   int i, newF = 0;
1727   Malloc (P->f, P->F, Vector);
1728   Malloc (P->ftype, P->F, int);
1729   Matalloc (P->incid, P->M, P->V, int);
1730   P->minr = 1 / fabs (tan (M_PI / P->n[P->hemi]) * tan (P->gamma[P->hemi]));
1731   for (i = P->M; --i>=0;) {
1732     int j;
1733     for (j = P->V; --j>=0;)
1734       P->incid[i][j] = -1;
1735   }
1736   for (i = 0; i < P->V; i++) {
1737     int j;
1738     for (j = 0; j < P->M; j++) {
1739       int i0, J;
1740       int pap=0;/* papillon edge type */
1741       if (P->incid[j][i] != -1)
1742         continue;
1743       P->incid[j][i] = newF;
1744       if (newF == P->F)
1745         Err("too many faces");
1746       P->f[newF] = pole(P->minr, P->v[i], P->v[P->adj[j][i]],
1747                         P->v[P->adj[mod(j + 1, P->M)][i]]);
1748       P->ftype[newF] = P->rot[mod(P->firstrot[i] + ((P->adj[0][i] <
1749                                                      P->adj[P->M - 1][i])
1750                                                     ?  j
1751                                                     : -j - 2),
1752                                   P->M)];
1753       if (P->onesided)
1754         pap = (P->firstrot[i] + j) % 2;
1755       i0 = i;
1756       J = j;
1757       for (;;) {
1758         int k;
1759         k = i0;
1760         if ((i0 = P->adj[J][k]) == i) break;
1761         for (J = 0; J < P->M && P->adj[J][i0] != k; J++)
1762           ;/* noop */
1763         if (J == P->M)
1764           Err("too many faces");
1765         if (P->onesided && (J + P->firstrot[i0]) % 2 == pap) {
1766           P->incid [J][i0] = newF;
1767           if (++J >= P->M)
1768             J = 0;
1769         } else {
1770           if (--J < 0)
1771             J = P->M - 1;
1772           P->incid [J][i0] = newF;
1773         }
1774       }
1775       newF++;
1776     }
1777   }
1778   Free(P->firstrot);
1779   Free(P->rot);
1780   Free(P->snub);
1781   return 1;
1782 }
1783
1784 /*
1785  * Compute edge list and graph polyhedron and dual.
1786  * If the polyhedron is of the "hemi" type, each edge has one finite vertex and
1787  * one ideal vertex. We make sure the latter is always the out-vertex, so that
1788  * the edge becomes a ray (half-line).  Each ideal vertex is represented by a
1789  * unit Vector, and the direction of the ray is either parallel or
1790  * anti-parallel this Vector. We flag this in the array P->anti[E].
1791  */
1792 static int
1793 edgelist(Polyhedron *P)
1794 {
1795   int i, j, *s, *t, *u;
1796   Matalloc(P->e, 2, P->E, int);
1797   Matalloc(P->dual_e, 2, P->E, int);
1798   s = P->e[0];
1799   t = P->e[1];
1800   for (i = 0; i < P->V; i++)
1801     for (j = 0; j < P->M; j++)
1802       if (i < P->adj[j][i]) {
1803         *s++ = i;
1804         *t++ = P->adj[j][i];
1805       }
1806   s = P->dual_e[0];
1807   t = P->dual_e[1];
1808   if (!P->hemi)
1809     P->anti = 0;
1810   else
1811     Malloc(P->anti, P->E, int);
1812   u = P->anti;
1813   for (i = 0; i < P->V; i++)
1814     for (j = 0; j < P->M; j++)
1815       if (i < P->adj[j][i])
1816         {
1817           if (!u) {
1818             *s++ = P->incid[mod(j-1,P->M)][i];
1819             *t++ = P->incid[j][i];
1820           } else {
1821             if (P->ftype[P->incid[j][i]]) {
1822               *s = P->incid[j][i];
1823               *t = P->incid[mod(j-1,P->M)][i];
1824             } else {
1825               *s = P->incid[mod(j-1,P->M)][i];
1826               *t = P->incid[j][i];
1827             }
1828             *u++ = dot(P->f[*s++], P->f[*t++]) > 0;
1829           }
1830         }
1831   return 1;
1832 }
1833
1834
1835 static char *
1836 sprintfrac(double x)
1837 {
1838   char *s;
1839   frac (x);
1840   if (!frax.d) {
1841     Malloc(s, sizeof ("infinity"), char);
1842     strcpy(s, "infinity");
1843   } else if (frax.d == 1) {
1844     char n[MAXDIGITS + 1];
1845     sprintf(n, "%ld", frax.n);
1846     Malloc(s, strlen (n) +  1, char);
1847     strcpy(s, n);
1848   } else {
1849     char n[MAXDIGITS + 1], d[MAXDIGITS + 1];
1850     sprintf(n, "%ld", frax.n);
1851     sprintf(d, "%ld", frax.d);
1852     Malloc(s, strlen (n) + strlen (d) + 2, char);
1853     sprintf(s, "%s/%s", n, d);
1854   }
1855   return s;
1856 }
1857
1858 static double
1859 dot(Vector a, Vector b)
1860 {
1861   return a.x * b.x + a.y * b.y + a.z * b.z;
1862 }
1863
1864 static Vector
1865 scale(double k, Vector a)
1866 {
1867   a.x *= k;
1868   a.y *= k;
1869   a.z *= k;
1870   return a;
1871 }
1872
1873 static Vector
1874 diff(Vector a, Vector b)
1875 {
1876   a.x -= b.x;
1877   a.y -= b.y;
1878   a.z -= b.z;
1879   return a;
1880 }
1881
1882 static Vector
1883 cross(Vector a, Vector b)
1884 {
1885   Vector p;
1886   p.x = a.y * b.z - a.z * b.y;
1887   p.y = a.z * b.x - a.x * b.z;
1888   p.z = a.x * b.y - a.y * b.x;
1889   return p;
1890 }
1891
1892 static Vector
1893 sum(Vector a, Vector b)
1894 {
1895   a.x += b.x;
1896   a.y += b.y;
1897   a.z += b.z;
1898   return a;
1899 }
1900
1901 static Vector
1902 sum3(Vector a, Vector b, Vector c)
1903 {
1904   a.x += b.x + c.x;
1905   a.y += b.y + c.y;
1906   a.z += b.z + c.z;
1907   return a;
1908 }
1909
1910 static Vector
1911 rotate(Vector vertex, Vector axis, double angle)
1912 {
1913   Vector p;
1914   p = scale(dot (axis, vertex), axis);
1915   return sum3(p, scale(cos(angle), diff(vertex, p)),
1916               scale(sin(angle), cross(axis, vertex)));
1917 }
1918
1919 Vector x, y, z;
1920
1921 /*
1922  * rotate the standard frame
1923  */
1924 static void
1925 rotframe(double azimuth, double elevation, double angle)
1926 {
1927   static Vector X = {1,0,0}, Y = {0,1,0}, Z = {0,0,1};
1928   Vector axis;
1929
1930   axis = rotate(rotate (X, Y, elevation), Z, azimuth);
1931   x = rotate(X, axis, angle);
1932   y = rotate(Y, axis, angle);
1933   z = rotate(Z, axis, angle);
1934 }
1935
1936 /*
1937  * rotate an array of n Vectors
1938  */
1939 static void
1940 rotarray(Vector *new, Vector *old, int n)
1941 {
1942   while (n--) {
1943     *new++ = sum3(scale(old->x, x), scale(old->y, y), scale(old->z, z));
1944     old++;
1945   }
1946 }
1947
1948 static int
1949 same(Vector a, Vector b, double epsilon)
1950 {
1951   return fabs(a.x - b.x) < epsilon && fabs(a.y - b.y) < epsilon
1952     && fabs(a.z - b.z) < epsilon;
1953 }
1954
1955 /*
1956  * Compute the polar reciprocal of the plane containing a, b and c:
1957  *
1958  * If this plane does not contain the origin, return p such that
1959  *      dot(p,a) = dot(p,b) = dot(p,b) = r.
1960  *
1961  * Otherwise, return p such that
1962  *      dot(p,a) = dot(p,b) = dot(p,c) = 0
1963  * and
1964  *      dot(p,p) = 1.
1965  */
1966 static Vector
1967 pole(double r, Vector a, Vector b, Vector c)
1968 {
1969   Vector p;
1970   double k;
1971   p = cross(diff(b, a), diff(c, a));
1972   k = dot(p, a);
1973   if (fabs(k) < 1e-6)
1974     return scale(1 / sqrt(dot(p, p)), p);
1975   else
1976     return scale(r/ k , p);
1977 }
1978
1979 \f
1980 /* output */
1981
1982
1983
1984
1985 static void rotframe(double azimuth, double elevation, double angle);
1986 static void rotarray(Vector *new, Vector *old, int n);
1987 static int mod (int i, int j);
1988
1989
1990 static void
1991 push_point (polyhedron *p, Vector v)
1992 {
1993   p->points[p->npoints].x = v.x;
1994   p->points[p->npoints].y = v.y;
1995   p->points[p->npoints].z = v.z;
1996   p->npoints++;
1997 }
1998
1999 static void
2000 push_face3 (polyhedron *p, int x, int y, int z)
2001 {
2002   p->faces[p->nfaces].npoints = 3;
2003   Malloc (p->faces[p->nfaces].points, 3, int);
2004   p->faces[p->nfaces].points[0] = x;
2005   p->faces[p->nfaces].points[1] = y;
2006   p->faces[p->nfaces].points[2] = z;
2007   p->nfaces++;
2008 }
2009
2010 static void
2011 push_face4 (polyhedron *p, int x, int y, int z, int w)
2012 {
2013   p->faces[p->nfaces].npoints = 4;
2014   Malloc (p->faces[p->nfaces].points, 4, int);
2015   p->faces[p->nfaces].points[0] = x;
2016   p->faces[p->nfaces].points[1] = y;
2017   p->faces[p->nfaces].points[2] = z;
2018   p->faces[p->nfaces].points[3] = w;
2019   p->nfaces++;
2020 }
2021
2022
2023
2024
2025 static polyhedron *
2026 construct_polyhedron (Polyhedron *P, Vector *v, int V, Vector *f, int F,
2027                       char *name, char *dual, char *class, char *star,
2028                       double azimuth, double elevation, double freeze)
2029 {
2030   int i, j, k=0, l, ll, ii, *hit=0, facelets;
2031
2032   polyhedron *result;
2033   Vector *temp;
2034
2035   Malloc (result, 1, polyhedron);
2036   memset (result, 0, sizeof(*result));
2037
2038   /*
2039    * Rotate polyhedron
2040    */
2041   rotframe(azimuth, elevation, freeze);
2042   Malloc(temp, V, Vector);
2043   rotarray(temp, v, V);
2044   v = temp;
2045   Malloc(temp, F, Vector);
2046   rotarray(temp, f, F);
2047   f = temp;
2048
2049   result->number = P->index + 1;
2050   result->name = strdup (name);
2051   result->dual = strdup (dual);
2052   result->wythoff = strdup (P->polyform);
2053   result->config = strdup (P->config);
2054   result->group = strdup (P->group);
2055   result->class = strdup (class);
2056
2057   /*
2058    * Vertex list
2059    */
2060   Malloc (result->points, V + F * 13, point);
2061   result->npoints = 0;
2062
2063   result->nedges = P->E;
2064   result->logical_faces = F;
2065   result->logical_vertices = V;
2066   result->density = P->D;
2067   result->chi = P->chi;
2068
2069   for (i = 0; i < V; i++)
2070     push_point (result, v[i]);
2071
2072   /*
2073    * Auxiliary vertices (needed because current VRML browsers cannot handle
2074    * non-simple polygons, i.e., ploygons with self intersections): Each
2075    * non-simple face is assigned an auxiliary vertex. By connecting it to the
2076    * rest of the vertices the face is triangulated. The circum-center is used
2077    * for the regular star faces of uniform polyhedra. The in-center is used for
2078    * the pentagram (#79) and hexagram (#77) of the high-density snub duals, and
2079    * for the pentagrams (#40, #58) and hexagram (#52) of the stellated duals
2080    * with configuration (....)/2. Finally, the self-intersection of the crossed
2081    * parallelogram is used for duals with form p q r| with an even denominator.
2082    *
2083    * This method do not work for the hemi-duals, whose faces are not
2084    * star-shaped and have two self-intersections each.
2085    *
2086    * Thus, for each face we need six auxiliary vertices: The self intersections
2087    * and the terminal points of the truncations of the infinite edges. The
2088    * ideal vertices are listed, but are not used by the face-list.
2089    *
2090    * Note that the face of the last dual (#80) is octagonal, and constists of
2091    * two quadrilaterals of the infinite type.
2092    */
2093
2094   if (*star && P->even != -1)
2095     Malloc(hit, F, int);
2096   for (i = 0; i < F; i++)
2097     if ((!*star &&
2098          (frac(P->n[P->ftype[i]]), frax.d != 1 && frax.d != frax.n - 1)) ||
2099         (*star &&
2100          P->K == 5 &&
2101          (P->D > 30 ||
2102           denominator (P->m[0]) != 1))) {
2103       /* find the center of the face */
2104       double h;
2105       if (!*star && P->hemi && !P->ftype[i])
2106         h = 0;
2107       else
2108         h = P->minr / dot(f[i],f[i]);
2109       push_point(result, scale (h, f[i]));
2110
2111     } else if (*star && P->even != -1) {
2112       /* find the self-intersection of a crossed parallelogram.
2113        * hit is set if v0v1 intersects v2v3*/
2114       Vector v0, v1, v2, v3, c0, c1, p;
2115       double d0, d1;
2116       v0 = v[P->incid[0][i]];
2117       v1 = v[P->incid[1][i]];
2118       v2 = v[P->incid[2][i]];
2119       v3 = v[P->incid[3][i]];
2120       d0 = sqrt(dot(diff(v0, v2), diff(v0, v2)));
2121       d1 = sqrt(dot (diff(v1, v3), diff(v1, v3)));
2122       c0 = scale(d1, sum(v0, v2));
2123       c1 = scale(d0, sum(v1, v3));
2124       p = scale(0.5 / (d0 + d1), sum(c0, c1));
2125       push_point (result, p);
2126       p = cross(diff(p, v2), diff(p, v3));
2127       hit[i] = (dot(p, p) < 1e-6);
2128     } else if (*star && P->hemi && P->index != last_uniform - 1) {
2129       /* find the terminal points of the truncation and the
2130        * self-intersections.
2131        *  v23       v0       v21
2132        *  |  \     /  \     /  |
2133        *  |   v0123    v0321   |
2134        *  |  /     \  /     \  |
2135        *  v01       v2       v03
2136        */
2137       Vector v0, v1, v2, v3, v01, v03, v21, v23, v0123, v0321 ;
2138       Vector u;
2139       double t = 1.5;/* truncation adjustment factor */
2140       j = !P->ftype[P->incid[0][i]];
2141       v0 = v[P->incid[j][i]];/* real vertex */
2142       v1 = v[P->incid[j+1][i]];/* ideal vertex (unit vector) */
2143       v2 = v[P->incid[j+2][i]];/* real */
2144       v3 = v[P->incid[(j+3)%4][i]];/* ideal */
2145       /* compute intersections
2146        * this uses the following linear algebra:
2147        * v0123 = v0 + a v1 = v2 + b v3
2148        * v0 x v3 + a (v1 x v3) = v2 x v3
2149        * a (v1 x v3) = (v2 - v0) x v3
2150        * a (v1 x v3) . (v1 x v3) = (v2 - v0) x v3 . (v1 x v3)
2151        */
2152       u = cross(v1, v3);
2153       v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
2154                             v1));
2155       v0321 = sum(v0, scale(dot(cross(diff(v0, v2), v1), u) / dot(u,u),
2156                             v3));
2157       /* compute truncations */
2158       v01 = sum(v0 , scale(t, diff(v0123, v0)));
2159       v23 = sum(v2 , scale(t, diff(v0123, v2)));
2160       v03 = sum(v0 , scale(t, diff(v0321, v0)));
2161       v21 = sum(v2 , scale(t, diff(v0321, v2)));
2162
2163       push_point(result, v01);
2164       push_point(result, v23);
2165       push_point(result, v0123);
2166       push_point(result, v03);
2167       push_point(result, v21);
2168       push_point(result, v0321);
2169
2170     } else if (*star && P->index == last_uniform - 1) {
2171       /* find the terminal points of the truncation and the
2172        * self-intersections.
2173        *  v23       v0       v21
2174        *  |  \     /  \     /  |
2175        *  |   v0123    v0721   |
2176        *  |  /     \  /     \  |
2177        *  v01       v2       v07
2178        *
2179        *  v65       v4       v67
2180        *  |  \     /  \     /  |
2181        *  |   v4365    v4567   |
2182        *  |  /     \  /     \  |
2183        *  v43       v6       v45
2184        */
2185       Vector v0, v1, v2, v3, v4, v5, v6, v7, v01, v07, v21, v23;
2186       Vector v43, v45, v65, v67, v0123, v0721, v4365, v4567;
2187       double t = 1.5;/* truncation adjustment factor */
2188       Vector u;
2189       for (j = 0; j < 8; j++)
2190         if (P->ftype[P->incid[j][i]] == 3)
2191           break;
2192       v0 = v[P->incid[j][i]];/* real {5/3} */
2193       v1 = v[P->incid[(j+1)%8][i]];/* ideal */
2194       v2 = v[P->incid[(j+2)%8][i]];/* real {3} */
2195       v3 = v[P->incid[(j+3)%8][i]];/* ideal */
2196       v4 = v[P->incid[(j+4)%8][i]];/* real {5/2} */
2197       v5 = v[P->incid[(j+5)%8][i]];/* ideal */
2198       v6 = v[P->incid[(j+6)%8][i]];/* real {3/2} */
2199       v7 = v[P->incid[(j+7)%8][i]];/* ideal */
2200       /* compute intersections */
2201       u = cross(v1, v3);
2202       v0123 = sum(v0, scale(dot(cross(diff(v2, v0), v3), u) / dot(u,u),
2203                             v1));
2204       u = cross(v7, v1);
2205       v0721 = sum(v0, scale(dot(cross(diff(v2, v0), v1), u) / dot(u,u),
2206                             v7));
2207       u = cross(v5, v7);
2208       v4567 = sum(v4, scale(dot(cross(diff(v6, v4), v7), u) / dot(u,u),
2209                             v5));
2210       u = cross(v3, v5);
2211       v4365 = sum(v4, scale(dot(cross(diff(v6, v4), v5), u) / dot(u,u),
2212                             v3));
2213       /* compute truncations */
2214       v01 = sum(v0 , scale(t, diff(v0123, v0)));
2215       v23 = sum(v2 , scale(t, diff(v0123, v2)));
2216       v07 = sum(v0 , scale(t, diff(v0721, v0)));
2217       v21 = sum(v2 , scale(t, diff(v0721, v2)));
2218       v45 = sum(v4 , scale(t, diff(v4567, v4)));
2219       v67 = sum(v6 , scale(t, diff(v4567, v6)));
2220       v43 = sum(v4 , scale(t, diff(v4365, v4)));
2221       v65 = sum(v6 , scale(t, diff(v4365, v6)));
2222
2223       push_point(result, v01);
2224       push_point(result, v23);
2225       push_point(result, v0123);
2226       push_point(result, v07);
2227       push_point(result, v21);
2228       push_point(result, v0721);
2229       push_point(result, v45);
2230       push_point(result, v67);
2231       push_point(result, v4567);
2232       push_point(result, v43);
2233       push_point(result, v65);
2234       push_point(result, v4365);
2235     }
2236
2237   /*
2238    * Face list:
2239    * Each face is printed in a separate line, by listing the indices of its
2240    * vertices. In the non-simple case, the polygon is represented by the
2241    * triangulation, each triangle consists of two polyhedron vertices and one
2242    * auxiliary vertex.
2243    */
2244   Malloc (result->faces, F * 10, face);
2245   result->nfaces = 0;
2246
2247   ii = V;
2248   facelets = 0;
2249   for (i = 0; i < F; i++) {
2250     if (*star) {
2251       if (P->K == 5 &&
2252           (P->D > 30 ||
2253            denominator (P->m[0]) != 1)) {
2254         for (j = 0; j < P->M - 1; j++) {
2255           push_face3 (result, P->incid[j][i], P->incid[j+1][i],  ii);
2256           facelets++;
2257         }
2258
2259         push_face3 (result, P->incid[j][i], P->incid[0][i],  ii++);
2260         facelets++;
2261
2262       } else if (P->even != -1) {
2263         if (hit[i]) {
2264           push_face3 (result, P->incid[3][i], P->incid[0][i],  ii);
2265           push_face3 (result, P->incid[1][i], P->incid[2][i],  ii);
2266         } else {
2267           push_face3 (result, P->incid[0][i], P->incid[1][i],  ii);
2268           push_face3 (result, P->incid[2][i], P->incid[3][i],  ii);
2269         }
2270         ii++;
2271         facelets += 2;
2272
2273       } else if (P->hemi && P->index != last_uniform - 1) {
2274         j = !P->ftype[P->incid[0][i]];
2275
2276         push_face3 (result, ii, ii + 1, ii + 2);
2277         push_face4 (result, P->incid[j][i], ii + 2, P->incid[j+2][i], ii + 5);
2278         push_face3 (result, ii + 3, ii + 4, ii + 5);
2279         ii += 6;
2280         facelets += 3;
2281       } else if (P->index == last_uniform - 1) {
2282         for (j = 0; j < 8; j++)
2283           if (P->ftype[P->incid[j][i]] == 3)
2284             break;
2285         push_face3 (result, ii, ii + 1, ii + 2);
2286         push_face4 (result,
2287                     P->incid[j][i], ii + 2, P->incid[(j+2)%8][i], ii + 5);
2288         push_face3 (result, ii + 3, ii + 4, ii + 5);
2289
2290         push_face3 (result, ii + 6, ii + 7, ii + 8);
2291         push_face4 (result,
2292                     P->incid[(j+4)%8][i], ii + 8, P->incid[(j+6)%8][i],
2293                     ii + 11);
2294         push_face3 (result, ii + 9, ii + 10, ii + 11);
2295         ii += 12;
2296         facelets += 6;
2297       } else {
2298
2299         result->faces[result->nfaces].npoints = P->M;
2300         Malloc (result->faces[result->nfaces].points, P->M, int);
2301         for (j = 0; j < P->M; j++)
2302           result->faces[result->nfaces].points[j] = P->incid[j][i];
2303         result->nfaces++;
2304         facelets++;
2305       }
2306     } else {
2307       int split = (frac(P->n[P->ftype[i]]),
2308                    frax.d != 1 && frax.d != frax.n - 1);
2309       for (j = 0; j < V; j++) {
2310         for (k = 0; k < P->M; k++)
2311           if (P->incid[k][j] == i)
2312             break;
2313         if (k != P->M)
2314           break;
2315       }
2316       if (split) {
2317         ll = j;
2318         for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
2319           for (k = 0; k < P->M; k++)
2320             if (P->incid[k][l] == i)
2321               break;
2322           if (P->adj[k][l] == ll)
2323             k = mod(k + 1 , P->M);
2324           push_face3 (result, ll, l, ii);
2325           facelets++;
2326           ll = l;
2327         }
2328         push_face3 (result, ll, j, ii++);
2329         facelets++;
2330
2331       } else {
2332
2333         int *pp;
2334         int pi = 0;
2335         Malloc (pp, 100, int);
2336
2337         pp[pi++] = j;
2338         ll = j;
2339         for (l = P->adj[k][j]; l != j; l = P->adj[k][l]) {
2340           for (k = 0; k < P->M; k++)
2341             if (P->incid[k][l] == i)
2342               break;
2343           if (P->adj[k][l] == ll)
2344             k = mod(k + 1 , P->M);
2345           pp[pi++] = l;
2346           ll = l;
2347         }
2348         result->faces[result->nfaces].npoints = pi;
2349         result->faces[result->nfaces].points = pp;
2350         result->nfaces++;
2351         facelets++;
2352       }
2353     }
2354   }
2355
2356   /*
2357    * Face color indices - for polyhedra with multiple face types
2358    * For non-simple faces, the index is repeated as many times as needed by the
2359    * triangulation.
2360    */
2361   {
2362     int ff = 0;
2363     if (!*star && P->N != 1) {
2364       for (i = 0; i < F; i++)
2365         if (frac(P->n[P->ftype[i]]), frax.d == 1 || frax.d == frax.n - 1)
2366           result->faces[ff++].color = P->ftype[i];
2367         else
2368           for (j = 0; j < frax.n; j++)
2369             result->faces[ff++].color = P->ftype[i];
2370     } else {
2371       for (i = 0; i < facelets; i++)
2372         result->faces[ff++].color = 0;
2373     }
2374   }
2375
2376   if (*star && P->even != -1)
2377     free(hit);
2378   free(v);
2379   free(f);
2380
2381   return result;
2382 }
2383
2384
2385 \f
2386 /* External interface (jwz)
2387  */
2388
2389 void
2390 free_polyhedron (polyhedron *p)
2391 {
2392   if (!p) return;
2393   Free (p->wythoff);
2394   Free (p->name);
2395   Free (p->group);
2396   Free (p->class);
2397   if (p->faces)
2398     {
2399       int i;
2400       for (i = 0; i < p->nfaces; i++)
2401         Free (p->faces[i].points);
2402       Free (p->faces);
2403     }
2404   Free (p);
2405 }
2406
2407
2408 int
2409 construct_polyhedra (polyhedron ***polyhedra_ret)
2410 {
2411   double freeze = 0;
2412   double azimuth = AZ;
2413   double elevation = EL;
2414   int index = 0;
2415
2416   int count = 0;
2417   polyhedron **result;
2418   Malloc (result, last_uniform * 2 + 1, polyhedron*);
2419
2420   while (index < last_uniform) {
2421     static char sym[4];
2422     Polyhedron *P;
2423
2424     sprintf(sym, "#%d", index + 1);
2425     if (!(P = kaleido(sym, 1, 0, 0, 0))) {
2426       Err (strerror(errno));
2427     }
2428
2429     result[count++] = construct_polyhedron (P, P->v, P->V, P->f, P->F,
2430                                             P->name, P->dual_name,
2431                                             P->class, "",
2432                                             azimuth, elevation, freeze);
2433
2434     result[count++] = construct_polyhedron (P, P->f, P->F, P->v, P->V,
2435                                             P->dual_name, P->name,
2436                                             P->dual_class, "*",
2437                                             azimuth, elevation, freeze);
2438     polyfree(P);
2439     index++;
2440   }
2441
2442   *polyhedra_ret = result;
2443   return count;
2444 }