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