From http://www.jwz.org/xscreensaver/xscreensaver-5.37.tar.gz
[xscreensaver] / hacks / glx / glhanoi.c
index 6ac0d1bcfad120a64e622316f47f17f441729872..e390f688c0663fa4f72a1f8281fee88b2968b687 100644 (file)
@@ -1,6 +1,7 @@
 /* -*- Mode: C; tab-width: 4 -*- */
-/* glhanoi, Copyright (c) 2005 Dave Atkinson <dave.atkinson@uwe.ac.uk>
+/* glhanoi, Copyright (c) 2005, 2009 Dave Atkinson <da@davea.org.uk>
  * except noise function code Copyright (c) 2002 Ken Perlin
+ * Modified by Lars Huttar (c) 2010, to generalize to 4 or more poles
  *
  * Permission to use, copy, modify, distribute, and sell this software and its
  * documentation for any purpose is hereby granted without fee, provided that
 #define DEF_LIGHT     "True"
 #define DEF_FOG       "False"
 #define DEF_TEXTURE   "True"
+#define DEF_POLES     "0"   /* choose random value */
+#define DEF_SPEED        "1"
+#define DEF_TRAILS       "2"
 
 #define DEFAULTS "*delay:     15000\n" \
                                 "*count:     0\n" \
                                 "*showFPS:   False\n" \
                                 "*wireframe: False\n"
-
+                                
 # define refresh_glhanoi 0
+# define release_glhanoi 0
 
+/* polygon resolution of poles and disks */
 #define NSLICE 32
 #define NLOOPS 1
+
+/* How long to wait at start and finish (seconds). */
 #define START_DURATION 1.0
 #define FINISH_DURATION 1.0
 #define BASE_LENGTH 30.0
 #define BOARD_SQUARES 8
 
+/* Don't draw trail lines till they're this old (sec).
+       Helps trails not be "attached" to the disks. */
+#define TRAIL_START_DELAY 0.1
+
 #define MAX_CAMERA_RADIUS 250.0
 #define MIN_CAMERA_RADIUS 75.0
 
 #define MARBLE_SCALE 1.01
 
 #undef BELLRAND
+/* Return a double precision number in [0...n], with bell curve distribution. */
 #define BELLRAND(n) ((frand((n)) + frand((n)) + frand((n))) / 3)
 
 enum {
@@ -87,11 +100,14 @@ typedef struct {
        GLfloat base0;
        GLfloat base1;
        GLfloat height;
-       GLfloat xmin, xmax, ymin;
+       GLfloat xmin, xmax, ymin, zmin, zmax;
        GLfloat u1, u2;
        GLfloat t1, t2;
        GLfloat ucostheta, usintheta;
-       GLdouble rotAngle;
+       GLfloat dx, dz;
+       GLdouble rotAngle; /* degree of "flipping" so far, during travel */
+       GLdouble phi; /* angle of motion in xz plane */
+       GLfloat speed;
     int polys;
 } Disk;
 
@@ -99,7 +115,22 @@ typedef struct {
        Disk **data;
        int count;
        int size;
-} Stack;
+       GLfloat position[3];
+} Pole;
+
+/* A SubProblem is a recursive subdivision of the problem, and means
+  "Move nDisks disks from src pole to dst pole, using the poles indicated in 'available'." */
+typedef struct {
+       int nDisks;
+       int src, dst;
+       unsigned long available; /* a bitmask of poles that have no smaller disks on them */
+} SubProblem;
+
+typedef struct {
+       GLfloat position[3];
+       double startTime, endTime;
+       Bool isEnd;
+} TrailPoint;
 
 typedef struct {
        GLXContext *glx_context;
@@ -107,22 +138,30 @@ typedef struct {
        Bool wire;
        Bool fog;
        Bool light;
+       Bool layoutLinear;
+       GLfloat trailDuration;
        double startTime;
        double lastTime;
        double duration;
        int numberOfDisks;
+       int numberOfPoles;
        int numberOfMoves;
        int maxDiskIdx;
        int magicNumber;
        Disk *currentDisk;
        int move;
+       /* src, tmp, dst: index of pole that is source / storage / destination for
+               current move */
        int src;
        int tmp;
        int dst;
        int oldsrc;
        int oldtmp;
        int olddst;
-       Stack pole[3];
+       GLfloat speed; /* coefficient for how fast the disks move */
+       SubProblem *solveStack;
+       int solveStackSize, solveStackIdx;
+       Pole *pole;
        float boardSize;
        float baseLength;
        float baseWidth;
@@ -130,13 +169,18 @@ typedef struct {
        float poleRadius;
        float poleHeight;
        float poleOffset;
+       float poleDist; /* distance of poles from center, for round layout */
        float diskHeight;
+       float maxDiskRadius;
        float *diskPos;                         /* pre-computed disk positions on rods */
        Disk *disk;
        GLint floorList;
        GLint baseList;
        GLint poleList;
     int floorpolys, basepolys, polepolys;
+       int trailQSize;
+       TrailPoint *trailQ;
+       int trailQFront, trailQBack;
        GLfloat camera[3];
        GLfloat centre[3];
        rotator *the_rotator;
@@ -146,15 +190,16 @@ typedef struct {
        int drag_x;
        int drag_y;
        int noise_initted;
-
-        int p[512];
-
+       int p[512];
 } glhcfg;
 
 static glhcfg *glhanoi_cfg = NULL;
 static Bool fog;
 static Bool light;
 static Bool texture;
+static GLfloat trails;
+static int poles;
+static GLfloat speed;
 
 static XrmOptionDescRec opts[] = {
        {"-light", ".glhanoi.light", XrmoptionNoArg, "true"},
@@ -162,19 +207,28 @@ static XrmOptionDescRec opts[] = {
        {"-fog", ".glhanoi.fog", XrmoptionNoArg, "true"},
        {"+fog", ".glhanoi.fog", XrmoptionNoArg, "false"},
        {"-texture", ".glhanoi.texture", XrmoptionNoArg, "true"},
-       {"+texture", ".glhanoi.texture", XrmoptionNoArg, "false"}
+       {"+texture", ".glhanoi.texture", XrmoptionNoArg, "false"},
+       {"-trails", ".glhanoi.trails", XrmoptionSepArg, 0},
+       {"-poles", ".glhanoi.poles", XrmoptionSepArg, 0 },
+       {"-speed", ".glhanoi.speed", XrmoptionSepArg, 0 }
 };
 
 static argtype vars[] = {
        {&light, "light", "Light", DEF_LIGHT, t_Bool},
        {&fog, "fog", "Fog", DEF_FOG, t_Bool},
-       {&texture, "texture", "Texture", DEF_TEXTURE, t_Bool}
+       {&texture, "texture", "Texture", DEF_TEXTURE, t_Bool},
+       {&trails, "trails", "Trails", DEF_TRAILS, t_Float},
+       {&poles, "poles", "Poles", DEF_POLES, t_Int},
+       {&speed, "speed", "Speed", DEF_SPEED, t_Float}
 };
 
 static OptionStruct desc[] = {
        {"+/-light", "whether to light the scene"},
        {"+/-fog", "whether to apply fog to the scene"},
-       {"+/-texture", "whether to apply texture to the scene"}
+       {"+/-texture", "whether to apply texture to the scene"},
+       {"-trails t", "how long of disk trails to show (sec.)"},
+       {"-poles r", "number of poles to move disks between"},
+       {"-speed s", "speed multiplier"}
 };
 
 ENTRYPOINT ModeSpecOpt glhanoi_opts = { countof(opts), opts, countof(vars), vars, desc };
@@ -182,7 +236,7 @@ ENTRYPOINT ModeSpecOpt glhanoi_opts = { countof(opts), opts, countof(vars), vars
 #ifdef USE_MODULES
 
 ModStruct glhanoi_description = {
-       "glhanoi", "init_glhanoi", "draw_glhanoi", "release_glhanoi",
+       "glhanoi", "init_glhanoi", "draw_glhanoi", NULL,
        "draw_glhanoi", "init_glhanoi", NULL, &glhanoi_opts,
        1000, 1, 2, 1, 4, 1.0, "",
        "Towers of Hanoi", 0, NULL
@@ -194,7 +248,9 @@ static const GLfloat cBlack[] = { 0.0, 0.0, 0.0, 1.0 };
 static const GLfloat cWhite[] = { 1.0, 1.0, 1.0, 1.0 };
 static const GLfloat poleColor[] = { 0.545, 0.137, 0.137 };
 static const GLfloat baseColor[] = { 0.34, 0.34, 0.48 };
+/* static const GLfloat baseColor[] = { 0.545, 0.137, 0.137 }; */
 static const GLfloat fogcolor[] = { 0.5, 0.5, 0.5 };
+static GLfloat trailColor[] = { 1.0, 1.0, 1.0, 0.5 };
 
 static const float left[] = { 1.0, 0.0, 0.0 };
 static const float up[] = { 0.0, 1.0, 0.0 };
@@ -215,14 +271,24 @@ static const GLfloat spc1[4] = { 1.0, 1.0, 1.0, 1.0 };
 
 static float g = 3.0 * 9.80665;                /* hmm, looks like we need more gravity, Scotty... */
 
+static void checkAllocAndExit(Bool item, char *descr) {
+       if (!item) {
+               fprintf(stderr, "%s: unable to allocate memory for %s\n",
+                               progname, descr);
+               exit(EXIT_FAILURE);
+       }
+}
+
 #define DOPUSH(X, Y)   (((X)->count) >= ((X)->size)) ? NULL : ((X)->data[(X)->count++] = (Y))
 #define DOPOP(X)       (X)->count <= 0 ? NULL : ((X)->data[--((X)->count)])
 
+/* push disk d onto pole idx */
 static Disk *push(glhcfg *glhanoi, int idx, Disk * d)
 {
        return DOPUSH(&glhanoi->pole[idx], d);
 }
 
+/* pop the top disk from pole idx */
 static Disk *pop(glhcfg *glhanoi, int idx)
 {
        return DOPOP(&glhanoi->pole[idx]);
@@ -237,6 +303,7 @@ static inline void swap(int *x, int *y)
 
 /*
  * magic - it's magic...
+ *  Return 1 if the number of trailing zeroes on i is even, unless i is 1 or 0.
  */
 static int magic(int i)
 {
@@ -250,7 +317,6 @@ static int magic(int i)
        return count % 2 == 0;
 }
 
-#if 0
 static float distance(float *p0, float *p1)
 {
        float x, y, z;
@@ -259,13 +325,18 @@ static float distance(float *p0, float *p1)
        z = p1[2] - p0[2];
        return (float)sqrt(x * x + y * y + z * z);
 }
-#endif
 
+/* What is this for?
+  = c / (a b - 0.25 (a^2 + 2 a b + b^2) )
+  = c / (-0.25 (a^2 - 2 a b + b^2) )
+  = c / (-0.25 ((a - b)(a - b)))
+  = -4 c / (a - b)^2
 static GLfloat A(GLfloat a, GLfloat b, GLfloat c)
 {
        GLfloat sum = a + b;
        return c / (a * b - 0.25 * sum * sum);
 }
+*/
 
 static void moveSetup(glhcfg *glhanoi, Disk * disk)
 {
@@ -275,35 +346,55 @@ static void moveSetup(glhcfg *glhanoi, Disk * disk)
        int dst = glhanoi->dst;
        GLfloat theta;
        GLfloat sintheta, costheta;
+       double dh;
+       double dx, dz; /* total x and z distances from src to dst */
+       Pole *poleSrc, *poleDst;
 
-       if(glhanoi->state != FINISHED && random() % 6 == 0) {
-               disk->rotAngle =
-                       -180.0 * (2 - 2 * random() % 2) * (random() % 3 + 1);
+       poleSrc = &(glhanoi->pole[src]);
+       poleDst = &(glhanoi->pole[dst]);
+
+       disk->xmin = poleSrc->position[0];
+       /* glhanoi->poleOffset * (src - (glhanoi->numberOfPoles - 1.0f) * 0.5); */
+       disk->xmax = poleDst->position[0];
+       /* disk->xmax = glhanoi->poleOffset * (dst - (glhanoi->numberOfPoles - 1.0f) * 0.5); */
+       disk->ymin = glhanoi->poleHeight;
+       disk->zmin = poleSrc->position[2];
+       disk->zmax = poleDst->position[2];
+       
+       dx = disk->xmax - disk->xmin;
+       dz = disk->zmax - disk->zmin;
+
+       if(glhanoi->state != FINISHED) {
+               double xxx = ((dx < 0) ? 180.0 : -180.0);
+               if(random() % 6 == 0) {
+                       disk->rotAngle = xxx * (2 - 2 * random() % 2) * (random() % 3 + 1);
+               } else {
+                       disk->rotAngle = xxx;
+               }
+               if(random() % 4 == 0) {
+                       /* Backflip */
+                       disk->rotAngle = -disk->rotAngle;
+               }
        } else {
                disk->rotAngle = -180.0;
        }
 
-       disk->base0 = glhanoi->diskPos[glhanoi->pole[src].count];
-       disk->base1 =
-               glhanoi->state ==
-               FINISHED ? disk->base0 : glhanoi->diskPos[glhanoi->pole[dst].
-                                                                                                 count];
-
-       disk->xmin = glhanoi->poleOffset * (src - 1);
-       disk->xmax = glhanoi->poleOffset * (dst - 1);
-       disk->ymin = glhanoi->poleHeight;
-       ymax =
-               glhanoi->poleHeight + fabs(disk->xmax -
-                                                                  disk->xmin) * (glhanoi->state ==
-                                                                                                 FINISHED ? 1.0 +
-                                                                                                 (double)(glhanoi->
-                                                                                                                  numberOfDisks -
-                                                                                                                  disk->id) /
-                                                                                                 (double)glhanoi->
-                                                                                                 numberOfDisks : 1.0);
-
+       disk->base0 = glhanoi->diskPos[poleSrc->count];
+       disk->base1 = (glhanoi->state == FINISHED) ?
+               disk->base0 : glhanoi->diskPos[poleDst->count];
+
+       /* horizontal distance to travel? */
+       /* was: absx = sqrt(disk->xmax - disk->xmin); */
+       dh = distance(poleSrc->position, poleDst->position);
+       /* absx = sqrt(dh); */
+       ymax = glhanoi->poleHeight + dh;
+       if(glhanoi->state == FINISHED) {
+               ymax += dh * (double)(glhanoi->numberOfDisks - disk->id);
+       }
        h = ymax - disk->ymin;
-       theta = atan((disk->xmin - disk->xmax) * A(disk->xmin, disk->xmax, h));
+       /* A(a, b, c) = -4 c / (a - b)^2 */
+       /* theta = atan(4 h / (b - a)) */
+       theta = atan(4 * h / dh);
        if(theta < 0.0)
                theta += M_PI;
        costheta = cos(theta);
@@ -311,33 +402,177 @@ static void moveSetup(glhcfg *glhanoi, Disk * disk)
        u = (float)
                sqrt(fabs
                         (-g /
-                         (2.0 * A(disk->xmin, disk->xmax, h) * costheta * costheta)));
+                         /* (2.0 * A(disk->xmin, disk->xmax, h) * costheta * costheta))); */
+                         (2.0 * -4 * h / (dh * dh) * costheta * costheta)));
        disk->usintheta = u * sintheta;
        disk->ucostheta = u * costheta;
+       /* Not to be confused: disk->dx is the per-time-unit portion of dx */
+       disk->dx = disk->ucostheta * dx / dh;
+       disk->dz = disk->ucostheta * dz / dh;
        disk->t1 =
                (-u + sqrt(u * u + 2.0 * g * fabs(disk->ymin - disk->base0))) / g;
        disk->u1 = u + g * disk->t1;
        disk->t2 = 2.0 * disk->usintheta / g;
        disk->u2 = disk->usintheta - g * disk->t2;
+
+       /* Compute direction of travel, in the XZ plane. */
+       disk->phi = atan(dz / dx);
+       disk->phi *= 180.0 / M_PI; /* convert radians to degrees */
+}
+
+/* For debugging: show a value as a string of ones and zeroes
+static const char *byteToBinary(int x) {
+    static char b[9];
+       int i, z;
+
+    for (z = 128, i = 0; z > 0; z >>= 1, i++) {
+               b[i] = ((x & z) == z) ? '1' : '0';
+    }
+       b[i] = '\0';
+
+    return b;
+}
+*/
+
+static void pushMove(glhcfg *glhanoi, int n, int src, int dst, int avail) {
+       SubProblem *sp = &(glhanoi->solveStack[glhanoi->solveStackIdx++]);
+       if (glhanoi->solveStackIdx > glhanoi->solveStackSize) {
+               fprintf(stderr, "solveStack overflow: pushed index %d: %d from %d to %d, using %d\n",
+               glhanoi->solveStackIdx, n, src, dst, avail);
+               exit(EXIT_FAILURE);
+       }
+
+       sp->nDisks = n;
+       sp->src = src;
+       sp->dst = dst;
+       sp->available = avail & ~((unsigned long)(1 << src))
+               & ~((unsigned long)(1 << dst));
+       /*
+       fprintf(stderr, "Debug: >  pushed solveStack %d: %d from %d to %d, using %s\n",
+               glhanoi->solveStackIdx - 1, n, src, dst, byteToBinary(sp->available));
+       */
+}
+
+static Bool solveStackEmpty(glhcfg *glhanoi) {
+       return (glhanoi->solveStackIdx < 1);
 }
 
+static SubProblem *popMove(glhcfg *glhanoi) {
+       SubProblem *sp;
+       if (solveStackEmpty(glhanoi)) return (SubProblem *)NULL;
+       sp = &(glhanoi->solveStack[--glhanoi->solveStackIdx]);
+       /* fprintf(stderr, "Debug:  < popped solveStack %d: %d from %d to %d, using %s\n",
+                       glhanoi->solveStackIdx, sp->nDisks, sp->src, sp->dst, byteToBinary(sp->available)); */
+       return sp;
+}
+
+/* Return number of bits set in b */
+static int numBits(unsigned long b) {
+   int count = 0;
+   while (b) {
+      count += b & 0x1u;
+      b >>= 1;
+   }
+   return count;
+}
+
+/* Return index (power of 2) of least significant 1 bit. */
+static int bitScan(unsigned long b) {
+       int count;
+       for (count = 0; b; count++, b >>= 1) {
+               if (b & 0x1u) return count;
+       }
+       return -1;
+}
+
+/* A bit pattern representing all poles */
+#define ALL_POLES ((1 << glhanoi->numberOfPoles) - 1)
+
+#define REMOVE_BIT(a, b)  ((a) & ~(1 << (b)))
+#define    ADD_BIT(a, b)  ((a) |  (1 << (b)))
+
 static void makeMove(glhcfg *glhanoi)
 {
-       int fudge = glhanoi->move + 2;
-       int magicNumber = magic(fudge);
+       if (glhanoi->numberOfPoles == 3) {
+               int fudge = glhanoi->move + 2;
+               int magicNumber = magic(fudge);
 
-       glhanoi->currentDisk = pop(glhanoi, glhanoi->src);
-       moveSetup(glhanoi, glhanoi->currentDisk);
-       push(glhanoi, glhanoi->dst, glhanoi->currentDisk);
-       fudge = fudge % 2;
 
-       if(fudge == 1 || magicNumber) {
-               swap(&glhanoi->src, &glhanoi->tmp);
-       }
-       if(fudge == 0 || glhanoi->magicNumber) {
-               swap(&glhanoi->dst, &glhanoi->tmp);
+               glhanoi->currentDisk = pop(glhanoi, glhanoi->src);
+               moveSetup(glhanoi, glhanoi->currentDisk);
+               push(glhanoi, glhanoi->dst, glhanoi->currentDisk);
+               
+               fudge = fudge % 2;
+
+               if(fudge == 1 || magicNumber) {
+                       swap(&glhanoi->src, &glhanoi->tmp);
+               }
+               if(fudge == 0 || glhanoi->magicNumber) {
+                       swap(&glhanoi->dst, &glhanoi->tmp);
+               }
+               glhanoi->magicNumber = magicNumber;
+       } else {
+               SubProblem sp;
+               int tmp = 0;
+               
+               if (glhanoi->move == 0) {
+                       /* Initialize the solution stack. Original problem:
+                               move all disks from pole 0 to furthest pole,
+                               using all other poles. */
+                       pushMove(glhanoi, glhanoi->numberOfDisks, 0,
+                               glhanoi->numberOfPoles - 1,
+                               REMOVE_BIT(REMOVE_BIT(ALL_POLES, 0), glhanoi->numberOfPoles - 1));
+               }
+               
+               while (!solveStackEmpty(glhanoi)) {
+                       int k, numAvail;
+                       sp = *popMove(glhanoi);
+
+                       if (sp.nDisks == 1) {
+                               /* We have a single, concrete move to do. */
+                               /* moveSetup uses glhanoi->src, dst. */
+                               glhanoi->src = sp.src;
+                               glhanoi->dst = sp.dst;
+                               glhanoi->tmp = tmp; /* Probably unnecessary */
+
+                               glhanoi->currentDisk = pop(glhanoi, sp.src);
+                               moveSetup(glhanoi, glhanoi->currentDisk);
+                               push(glhanoi, sp.dst, glhanoi->currentDisk);
+
+                               return;
+                       } else {
+                               /* Divide and conquer, using Frame-Stewart algorithm, until we get to base case */
+                               if (sp.nDisks == 1) break;
+                               
+                               numAvail = numBits(sp.available);
+                               if (numAvail < 2) k = sp.nDisks - 1;
+                               else if(numAvail >= sp.nDisks - 2) k = 1;
+                               /* heuristic for optimal k: sqrt(2n) (see http://www.cs.wm.edu/~pkstoc/boca.pdf) */
+                               else k = (int)(sqrt(2 * sp.nDisks));
+                               
+                               if (k >= sp.nDisks) k = sp.nDisks - 1;
+                               else if (k < 1) k = 1;
+                               
+                               tmp = bitScan(sp.available);
+                               /* fprintf(stderr, "Debug: k is %d, tmp is %d\n", k, tmp); */
+                               if (tmp == -1) {
+                                       fprintf(stderr, "Error: n > 1 (%d) and no poles available\n",
+                                               sp.nDisks);
+                               }
+                               
+                               /* Push on moves in reverse order, since this is a stack. */
+                               pushMove(glhanoi, k, tmp, sp.dst,
+                                       REMOVE_BIT(ADD_BIT(sp.available, sp.src), tmp));
+                               pushMove(glhanoi, sp.nDisks - k, sp.src, sp.dst,
+                                       REMOVE_BIT(sp.available, tmp));
+                               pushMove(glhanoi, k, sp.src, tmp,
+                                       REMOVE_BIT(ADD_BIT(sp.available, sp.dst), tmp));
+                               
+                               /* Repeat until we've found a move we can make. */
+                       }
+               }
        }
-       glhanoi->magicNumber = magicNumber;
 }
 
 static double lerp(double alpha, double start, double end)
@@ -349,28 +584,67 @@ static void upfunc(GLdouble t, Disk * d)
 {
        d->position[0] = d->xmin;
        d->position[1] = d->base0 + (d->u1 - 0.5 * g * t) * t;
+       d->position[2] = d->zmin;
 
        d->rotation[1] = 0.0;
 }
 
 static void parafunc(GLdouble t, Disk * d)
 {
-       d->position[0] = d->xmin + d->ucostheta * t;
+       /* ##was: d->position[0] = d->xmin + d->ucostheta * t; */
+       d->position[0] = d->xmin + d->dx * t;
+       d->position[2] = d->zmin + d->dz * t;
        d->position[1] = d->ymin + (d->usintheta - 0.5 * g * t) * t;
-
-       d->rotation[1] =
-               d->rotAngle * (d->position[0] - d->xmin) / (d->xmax - d->xmin);
+       
+       d->rotation[1] = d->rotAngle * t / d->t2;
+               /* d->rotAngle * (d->position[0] - d->xmin) / (d->xmax - d->xmin); */
 }
 
 static void downfunc(GLdouble t, Disk * d)
 {
        d->position[0] = d->xmax;
        d->position[1] = d->ymin + (d->u2 - 0.5 * g * t) * t;
-
+       d->position[2] = d->zmax;
+       
        d->rotation[1] = 0.0;
 }
 
-static Bool computePosition(GLfloat t, Disk * d)
+#define normalizeQ(i) ((i) >= glhanoi->trailQSize ? (i) - glhanoi->trailQSize : (i))
+#define normalizeQNeg(i) ((i) < 0 ? (i) + glhanoi->trailQSize : (i))
+
+/* Add trail point at position posn at time t onto back of trail queue.
+       Removes old trails if necessary to make room. */
+static void enQTrail(glhcfg *glhanoi, GLfloat *posn)
+{      
+       if (glhanoi->trailQSize && glhanoi->state != MONEY_SHOT)  {
+               TrailPoint *tp = &(glhanoi->trailQ[glhanoi->trailQBack]);
+               double t = getTime();
+               
+               tp->position[0] = posn[0];
+               tp->position[1] = posn[1] + glhanoi->diskHeight;
+               /* Slight jitter to prevent clashing with other trails */
+               tp->position[2] = posn[2] + (glhanoi->move % 23) * 0.01;
+               tp->startTime = t + TRAIL_START_DELAY;
+               tp->endTime = t + TRAIL_START_DELAY + glhanoi->trailDuration;
+               tp->isEnd = False;
+               
+               /* Update queue back/front indices */
+               glhanoi->trailQBack = normalizeQ(glhanoi->trailQBack + 1);
+               if (glhanoi->trailQBack == glhanoi->trailQFront)
+                       glhanoi->trailQFront = normalizeQ(glhanoi->trailQFront + 1);    
+       }
+}
+
+/* Mark last trailpoint in queue as the end of a trail. */
+/* was: #define endTrail(glh) ((glh)->trailQ[(glh)->trailQBack].isEnd = True) */
+static void endTrail(glhcfg *glhanoi) {
+       if (glhanoi->trailQSize)        
+               glhanoi->trailQ[normalizeQNeg(glhanoi->trailQBack - 1)].isEnd = True;
+}
+
+/* Update disk d's position and rotation based on time t.
+       Returns true iff move is finished. */
+static Bool computePosition(glhcfg *glhanoi, GLfloat t, Disk * d)
 {
        Bool finished = False;
 
@@ -378,11 +652,13 @@ static Bool computePosition(GLfloat t, Disk * d)
                upfunc(t, d);
        } else if(t < d->t1 + d->t2) {
                parafunc(t - d->t1, d);
+               enQTrail(glhanoi, d->position);
        } else {
                downfunc(t - d->t1 - d->t2, d);
                if(d->position[1] <= d->base1) {
                        d->position[1] = d->base1;
                        finished = True;
+                       endTrail(glhanoi);
                }
        }
        return finished;
@@ -400,6 +676,7 @@ static void updateView(glhcfg *glhanoi)
        longitude += glhanoi->camera[0];
        latitude += glhanoi->camera[1];
        radius += glhanoi->camera[2];
+       /* FUTURE: tweak this to be smooth: */
        longitude = longitude - floor(longitude);
        latitude = latitude - floor(latitude);
        radius = radius - floor(radius);
@@ -409,7 +686,7 @@ static void updateView(glhcfg *glhanoi)
        if(radius > 0.5) {
                radius = 1.0 - radius;
        }
-
+       
        b = glhanoi->centre[1];
        c = (MIN_CAMERA_RADIUS +
                 radius * (MAX_CAMERA_RADIUS - MIN_CAMERA_RADIUS));
@@ -432,6 +709,13 @@ static void changeState(glhcfg *glhanoi, State state)
        glhanoi->startTime = getTime();
 }
 
+static Bool finishedHanoi(glhcfg *glhanoi) {
+       /* use different criteria depending on algorithm */
+       return (glhanoi->numberOfPoles == 3 ?
+               glhanoi->move >= glhanoi->numberOfMoves :
+               solveStackEmpty(glhanoi));
+}
+
 static void update_glhanoi(glhcfg *glhanoi)
 {
        double t = getTime() - glhanoi->startTime;
@@ -453,13 +737,14 @@ static void update_glhanoi(glhcfg *glhanoi)
                break;
 
        case MOVE_DISK:
-               if(computePosition(t, glhanoi->currentDisk)) {
+               if(computePosition(glhanoi, t * glhanoi->currentDisk->speed, glhanoi->currentDisk)) {
                        changeState(glhanoi, MOVE_FINISHED);
                }
                break;
 
        case MOVE_FINISHED:
-               if(++glhanoi->move < glhanoi->numberOfMoves) {
+               ++glhanoi->move;
+               if(!finishedHanoi(glhanoi)) {
                        makeMove(glhanoi);
                        changeState(glhanoi, MOVE_DISK);
                } else {
@@ -469,9 +754,8 @@ static void update_glhanoi(glhcfg *glhanoi)
                break;
 
        case FINISHED:
-               while(t < glhanoi->duration) {
+               if (t < glhanoi->duration)
                        break;
-               }
                glhanoi->src = glhanoi->olddst;
                glhanoi->dst = glhanoi->oldsrc;
                for(i = 0; i < glhanoi->numberOfDisks; ++i) {
@@ -496,7 +780,7 @@ static void update_glhanoi(glhcfg *glhanoi)
                                continue;
                        }
 
-                       finished = computePosition(t - delay, &glhanoi->disk[i]);
+                       finished = computePosition(glhanoi, t - delay, &glhanoi->disk[i]);
                        glhanoi->disk[i].rotation[1] = 0.0;
 
                        if(!finished) {
@@ -605,7 +889,6 @@ static int drawTube(GLdouble bottomRadius, GLdouble topRadius,
        GLint lastSlice = nSlice - 1;
        GLfloat radius;
        GLfloat innerRadius;
-       GLfloat maxRadius;
 
        if(bottomThickness > bottomRadius) {
                bottomThickness = bottomRadius;
@@ -619,11 +902,11 @@ static int drawTube(GLdouble bottomRadius, GLdouble topRadius,
        if(topThickness < 0.0) {
                topThickness = 0.0;
        }
-       if(topRadius >= bottomRadius) {
+/*     if(topRadius >= bottomRadius) {
                maxRadius = topRadius;
        } else {
                maxRadius = bottomRadius;
-       }
+       } */
 
        /* bottom */
        y = 0.0;
@@ -745,6 +1028,7 @@ static int drawDisk3D(GLdouble inner_radius, GLdouble outer_radius,
                   outer_radius - inner_radius, height, NSLICE, NLOOPS);
 }
 
+/* used for drawing base */
 static int drawCuboid(GLfloat length, GLfloat width, GLfloat height)
 {
        GLfloat xmin = -length / 2.0f;
@@ -802,6 +1086,137 @@ static int drawCuboid(GLfloat length, GLfloat width, GLfloat height)
     return polys;
 }
 
+/* Set normal vector in xz plane, based on rotation around center. */
+static void setNormalV(glhcfg *glhanoi, GLfloat theta, int y1, int y2, int r1) {
+       if (y1 == y2) /* up/down */
+               glNormal3f(0.0, y1 ? 1.0 : -1.0, 0.0);
+       else if (!r1) /* inward */
+               glNormal3f(-cos(theta), 0.0, -sin(theta));
+       else /* outward */
+               glNormal3f(cos(theta), 0.0, sin(theta));        
+}
+
+/* y1, r1, y2, r2 are indices into y, r, beg, end */
+static int drawBaseStrip(glhcfg *glhanoi, int y1, int r1, int y2, int r2,
+               GLfloat y[2], GLfloat r[2], GLfloat beg[2][2], GLfloat end[2][2]) {
+       int i;
+       GLfloat theta, costh, sinth, x[2], z[2];
+       GLfloat theta1 = (M_PI * 2) / (glhanoi->numberOfPoles + 1);
+       
+       glBegin(GL_QUAD_STRIP);
+       
+       /* beginning edge */
+       glVertex3f(beg[r1][0], y[y1], beg[r1][1]);
+       glVertex3f(beg[r2][0], y[y2], beg[r2][1]);
+       setNormalV(glhanoi, theta1, y1, y2, r1);
+       
+       for (i = 1; i < glhanoi->numberOfPoles; i++) {
+               theta = theta1 * (i + 0.5);
+               costh = cos(theta);
+               sinth = sin(theta);
+               x[0] = costh * r[0];
+               x[1] = costh * r[1];
+               z[0] = sinth * r[0];
+               z[1] = sinth * r[1];
+               
+               glVertex3f(x[r1], y[y1], z[r1]);
+               glVertex3f(x[r2], y[y2], z[r2]);
+               
+               setNormalV(glhanoi, theta1 * (i + 1), y1, y2, r1);
+       }
+
+       /* end edge */
+       glVertex3f(end[r1][0], y[y1], end[r1][1]);
+       glVertex3f(end[r2][0], y[y2], end[r2][1]);
+       setNormalV(glhanoi, glhanoi->numberOfPoles, y1, y2, r1);
+
+       glEnd();
+    return glhanoi->numberOfPoles;
+}
+
+/* Draw base such that poles are distributed around a regular polygon. */
+static int drawRoundBase(glhcfg *glhanoi) {
+       int polys = 0;
+       GLfloat theta, sinth, costh;
+       
+       /* 
+               r[0] = (minimum) inner radius of base at vertices
+               r[1] = (minimum) outer radius of base at vertices
+               y[0] = bottom of base
+               y[1] = top of base */
+       GLfloat r[2], y[2];
+       /* positions of end points: beginning, end.
+               beg[0] is inner corner of beginning of base, beg[1] is outer corner.
+               beg[i][0] is x, [i][1] is z. */
+       GLfloat beg[2][2], end[2][2], begNorm, endNorm;
+       /* ratio of radius at base vertices to ratio at poles */
+       GLfloat longer = 1.0 / cos(M_PI / (glhanoi->numberOfPoles + 1));
+
+       r[0] = (glhanoi->poleDist - glhanoi->maxDiskRadius) * longer;
+       r[1] = (glhanoi->poleDist + glhanoi->maxDiskRadius) * longer;
+       y[0] = 0;
+       y[1] = glhanoi->baseHeight;
+
+       /* compute beg, end. Make the ends square. */
+       theta = M_PI * 2 / (glhanoi->numberOfPoles + 1);
+
+       costh = cos(theta);
+       sinth = sin(theta);
+       beg[0][0] = (glhanoi->poleDist - glhanoi->maxDiskRadius) * costh +
+               glhanoi->maxDiskRadius * sinth;
+       beg[1][0] = (glhanoi->poleDist + glhanoi->maxDiskRadius) * costh +
+               glhanoi->maxDiskRadius * sinth;
+       beg[0][1] = (glhanoi->poleDist - glhanoi->maxDiskRadius) * sinth -
+               glhanoi->maxDiskRadius * costh;
+       beg[1][1] = (glhanoi->poleDist + glhanoi->maxDiskRadius) * sinth -
+               glhanoi->maxDiskRadius * costh;
+       begNorm = theta - M_PI * 0.5;
+       
+       theta = M_PI * 2 * glhanoi->numberOfPoles / (glhanoi->numberOfPoles + 1);
+
+       costh = cos(theta);
+       sinth = sin(theta);
+       end[0][0] = (glhanoi->poleDist - glhanoi->maxDiskRadius) * costh -
+               glhanoi->maxDiskRadius * sinth;
+       end[1][0] = (glhanoi->poleDist + glhanoi->maxDiskRadius) * costh -
+               glhanoi->maxDiskRadius * sinth;
+       end[0][1] = (glhanoi->poleDist - glhanoi->maxDiskRadius) * sinth +
+               glhanoi->maxDiskRadius * costh;
+       end[1][1] = (glhanoi->poleDist + glhanoi->maxDiskRadius) * sinth +
+               glhanoi->maxDiskRadius * costh;
+       endNorm = theta + M_PI * 0.5;
+       
+       /* bottom: never seen
+       polys  = drawBaseStrip(glhanoi, 0, 0, 0, 1, y, r, beg, end); */
+       /* outside edge */
+       polys += drawBaseStrip(glhanoi, 0, 1, 1, 1, y, r, beg, end);
+       /* top */
+       polys += drawBaseStrip(glhanoi, 1, 1, 1, 0, y, r, beg, end);    
+       /* inside edge */
+       polys += drawBaseStrip(glhanoi, 1, 0, 0, 0, y, r, beg, end);
+       
+       /* Draw ends */
+       glBegin(GL_QUADS);
+
+       glVertex3f(beg[0][0], y[1], beg[0][1]);
+       glVertex3f(beg[1][0], y[1], beg[1][1]);
+       glVertex3f(beg[1][0], y[0], beg[1][1]);
+       glVertex3f(beg[0][0], y[0], beg[0][1]);
+       glNormal3f(cos(begNorm), 0, sin(begNorm));
+
+       glVertex3f(end[0][0], y[0], end[0][1]);
+       glVertex3f(end[1][0], y[0], end[1][1]);
+       glVertex3f(end[1][0], y[1], end[1][1]);
+       glVertex3f(end[0][0], y[1], end[0][1]);
+       glNormal3f(cos(endNorm), 0, sin(endNorm));
+       
+       polys += 2;
+       
+       glEnd();
+
+       return polys;
+}
+
 static int drawDisks(glhcfg *glhanoi)
 {
        int i;
@@ -818,7 +1233,12 @@ static int drawDisks(glhcfg *glhanoi)
                glTranslatef(pos[0], pos[1], pos[2]);
                if(rot[1] != 0.0) {
                        glTranslatef(0.0, glhanoi->diskHeight / 2.0, 0.0);
+                       /* rotate around different axis depending on phi */
+                       if (disk->phi != 0.0)
+                               glRotatef(-disk->phi, 0.0, 1.0, 0.0);
                        glRotatef(rot[1], 0.0, 0.0, 1.0);
+                       if (disk->phi != 0.0)
+                               glRotatef(disk->phi, 0.0, 1.0, 0.0);
                        glTranslatef(0.0, -glhanoi->diskHeight / 2.0, 0.0);
                }
                glCallList(disk->displayList);
@@ -831,48 +1251,73 @@ static int drawDisks(glhcfg *glhanoi)
 
 static GLfloat getDiskRadius(glhcfg *glhanoi, int i)
 {
-       return ((GLfloat) i + 3.0) * glhanoi->poleRadius;
+       GLfloat retVal = glhanoi->maxDiskRadius *
+               ((GLfloat) i + 3.0) / (glhanoi->numberOfDisks + 3.0);
+       return retVal;
 }
 
 static void initData(glhcfg *glhanoi)
 {
-       GLfloat maxDiskRadius;
        int i;
+       GLfloat sinPiOverNP;
 
        glhanoi->baseLength = BASE_LENGTH;
-       glhanoi->poleRadius = glhanoi->baseLength /
-               (2.0 * (3 * glhanoi->numberOfDisks + 7.0));
-       maxDiskRadius = getDiskRadius(glhanoi, glhanoi->numberOfDisks);
-       glhanoi->baseWidth = 2.0 * maxDiskRadius;
+       if (glhanoi->layoutLinear) {
+               glhanoi->maxDiskRadius = glhanoi->baseLength /
+                       (2 * 0.95 * glhanoi->numberOfPoles);
+       } else {
+           sinPiOverNP = sin(M_PI / (glhanoi->numberOfPoles + 1));
+               glhanoi->maxDiskRadius = (sinPiOverNP * glhanoi->baseLength * 0.5 * 0.95) / (1 + sinPiOverNP);
+       }
+
+       glhanoi->poleDist = glhanoi->baseLength * 0.5 - glhanoi->maxDiskRadius;
+       glhanoi->poleRadius = glhanoi->maxDiskRadius / (glhanoi->numberOfDisks + 3.0);
+       /* fprintf(stderr, "Debug: baseL = %f, maxDiskR = %f, poleR = %f\n",
+               glhanoi->baseLength, glhanoi->maxDiskRadius, glhanoi->poleRadius); */
+       glhanoi->baseWidth  = 2.0 * glhanoi->maxDiskRadius;
        glhanoi->poleOffset = 2.0 * getDiskRadius(glhanoi, glhanoi->maxDiskIdx);
        glhanoi->diskHeight = 2.0 * glhanoi->poleRadius;
        glhanoi->baseHeight = 2.0 * glhanoi->poleRadius;
        glhanoi->poleHeight = glhanoi->numberOfDisks *
                glhanoi->diskHeight + glhanoi->poleRadius;
+       /* numberOfMoves only applies if numberOfPoles = 3 */
        glhanoi->numberOfMoves = (1 << glhanoi->numberOfDisks) - 1;
+       /* use golden ratio */
        glhanoi->boardSize = glhanoi->baseLength * 0.5 * (1.0 + sqrt(5.0));
 
-       for(i = 0; i < 3; i++) {
-               if((glhanoi->pole[i].data =
-                       calloc(glhanoi->numberOfDisks, sizeof(Disk *))) == NULL) {
-                       fprintf(stderr, "%s: out of memory creating stack %d\n",
-                                       progname, i);
-                       exit(1);
-               }
+       glhanoi->pole = (Pole *)calloc(glhanoi->numberOfPoles, sizeof(Pole));
+       checkAllocAndExit(!!glhanoi->pole, "poles");
+
+       for(i = 0; i < glhanoi->numberOfPoles; i++) {
+               checkAllocAndExit(
+                       !!(glhanoi->pole[i].data = calloc(glhanoi->numberOfDisks, sizeof(Disk *))),
+                       "disk stack");
                glhanoi->pole[i].size = glhanoi->numberOfDisks;
        }
-       if((glhanoi->diskPos =
-               calloc(glhanoi->numberOfDisks, sizeof(double))) == NULL) {
-               fprintf(stderr, "%s: out of memory creating diskPos\n", progname);
-               exit(1);
-       }
-
+       checkAllocAndExit(
+                       !!(glhanoi->diskPos = calloc(glhanoi->numberOfDisks, sizeof(double))),
+                       "diskPos");
+
+       if (glhanoi->trailQSize) {
+               glhanoi->trailQ = (TrailPoint *)calloc(glhanoi->trailQSize, sizeof(TrailPoint));
+               checkAllocAndExit(!!glhanoi->trailQ, "trail queue");
+       } else glhanoi->trailQ = (TrailPoint *)NULL;
+       glhanoi->trailQFront = glhanoi->trailQBack = 0;
+       
        glhanoi->the_rotator = make_rotator(0.1, 0.025, 0, 1, 0.005, False);
+       /* or glhanoi->the_rotator = make_rotator(0.025, 0.025, 0.025, 0.5, 0.005, False); */
        glhanoi->button_down_p = False;
 
        glhanoi->src = glhanoi->oldsrc = 0;
        glhanoi->tmp = glhanoi->oldtmp = 1;
-       glhanoi->dst = glhanoi->olddst = 2;
+       glhanoi->dst = glhanoi->olddst = glhanoi->numberOfPoles - 1;
+       
+       if (glhanoi->numberOfPoles > 3) {
+               glhanoi->solveStackSize = glhanoi->numberOfDisks + 2;
+               glhanoi->solveStack = (SubProblem *)calloc(glhanoi->solveStackSize, sizeof(SubProblem));
+               checkAllocAndExit(!!glhanoi->solveStack, "solving stack");
+               glhanoi->solveStackIdx = 0;
+       }
 }
 
 static void initView(glhcfg *glhanoi)
@@ -944,15 +1389,19 @@ static double improved_noise(glhcfg *glhanoi, double x, double y, double z)
        x -= floor(x);                          /* FIND RELATIVE X,Y,Z */
        y -= floor(y);                          /* OF POINT IN CUBE. */
        z -= floor(z);
-       u = fade(x),                            /* COMPUTE FADE CURVES */
-               v = fade(y),                    /* FOR EACH OF X,Y,Z. */
-               w = fade(z);
-       A = glhanoi->p[X] + Y, AA = glhanoi->p[A] + Z, AB = glhanoi->p[A + 1] + Z,      /* HASH COORDINATES OF */
-               B = glhanoi->p[X + 1] + Y, BA = glhanoi->p[B] + Z, BB = glhanoi->p[B + 1] + Z;  /* THE 8 CUBE CORNERS, */
-       return lerp(w, lerp(v, lerp(u, grad(glhanoi->p[AA], x, y, z),   /* AND ADD */
-                                                               grad(glhanoi->p[BA], x - 1, y, z)),     /* BLENDED */
-                                               lerp(u, grad(glhanoi->p[AB], x, y - 1, z),      /* RESULTS */
-                                                        grad(glhanoi->p[BB], x - 1, y - 1, z))),       /* FROM 8 CORNERS */
+       u  = fade(x),                           /* COMPUTE FADE CURVES */
+       v  = fade(y),                           /* FOR EACH OF X,Y,Z. */
+       w  = fade(z);
+       A  = glhanoi->p[X] + Y;
+       AA = glhanoi->p[A] + Z;
+       AB = glhanoi->p[A + 1] + Z,     /* HASH COORDINATES OF */
+       B  = glhanoi->p[X + 1] + Y;
+       BA = glhanoi->p[B] + Z;
+       BB = glhanoi->p[B + 1] + Z;     /* THE 8 CUBE CORNERS, */
+       return lerp(w, lerp(v, lerp(u, grad(glhanoi->p[AA], x, y, z),/* AND ADD */
+                                                               grad(glhanoi->p[BA], x - 1, y, z)),/* BLENDED */
+                                               lerp(u, grad(glhanoi->p[AB], x, y - 1, z),/* RESULTS */
+                                                        grad(glhanoi->p[BB], x - 1, y - 1, z))),/* FROM 8 CORNERS */
                                lerp(v, lerp(u, grad(glhanoi->p[AA + 1], x, y, z - 1), grad(glhanoi->p[BA + 1], x - 1, y, z - 1)),      /* OF CUBE */
                                         lerp(u, grad(glhanoi->p[AB + 1], x, y - 1, z - 1),
                                                  grad(glhanoi->p[BB + 1], x - 1, y - 1, z - 1))));
@@ -1003,16 +1452,25 @@ static GLubyte *makeTexture(glhcfg *glhanoi, int x_size, int y_size, int z_size,
        return textureData;
 }
 
-static tex_col_t makeMarbleColours(void)
+static void freeTexCols(tex_col_t*p)
 {
-       tex_col_t marbleColours;
+       free(p->colours);
+       free(p);
+}
+
+static tex_col_t *makeMarbleColours(void)
+{
+       tex_col_t *marbleColours;
        int ncols = 2;
 
-       marbleColours.colours = calloc(sizeof(GLuint), ncols);
-       marbleColours.ncols = ncols;
+       marbleColours = malloc(sizeof(tex_col_t));
+       if(marbleColours == NULL) return NULL;
+       marbleColours->colours = calloc(sizeof(GLuint), ncols);
+       if(marbleColours->colours == NULL) return NULL;
+       marbleColours->ncols = ncols;
 
-       marbleColours.colours[0] = 0x3f3f3f3f;
-       marbleColours.colours[1] = 0xffffffff;
+       marbleColours->colours[0] = 0x3f3f3f3f;
+       marbleColours->colours[1] = 0xffffffff;
 
        return marbleColours;
 }
@@ -1076,17 +1534,20 @@ static void setTexture(glhcfg *glhanoi, int n)
        glBindTexture(GL_TEXTURE_2D, glhanoi->textureNames[n]);
 }
 
+/* returns 1 on failure, 0 on success */
 static int makeTextures(glhcfg *glhanoi)
 {
        GLubyte *marbleTexture;
-       tex_col_t marbleColours;
+       tex_col_t *marbleColours;
 
        glGenTextures(N_TEXTURES, glhanoi->textureNames);
 
-       marbleColours = makeMarbleColours();
+       if((marbleColours = makeMarbleColours()) == NULL) {
+               return 1;
+       }
        if((marbleTexture =
                makeTexture(glhanoi, MARBLE_TEXTURE_SIZE, MARBLE_TEXTURE_SIZE, 1,
-                                       makeMarbleTexture, &marbleColours)) == NULL) {
+                                       makeMarbleTexture, marbleColours)) == NULL) {
                return 1;
        }
 
@@ -1099,6 +1560,7 @@ static int makeTextures(glhcfg *glhanoi)
                                 MARBLE_TEXTURE_SIZE, MARBLE_TEXTURE_SIZE, 0,
                                 GL_RGBA, GL_UNSIGNED_BYTE, marbleTexture);
        free(marbleTexture);
+       freeTexCols(marbleColours);
 
        return 0;
 }
@@ -1113,10 +1575,7 @@ static void initFloor(glhcfg *glhanoi)
        float texIncr = 1.0 / BOARD_SQUARES;
 
     glhanoi->floorpolys = 0;
-       if((glhanoi->floorList = glGenLists(1)) == 0) {
-               fprintf(stderr, "can't allocate memory for floor display list\n");
-               exit(EXIT_FAILURE);
-       }
+       checkAllocAndExit(!!(glhanoi->floorList = glGenLists(1)), "floor display list");
        glNewList(glhanoi->floorList, GL_COMPILE);
        x0 = -glhanoi->boardSize / 2.0;
        tx0 = 0.0f;
@@ -1161,38 +1620,56 @@ static void initFloor(glhcfg *glhanoi)
        glEndList();
 }
 
-static void initTowers(glhcfg *glhanoi)
+static void initBase(glhcfg *glhanoi)
 {
-       if((glhanoi->baseList = glGenLists(1)) == 0) {
-               fprintf(stderr, "can't allocate memory for towers display list\n");
-               exit(EXIT_FAILURE);
-       }
+       checkAllocAndExit(!!(glhanoi->baseList = glGenLists(1)), "tower bases display list");
+
        glNewList(glhanoi->baseList, GL_COMPILE);
        setMaterial(baseColor, cWhite, 50);
-    glhanoi->basepolys = drawCuboid(glhanoi->baseLength, glhanoi->baseWidth,
-                                    glhanoi->baseHeight);
+       if (glhanoi->layoutLinear) {
+               glhanoi->basepolys = drawCuboid(glhanoi->baseLength, glhanoi->baseWidth,
+                                                                               glhanoi->baseHeight);
+       } else {
+               glhanoi->basepolys = drawRoundBase(glhanoi);
+       }
        glEndList();
+}
 
+static void initTowers(glhcfg *glhanoi)
+{
+       int i;
+               
+       checkAllocAndExit(!!(glhanoi->poleList = glGenLists(1)), "poles display list\n");
 
-       if((glhanoi->poleList = glGenLists(1)) == 0) {
-               fprintf(stderr, "can't allocate memory for towers display list\n");
-               exit(EXIT_FAILURE);
-       }
        glNewList(glhanoi->poleList, GL_COMPILE);
-       glPushMatrix();
-       glTranslatef(0.0f, glhanoi->baseHeight, 0.0f);
+       /* glTranslatef(-glhanoi->poleOffset * (glhanoi->numberOfPoles - 1.0f) * 0.5f, glhanoi->baseHeight, 0.0f); */
        setMaterial(poleColor, cWhite, 50);
-    glhanoi->polepolys = drawPole(glhanoi->poleRadius, glhanoi->poleHeight);
-       glPushMatrix();
-       glTranslatef(-glhanoi->poleOffset, 0.0, 0.0);
-       glhanoi->polepolys += drawPole(glhanoi->poleRadius, glhanoi->poleHeight);
-       glPopMatrix();
-       glTranslatef(glhanoi->poleOffset, 0.0, 0.0);
-       glhanoi->polepolys += drawPole(glhanoi->poleRadius, glhanoi->poleHeight);
-       glPopMatrix();
+       for (i = 0; i < glhanoi->numberOfPoles; i++) {          
+               GLfloat *p = glhanoi->pole[i].position;
+               GLfloat rad = (M_PI * 2.0 * (i + 1)) / (glhanoi->numberOfPoles + 1);
+               
+               p[1] = glhanoi->baseHeight;
+
+               if (glhanoi->layoutLinear) {
+                       /* Linear: */
+                       p[0] = -glhanoi->poleOffset * ((glhanoi->numberOfPoles - 1) * 0.5f - i);
+                       p[2] = 0.0f;
+               } else {
+                       /* Circular layout: */
+                       p[0] = cos(rad) * glhanoi->poleDist;
+                       p[2] = sin(rad) * glhanoi->poleDist;
+               }
+               
+               glPushMatrix();
+               glTranslatef(p[0], p[1], p[2]);
+               glhanoi->polepolys = drawPole(glhanoi->poleRadius, glhanoi->poleHeight);
+               glPopMatrix();
+
+       }
        glEndList();
 }
 
+/* Parameterized hue based on input 0.0 - 1.0. */
 static double cfunc(double x)
 {
 #define COMP <
@@ -1200,6 +1677,7 @@ static double cfunc(double x)
                return (1.0 / 12.0) / (1.0 / 7.0) * x;
        }
        if(x < 3.0 / 7.0) {
+               /* (7x - 1) / 6 */
                return (1.0 + 1.0 / 6.0) * x - 1.0 / 6.0;
        }
        if(x < 4.0 / 7.0) {
@@ -1214,11 +1692,8 @@ static double cfunc(double x)
 static void initDisks(glhcfg *glhanoi)
 {
        int i;
-       if((glhanoi->disk =
-               (Disk *) calloc(glhanoi->numberOfDisks, sizeof(Disk))) == NULL) {
-               perror("initDisks");
-               exit(EXIT_FAILURE);
-       }
+       glhanoi->disk = (Disk *) calloc(glhanoi->numberOfDisks, sizeof(Disk));
+       checkAllocAndExit(!!glhanoi->disk, "disks");
 
        for(i = glhanoi->maxDiskIdx; i >= 0; i--) {
                GLfloat height = (GLfloat) (glhanoi->maxDiskIdx - i);
@@ -1228,29 +1703,33 @@ static void initDisks(glhcfg *glhanoi)
                Disk *disk = &glhanoi->disk[i];
 
                disk->id = i;
-               disk->position[0] = -glhanoi->poleOffset;
+               disk->position[0] = glhanoi->pole[0].position[0]; /* -glhanoi->poleOffset * (glhanoi->numberOfPoles - 1.0f) * 0.5; */
                disk->position[1] = glhanoi->diskHeight * height;
-               disk->position[2] = 0.0;
+               disk->position[2] = glhanoi->pole[0].position[2];
                disk->rotation[0] = 0.0;
                disk->rotation[1] = 0.0;
                disk->rotation[2] = 0.0;
                disk->polys = 0;
 
+               /* make smaller disks move faster */
+               disk->speed = lerp(((double)glhanoi->numberOfDisks - i) / glhanoi->numberOfDisks,
+                       1.0, glhanoi->speed);
+               /* fprintf(stderr, "disk id: %d, alpha: %0.2f, speed: %0.2f\n", disk->id,
+                       ((double)(glhanoi->maxDiskIdx - i)) / glhanoi->numberOfDisks, disk->speed); */
+
                color[0] = diskColor;
                color[1] = 1.0f;
                color[2] = 1.0f;
                HSVtoRGBv(color, color);
 
-               if((disk->displayList = glGenLists(1)) == 0) {
-                       fprintf(stderr,
-                                       "can't allocate memory for disk %d display list\n", i);
-                       exit(EXIT_FAILURE);
-               }
+               checkAllocAndExit(!!(disk->displayList = glGenLists(1)), "disk display list");
                glNewList(disk->displayList, GL_COMPILE);
                setMaterial(color, cWhite, 100.0);
                disk->polys += drawDisk3D(glhanoi->poleRadius, 
                                   getDiskRadius(glhanoi, i),
                                   glhanoi->diskHeight);
+               /*fprintf(stderr, "Debug: disk %d has radius %f\n", i,
+                       getDiskRadius(glhanoi, i)); */
                glEndList();
        }
        for(i = glhanoi->maxDiskIdx; i >= 0; --i) {
@@ -1295,10 +1774,81 @@ static int drawTowers(glhcfg *glhanoi)
     return glhanoi->basepolys + glhanoi->polepolys;
 }
 
+static int drawTrails1(glhcfg *glhanoi, double t, double thickness, double alpha) {
+       int i, prev = -1, lines = 0;
+       Bool fresh = False;
+       GLfloat trailDurInv = 1.0f / glhanoi->trailDuration;
+
+       glLineWidth(thickness);
+
+       glBegin(GL_LINES);
+       
+       for (i = glhanoi->trailQFront;
+                       i != glhanoi->trailQBack;
+                       i = normalizeQ(i + 1)) {
+               TrailPoint *tqi = &(glhanoi->trailQ[i]);
+               
+               if (!fresh && t > tqi->endTime) {
+                       glhanoi->trailQFront = normalizeQ(i + 1);
+               } else {
+                       if (tqi->startTime > t) break;
+                       /* Found trails that haven't timed out. */
+                       if (!fresh) fresh = True;
+                       if (prev > -1) {
+                               /* Fade to invisible with age */
+                               trailColor[3] = alpha * (tqi->endTime - t) * trailDurInv;
+                               /* Can't use setMaterial(trailColor, cBlack, 0) because our color needs an alpha value. */
+                               glColor4fv(trailColor);
+                               glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, trailColor);
+                               /* FUTURE: to really do this right, trails should be drawn in back-to-front
+                                       order, so that blending is done correctly.
+                                       Currently it looks poor when a faded trail is in front of, or coincident with,
+                                       a bright trail but is drawn first.
+                                       I think for now it's good enough to recommend shorter trails so they
+                                       never/rarely overlap.
+                                       A jitter per trail arc would also mitigate this problem, to a lesser degree. */
+                               glVertex3fv(glhanoi->trailQ[prev].position);
+                               glVertex3fv(glhanoi->trailQ[i].position);
+                               lines++;
+                       }
+                       if (glhanoi->trailQ[i].isEnd)
+                               prev = -1;
+                       else
+                               prev = i;
+               }               
+       }
+       
+       glEnd();
+       
+       return lines;
+}
+
+static int drawTrails(glhcfg *glhanoi) {
+       int lines = 0;
+       double t = getTime();
+
+       glEnable (GL_BLEND);
+       glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+       glMaterialfv(GL_FRONT, GL_SPECULAR, cBlack);
+       glMateriali(GL_FRONT, GL_SHININESS, 0);
+
+       /* Draw them twice, with different widths and opacities, to make them smoother. */
+       lines = drawTrails1(glhanoi, t, 1.0, 0.75);
+       lines += drawTrails1(glhanoi, t, 2.5, 0.5);
+
+       glDisable (GL_BLEND);
+
+       /* fprintf(stderr, "Drew trails: %d lines\n", lines); */
+       return lines;
+}
+
 /* Window management, etc
  */
 ENTRYPOINT void reshape_glhanoi(ModeInfo * mi, int width, int height)
 {
+       glhcfg *glhanoi = &glhanoi_cfg[MI_SCREEN(mi)];
+       glXMakeCurrent(MI_DISPLAY(mi), MI_WINDOW(mi), *(glhanoi->glx_context));
+
        glViewport(0, 0, (GLint) width, (GLint) height);
 
        glMatrixMode(GL_PROJECTION);
@@ -1312,37 +1862,51 @@ ENTRYPOINT void reshape_glhanoi(ModeInfo * mi, int width, int height)
        glClear(GL_COLOR_BUFFER_BIT);
 }
 
+static void free_glhanoi(ModeInfo * mi);
+
 ENTRYPOINT void init_glhanoi(ModeInfo * mi)
 {
        glhcfg *glhanoi;
-       if(!glhanoi_cfg) {
-               glhanoi_cfg =
-                       (glhcfg *) calloc(MI_NUM_SCREENS(mi), sizeof(glhcfg));
-               if(!glhanoi_cfg) {
-                       fprintf(stderr, "%s: out of memory creating configs\n",
-                                       progname);
-                       exit(1);
-               }
-       }
+       MI_INIT(mi, glhanoi_cfg, free_glhanoi);
 
        glhanoi = &glhanoi_cfg[MI_SCREEN(mi)];
        glhanoi->glx_context = init_GL(mi);
        glhanoi->numberOfDisks = MI_BATCHCOUNT(mi);
-
+       
     if (glhanoi->numberOfDisks <= 1)
       glhanoi->numberOfDisks = 3 + (int) BELLRAND(9);
 
-    /* magicnumber is a bitfield, so we can't have more than 31 discs
-       on a system with 4-byte ints. */
-    if (glhanoi->numberOfDisks >= 8 * sizeof(int))
-      glhanoi->numberOfDisks = (8 * sizeof(int)) - 1;
+       /* magicnumber is a bitfield, so we can't have more than 31 discs
+          on a system with 4-byte ints. */
+       if (glhanoi->numberOfDisks >= 8 * sizeof(int))
+               glhanoi->numberOfDisks = (8 * sizeof(int)) - 1;
 
        glhanoi->maxDiskIdx = glhanoi->numberOfDisks - 1;
+
+       glhanoi->numberOfPoles = get_integer_resource(MI_DISPLAY(mi), "poles", "Integer");
+       /* Set a number of poles from 3 to numberOfDisks + 1, biased toward lower values,
+               with probability decreasing linearly. */
+    if (glhanoi->numberOfPoles <= 2)
+      glhanoi->numberOfPoles = 3 +
+               (int)((1 - sqrt(frand(1.0))) * (glhanoi->numberOfDisks - 1));
+       
        glhanoi->wire = MI_IS_WIREFRAME(mi);
+
+# ifdef HAVE_JWZGLES /* #### glPolygonMode other than GL_FILL unimplemented */
+    glhanoi->wire = 0;
+# endif
+
        glhanoi->light = light;
        glhanoi->fog = fog;
        glhanoi->texture = texture;
-
+       glhanoi->speed = speed;
+       glhanoi->trailDuration = trails;
+       /* set trailQSize based on 60 fps (a maximum, more or less) */
+       /* FUTURE: Should clamp framerate to 60 fps? See flurry.c's draw_flurry().
+               The only bad effect if we don't is that trail-ends could
+               show "unnatural" pauses at high fps. */
+       glhanoi->trailQSize = (int)(trails * 60.0);
+       
        reshape_glhanoi(mi, MI_WIDTH(mi), MI_HEIGHT(mi));
 
        if(glhanoi->wire) {
@@ -1352,14 +1916,15 @@ ENTRYPOINT void init_glhanoi(ModeInfo * mi)
        }
 
        initLights(!glhanoi->wire && glhanoi->light);
-       if(makeTextures(glhanoi) != 0) {
-               fprintf(stderr, "can't allocate memory for marble texture\n");
-               exit(EXIT_FAILURE);
-       }
+       checkAllocAndExit(!makeTextures(glhanoi), "textures\n");
 
+       /* Choose linear or circular layout. Could make this a user option. */
+       glhanoi->layoutLinear = (glhanoi->numberOfPoles == 3);
+       
        initData(glhanoi);
        initView(glhanoi);
        initFloor(glhanoi);
+       initBase(glhanoi);
        initTowers(glhanoi);
        initDisks(glhanoi);
 
@@ -1378,21 +1943,20 @@ ENTRYPOINT void init_glhanoi(ModeInfo * mi)
                glEnable(GL_FOG);
        }
 
-
        glhanoi->duration = START_DURATION;
        changeState(glhanoi, START);
 }
 
 ENTRYPOINT void draw_glhanoi(ModeInfo * mi)
 {
-        glhcfg *glhanoi = &glhanoi_cfg[MI_SCREEN(mi)];
+       glhcfg *glhanoi = &glhanoi_cfg[MI_SCREEN(mi)];
        Display *dpy = MI_DISPLAY(mi);
        Window window = MI_WINDOW(mi);
 
        if(!glhanoi->glx_context)
                return;
 
-        glXMakeCurrent(MI_DISPLAY(mi), MI_WINDOW(mi), *(glhanoi->glx_context));
+    glXMakeCurrent(MI_DISPLAY(mi), MI_WINDOW(mi), *(glhanoi->glx_context));
 
        glPolygonMode(GL_FRONT, glhanoi->wire ? GL_LINE : GL_FILL);
 
@@ -1400,10 +1964,20 @@ ENTRYPOINT void draw_glhanoi(ModeInfo * mi)
     mi->polygon_count = 0;
 
        glLoadIdentity();
+    glRotatef(current_device_rotation(), 0, 0, 1);
 
        update_glhanoi(glhanoi);
        updateView(glhanoi);
 
+# ifdef HAVE_MOBILE    /* Keep it the same relative size when rotated. */
+    {
+      GLfloat h = MI_HEIGHT(mi) / (GLfloat) MI_WIDTH(mi);
+      int o = (int) current_device_rotation();
+      if (o != 0 && o != 180 && o != -180)
+        glScalef (1/h, 1/h, 1/h);
+    }
+# endif
+
        if(!glhanoi->wire && glhanoi->texture) {
                glEnable(GL_TEXTURE_2D);
        }
@@ -1413,6 +1987,11 @@ ENTRYPOINT void draw_glhanoi(ModeInfo * mi)
        mi->polygon_count += drawTowers(glhanoi);
        mi->polygon_count += drawDisks(glhanoi);
 
+       if (glhanoi->trailQSize) {
+               /* No polygons, just lines. So ignore the return count. */
+               (void)drawTrails(glhanoi);
+       }
+       
        if(mi->fps_p) {
                do_fps(mi);
        }
@@ -1425,6 +2004,8 @@ ENTRYPOINT Bool glhanoi_handle_event(ModeInfo * mi, XEvent * event)
 {
        glhcfg *glhanoi = &glhanoi_cfg[MI_SCREEN(mi)];
 
+    /* #### this is all wrong on iOS -- should be using gltrackball. */
+
        if(event->xany.type == ButtonPress && event->xbutton.button == Button1) {
                glhanoi->button_down_p = True;
                glhanoi->drag_x = event->xbutton.x;
@@ -1461,34 +2042,37 @@ ENTRYPOINT Bool glhanoi_handle_event(ModeInfo * mi, XEvent * event)
 
                return True;
        }
+#if 0 /* #### doesn't work */
+  else if (screenhack_event_helper (MI_DISPLAY(mi), MI_WINDOW(mi), event))
+    {
+      changeState(glhanoi, START);
+      return True;
+    }
+#endif
        return False;
 }
 
-ENTRYPOINT void release_glhanoi(ModeInfo * mi)
+static void free_glhanoi(ModeInfo * mi)
 {
-       if(glhanoi_cfg != NULL) {
-               int screen;
-               for(screen = 0; screen < MI_NUM_SCREENS(mi); screen++) {
-                       int i;
-                       int j;
-                       glhcfg *glh = &glhanoi_cfg[screen];
-                       glDeleteLists(glh->floorList, 1);
-                       glDeleteLists(glh->baseList, 1);
-                       glDeleteLists(glh->poleList, 1);
-                        glDeleteLists(glh->textureNames[0], 2);
-                       for(j = 0; j < glh->numberOfDisks; ++j) {
-                               glDeleteLists(glh->disk[j].displayList, 1);
-                       }
-                       free(glh->disk);
-                       for(i = 0; i < 3; i++) {
-                               if(glh->pole[i].data != NULL) {
-                                       free(glh->pole[i].data);
-                               }
+       int i;
+       int j;
+       glhcfg *glh = &glhanoi_cfg[MI_SCREEN(mi)];
+       if (glh->glx_context) {
+               glXMakeCurrent(MI_DISPLAY(mi), MI_WINDOW(mi), *(glh->glx_context));
+               glDeleteLists(glh->floorList, 1);
+               glDeleteLists(glh->baseList, 1);
+               glDeleteLists(glh->poleList, 1);
+                glDeleteLists(glh->textureNames[0], 2);
+               for(j = 0; j < glh->numberOfDisks; ++j) {
+                       glDeleteLists(glh->disk[j].displayList, 1);
+               }
+               free(glh->disk);
+               for(i = 0; i < glh->numberOfPoles; i++) {
+                       if(glh->pole[i].data != NULL) {
+                               free(glh->pole[i].data);
                        }
                }
        }
-       free(glhanoi_cfg);
-       glhanoi_cfg = NULL;
 }
 
 XSCREENSAVER_MODULE ("GLHanoi", glhanoi)