/* Copyright (c) Mark J. Kilgard, 1999. */

/* This program is freely distributable without licensing fees  and is
   provided without guarantee or warrantee expressed or  implied. This
   program is -not- in the public domain. */

/* nvshadow.c - this example demonstrates stenciled shadow volume shadows,
   shadow volume visualization, soft shadows through stenciled shadow volumes,
   planar projected shadows, and a hybrid combination of shadow volumes and
   planar projected shadows.  Graphics hardware with a hardware stencil is
   DEFINITELY NEEDED for good performance.  Try it in 32-bit mode on your
   RIVA TNT! */

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <GL/glut.h>

/* Some <math.h> files do not define M_PI... */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef _WIN32
#define CALLBACK  /* compensate for Microsoft's stupidity */
#endif

/* This "fastTeapot" routine is faster than glutSolidTeapot. */
static void fastTeapot(GLint grid, GLdouble scale);

/* 2D points for the complex polygons making up the NVIDIA logo. */
GLdouble nvlogo0[][3] = {
  { -0.474465, -1.259490, 0 },
  { 0.115919, -1.113297, 0 },
  { 0.588227, -0.899634, 0 },
  { 0.942455, -0.652235, 0 },
  { 1.296687, -0.348609, 0 },
  { 1.690275, -0.033738, 0 },
  { 1.926431, 0.269888, 0 },
  { 2.123226, 0.494796, 0 },
  { 1.847713, 0.832160, 0 },
  { 1.532842, 1.124540, 0 },
  { 1.178611, 1.394430, 0 },
  { 0.706303, 1.709300, 0 },
  { 0.076562, 1.967940, 0 },
  { -0.395747, 2.080400, 0 },
  { -1.064784, 2.058940, 0 },
  { -1.064847, 2.811350, 0 },
  { 3.973113, 2.811350, 0 },
  { 3.973113, -2.811350, 0 },
  { -1.025490, -2.811350, 0 },
  { -1.025490, -2.159120, 0 },
  { -0.474465, -2.159120, 0 },
  { 0.155277, -2.102890, 0 },
  { 0.706303, -1.979190, 0 },
  { 1.178611, -1.833000, 0 },
  { 1.690275, -1.653080, 0 },
  { 2.201941, -1.450660, 0 },
  { 2.674248, -1.214510, 0 },
  { 3.107212, -0.955861, 0 },
  { 3.343357, -0.719707, 0 },
  { 2.438097, -0.179928, 0 },
  { 2.005147, -0.517290, 0 },
  { 1.690275, -0.820916, 0 },
  { 1.296687, -1.079560, 0 },
  { 0.863740, -1.338200, 0 },
  { 0.273356, -1.585600, 0 },
  { -0.198952, -1.709300, 0 },
  { -1.025490, -1.731790, 0 },
  { -1.025490, -1.248240, 0 },
};
GLdouble nvlogo1[][3] = {
  { -0.493508, 0.560265, 0 },
  { -0.233835, 0.218981, 0 },
  { -0.078033, -0.107463, 0 },
  { 0.545180, 0.441557, 0 },
  { 0.285509, 0.753164, 0 },
  { -0.129966, 1.005420, 0 },
  { -0.545442, 1.153800, 0 },
  { -1.034999, 1.167860, 0 },
  { -1.064784, 1.658310, 0 },
  { -0.233835, 1.598950, 0 },
  { 0.233576, 1.361540, 0 },
  { 0.649050, 1.094450, 0 },
  { 1.012591, 0.753164, 0 },
  { 1.324197, 0.426719, 0 },
  { 1.064524, 0.189305, 0 },
  { 0.804852, -0.166817, 0 },
  { 0.389378, -0.508100, 0 },
  { -0.078033, -0.745515, 0 },
  { -0.441573, -0.879060, 0 },
  { -1.013530, -0.889070, 0 },
  { -1.012851, 0.723487, 0 },
};
GLdouble nvlogo2[][3] = {
  { -1.025490, -2.159120, 0 },
  { -1.843800, -1.962260, 0 },
  { -2.415081, -1.635820, 0 },
  { -2.934425, -1.205510, 0 },
  { -3.297966, -0.760353, 0 },
  { -3.609571, -0.315201, 0 },
  { -3.869244, 0.204143, 0 },
  { -3.973113, 0.545426, 0 },
  { -3.505702, 0.960900, 0 },
  { -2.830556, 1.435730, 0 },
  { -2.051539, 1.851210, 0 },
  { -1.064784, 2.058940, 0 },
  { -1.064784, 1.658310, 0 },
  { -1.791868, 1.495080, 0 },
  { -2.363145, 1.183480, 0 },
  { -2.830556, 0.842190, 0 },
  { -3.194097, 0.471234, 0 },
  { -3.090228, 0.055759, 0 },
  { -2.830556, -0.315201, 0 },
  { -2.570884, -0.760353, 0 },
  { -2.103473, -1.220340, 0 },
  { -1.584129, -1.531950, 0 },
  { -1.025490, -1.731790, 0 },
};
GLdouble nvlogo3[][3] = {
  { -1.025490, -1.248240, 0 },
  { -1.472016, -1.099371, 0 },
  { -1.794030, -0.875934, 0 },
  { -2.047038, -0.606495, 0 },
  { -2.254046, -0.337056, 0 },
  { -2.415053, -0.047902, 0 },
  { -2.530060, 0.260968, 0 },
  { -2.392054, 0.536978, 0 },
  { -2.047038, 0.806418, 0 },
  { -1.633023, 1.016710, 0 },
  { -1.034999, 1.167860, 0 },
  { -1.012851, 0.723487, 0 },
  { -1.380012, 0.681555, 0 },
  { -1.610022, 0.530407, 0 },
  { -1.863033, 0.326685, 0 },
  { -1.909033, 0.076960, 0 },
  { -1.748027, -0.159620, 0 },
  { -1.541019, -0.448774, 0 },
  { -1.311010, -0.685355, 0 },
  { -1.013530, -0.889070, 0 },
};
GLdouble base[][3] = {
  { -5, -5, 0 },
  { -5,  5, 0 },
  {  5,  5, 0 },
  {  5, -5, 0 }
};

#define SIZE(a) (sizeof(a)/sizeof(a[0]))

enum {
  DL_BOGUS = 0,
  DL_CUT_OUT,
  DL_CUT_OUT_VOLUME,
  DL_TEAPOT,
  DL_SPHERE
};

enum {
  M_ANIMATE,
  M_VOLUME_SHADOWS,
  M_SOFT_SHADOWS,
  M_NO_SHADOWS,
  M_SHADOW_VOLUME,
  M_PLANAR_SHADOWS,
  M_HYBRID_SHADOWS,
  M_BENCHMARK
};

enum {
  X, Y, Z, W  /* Used to make code involving coordinates more readable. */
};

int doubleBuffer = 1;
int useTextures = 1;
int forceStencilHack = 0;
int fullscreen = 0, smallFullscreen = 0;
int moving = 0, beginx, beginy;
int windowWidth, windowHeight;

int teapotAngle = 0;
float time1 = 0.0;
float time2 = 0.0;
int animate = 0;
int animateMode = 1;
int logoAngle = 0;
int oneFrame = 0;
int displayMode = M_HYBRID_SHADOWS;
int numLightSamples = 4;
int useMipmaps = 1;
int linearFiltering = 1;

GLUtesselator *tess;

double L[4] = { 3,4,2, 1 };   /* Light location. */
double O[3] = { 3,2,2 };      /* Cut-out origin. */
double M[3], N[3];            /* Vectors to orient the cut-out. */
double Pbg[4];                /* Plane below the true ground plane. */
double Pg[4];                 /* Plane at the actual ground plane. */

/* Three vertices just below the ground plane. */
double Sbg[3] = { 0, -1.5, 0 };
double Tbg[3] = { 0, -1.5, 1 };
double Rbg[3] = { 1, -1.5, 0 };

/* Three vertices on the ground plane. */
double Sg[3] = { 0, -1, 0 };
double Tg[3] = { 0, -1, 1 };
double Rg[3] = { 1, -1, 0 };

/* Nice floor texture tiling pattern. */
static char *circles[] = {
  "....xxxx........",
  "..xxxxxxxx......",
  ".xxxxxxxxxx.....",
  ".xxx....xxx.....",
  "xxx......xxx....",
  "xxx......xxx....",
  "xxx......xxx....",
  "xxx......xxx....",
  ".xxx....xxx.....",
  ".xxxxxxxxxx.....",
  "..xxxxxxxx......",
  "....xxxx........",
  "................",
  "................",
  "................",
  "................",
};

static void
makeFloorTexture(void)
{
  GLubyte floorTexture[16][16][3];
  GLubyte *loc;
  int s, t;

  /* Setup RGB image for the texture. */
  loc = (GLubyte*) floorTexture;
  for (t = 0; t < 16; t++) {
    for (s = 0; s < 16; s++) {
      if (circles[t][s] == 'x') {
	/* Nice green. */
        loc[0] = 0x1f;
        loc[1] = 0x8f;
        loc[2] = 0x1f;
      } else {
	/* Light gray. */
        loc[0] = 0xaa;
        loc[1] = 0xaa;
        loc[2] = 0xaa;
      }
      loc += 3;
    }
  }

  glPixelStorei(GL_UNPACK_ALIGNMENT, 1);

  if (useMipmaps) {
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER,
      GL_LINEAR_MIPMAP_LINEAR);
    gluBuild2DMipmaps(GL_TEXTURE_2D, 3, 16, 16,
      GL_RGB, GL_UNSIGNED_BYTE, floorTexture);
  } else {
    if (linearFiltering) {
      glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
    } else {
      glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    }
    glTexImage2D(GL_TEXTURE_2D, 0, 3, 16, 16, 0,
      GL_RGB, GL_UNSIGNED_BYTE, floorTexture);
  }
}

/* Create a matrix that will project the desired shadow. */
void
shadowMatrix(GLfloat shadowMat[4][4],
  GLdouble groundplane[4],
  GLfloat lightpos[4])
{
  GLfloat dot;

  /* Find dot product between light position vector and ground plane normal. */
  dot = groundplane[X] * lightpos[X] +
    groundplane[Y] * lightpos[Y] +
    groundplane[Z] * lightpos[Z] +
    groundplane[W] * lightpos[W];

  shadowMat[0][0] = dot - lightpos[X] * groundplane[X];
  shadowMat[1][0] = 0.f - lightpos[X] * groundplane[Y];
  shadowMat[2][0] = 0.f - lightpos[X] * groundplane[Z];
  shadowMat[3][0] = 0.f - lightpos[X] * groundplane[W];

  shadowMat[X][1] = 0.f - lightpos[Y] * groundplane[X];
  shadowMat[1][1] = dot - lightpos[Y] * groundplane[Y];
  shadowMat[2][1] = 0.f - lightpos[Y] * groundplane[Z];
  shadowMat[3][1] = 0.f - lightpos[Y] * groundplane[W];

  shadowMat[X][2] = 0.f - lightpos[Z] * groundplane[X];
  shadowMat[1][2] = 0.f - lightpos[Z] * groundplane[Y];
  shadowMat[2][2] = dot - lightpos[Z] * groundplane[Z];
  shadowMat[3][2] = 0.f - lightpos[Z] * groundplane[W];

  shadowMat[X][3] = 0.f - lightpos[W] * groundplane[X];
  shadowMat[1][3] = 0.f - lightpos[W] * groundplane[Y];
  shadowMat[2][3] = 0.f - lightpos[W] * groundplane[Z];
  shadowMat[3][3] = dot - lightpos[W] * groundplane[W];
}

/* x86 generates tighter double-precision code. */
typedef GLdouble FPvalue;
#define constructShadowVolumeMatrix constructShadowVolumeMatrixd

#if ( _MSC_VER >= 800 )
  /* Assume no aliasing in the constructShadowVolumeMatrix. */
# pragma optimize("a", on)
  /* Visual C++ 5.0 does poor common subexpression folding.
     Define to ENABLE common expression folding. */
# define CEF
#endif

/* PERFORMANCE RESULTS on 350 Mhz Pentium II
     "/Ox -DCEF" -- 160,500 calls/sec
     "/O1 -DCEF" -- 124,000 calls/sec
     "/Ox      " -- 109,500 calls/sec
     "/O1      " --  76,000 calls/sec
*/

/* Construct a shadow volume matrix.  The code below was generated by Maple.
   The code implements an inlined 15x15 matrix inverse for the special case
   we are interested in.  No guarantees if this code is numerically stable,
   but it is good enough for nvshadow and pretty fast.

   Consult the "Reconstruction" chapter of Michael Penna and Richard
   Patterson's book _Projective Geometry and its Applications to Computer
   Graphics_ (Prentice-Hall, 1986) for more details. */
void
constructShadowVolumeMatrix(FPvalue P[5][3], FPvalue M[4][4])
{
#ifndef CEF
  /* Common sub-expressions NOT folded together via temporary variables (except for "r"). */
  FPvalue r;
  const FPvalue two = (FPvalue) 2.0;
  const FPvalue three = (FPvalue) 3.0;

  assert( (((char*)M) < ((char*)P)) || (((char*)M) >= (((char*)P)+sizeof(P))) );
  assert( (((char*)P) < ((char*)M)) || (((char*)P) >= (((char*)M)+sizeof(M))) );

  r = 1/(-P[4][2]*P[1][1]*P[2][0]+P[2][0]*P[1][1]*P[3][2]+P[4][1]*P[1][2]*P[2][0]-P[2][0]*
    P[3][1]*P[1][2]-P[2][0]*P[4][1]*P[3][2]+P[2][0]*P[3][1]*P[4][2]+P[2][2]*P[1][1]*P[4][0]+P[1][1]*P[4][2]*P[3][0]-P[1][1]*P[2][2]*P[3][0]-P[1][1]*P[3][2]*P[4][0]+P[2][1]*P[4][0]*P[3][2]+P[4][1]*P[1][0]*P[3][2]+
    P[2][1]*P[3][0]*P[1][2]+P[4][2]*P[2][1]*P[1][0]+P[3][1]*P[1][2]*P[4][0]+P[4][1]*P[2][2]*P[3][0]-P[2][2]*P[4][1]*P[1][0]-P[4][1]*P[3][0]*P[1][2]-P[3][1]*P[2][2]*P[4][0]+P[3][1]*P[2][2]*P[1][0]-P[2][1]*P[1][0]*
    P[3][2]-P[3][1]*P[4][2]*P[1][0]-P[2][1]*P[3][0]*P[4][2]-P[2][1]*P[1][2]*P[4][0]);

  M[0][0] = (-two*P[0][2]*P[2][0]*P[1][1]*P[3][0]+two*P[0][2]*P[2][0]*P[1][1]*P[4][0]-two*P[0][2]*P[2][0]*P[4][1]*P[1][0]-two*P[0][2]*P[2][0]*P[3][1]*
    P[4][0]+two*P[0][2]*P[2][0]*P[3][1]*P[1][0]+two*P[0][2]*P[2][0]*P[4][1]*P[3][0]-two*P[0][1]*P[2][0]*P[1][2]*P[4][0]+two*P[0][1]*P[2][0]*P[1][2]*P[3][0]-two*P[0][1]*P[2][0]*P[3][2]
    *P[1][0]+two*P[0][1]*P[2][0]*P[4][2]*P[1][0]+two*P[0][1]*P[2][0]*P[3][2]*P[4][0]-two*P[0][1]*P[2][0]*P[4][2]*P[3][0]+two*P[2][0]*P[3][1]*P[1][2]*P[4][0]-two*P[2][0]*P[4][1]*
    P[3][0]*P[1][2]-two*P[2][0]*P[1][1]*P[3][2]*P[4][0]+two*P[2][0]*P[1][1]*P[4][2]*P[3][0]+P[0][0]*P[2][1]*P[1][0]*P[3][2]+P[0][0]*P[3][1]*P[4][2]*P[1][0]+P[0][0]*P[2][1]*P[3][0]*P[4][2]+P[0][0]*P[2][1]
    *P[1][2]*P[4][0]+two*P[1][0]*P[2][0]*P[4][1]*P[3][2]-two*P[1][0]*P[2][0]*P[3][1]*P[4][2]+P[0][0]*P[4][1]*P[3][0]*P[1][2]+P[0][0]*P[3][1]*P[2][2]*P[4][0]-P[0][0]*P[3][1]*P[2][2]*P[1][0]-P[0][0]*
    P[2][0]*P[3][1]*P[1][2]-P[0][0]*P[2][0]*P[4][1]*P[3][2]+P[0][0]*P[2][0]*P[3][1]*P[4][2]-P[0][0]*P[2][2]*P[1][1]*P[4][0]-P[0][0]*P[1][1]*P[4][2]*P[3][0]+P[0][0]*P[1][1]*P[2][2]*P[3][0]+P[0][0]*P[1][1]*P[3][2]*
    P[4][0]-P[0][0]*P[2][1]*P[4][0]*P[3][2]-P[0][0]*P[4][1]*P[1][0]*P[3][2]-P[0][0]*P[2][1]*P[3][0]*P[1][2]-P[0][0]*P[4][2]*P[2][1]*P[1][0]-P[0][0]*P[3][1]*P[1][2]*P[4][0]-P[0][0]*P[4][1]*P[2][2]*P[3][0]+P[0][0]*
    P[2][2]*P[4][1]*P[1][0]-P[0][0]*P[4][2]*P[1][1]*P[2][0]+P[0][0]*P[2][0]*P[1][1]*P[3][2]+P[0][0]*P[4][1]*P[1][2]*P[2][0])*r;
  M[0][1] = -(-two*P[2][1]*P[1][0]*P[4][1]*P[3][2]+P[0][1]*P[2][1]*P[1][0]*P[3][2]-P[0][1]*P[3][1]*P[4][2]*P[1][0]+P[0][1]*P[2][1]*P[3][0]*P[4][2]+P[0][1]*P[2][1]*P[1][2]
    *P[4][0]+P[0][1]*P[4][1]*P[2][2]*P[3][0]-P[0][1]*P[4][1]*P[3][0]*P[1][2]-P[0][1]*P[3][1]*P[2][2]*P[4][0]+P[0][1]*P[3][1]*P[2][2]*P[1][0]-P[0][1]*P[2][2]*P[4][1]*P[1][0]-P[0][1]*P[2][1]*P[4][0]*P[3][2]+P[0][1]*
    P[4][1]*P[1][0]*P[3][2]-P[0][1]*P[2][1]*P[3][0]*P[1][2]-P[0][1]*P[4][2]*P[2][1]*P[1][0]+P[0][1]*P[3][1]*P[1][2]*P[4][0]+P[0][1]*P[2][0]*P[3][1]*P[4][2]+P[0][1]*P[2][2]*P[1][1]*P[4][0]+P[0][1]*P[1][1]*P[4][2]*
    P[3][0]-P[0][1]*P[1][1]*P[2][2]*P[3][0]-P[0][1]*P[1][1]*P[3][2]*P[4][0]+P[0][1]*P[2][0]*P[1][1]*P[3][2]+P[0][1]*P[4][1]*P[1][2]*P[2][0]-P[0][1]*P[2][0]*P[3][1]*P[1][2]-P[0][1]*P[2][0]*P[4][1]*P[3][2]-P[0][1]*
    P[4][2]*P[1][1]*P[2][0]-two*P[2][1]*P[4][0]*P[3][1]*P[1][2]+two*P[0][0]*P[2][1]*P[4][2]*P[1][1]+two*P[2][1]*P[4][0]*P[3][2]*P[1][1]+two*P[2][1]*P[3][0]*P[4][1]*P[1][2]+two*P[0][0]
    *P[2][1]*P[4][1]*P[3][2]-two*P[2][1]*P[3][0]*P[4][2]*P[1][1]+two*P[0][2]*P[2][1]*P[4][1]*P[1][0]+two*P[0][2]*P[2][1]*P[3][1]*P[4][0]-two*P[0][2]*P[2][1]*P[4][1]*P[3][0]+two*
    P[0][2]*P[2][1]*P[1][1]*P[3][0]-two*P[0][2]*P[2][1]*P[3][1]*P[1][0]-two*P[0][2]*P[2][1]*P[1][1]*P[4][0]+two*P[2][1]*P[1][0]*P[3][1]*P[4][2]-two*P[0][0]*P[2][1]*P[3][2]*P[1][1]-two
    *P[0][0]*P[2][1]*P[4][1]*P[1][2]+two*P[0][0]*P[2][1]*P[3][1]*P[1][2]-two*P[0][0]*P[2][1]*P[3][1]*P[4][2])*r;
  M[0][2] = -(P[0][2]*P[2][2]*P[4][1]*P[1][0]-P[0][2]*P[4][1]*P[3][0]*P[1][2]+P[0][2]*P[3][1]*P[2][2]*P[4][0]-P[0][2]*P[3][1]*P[2][2]*P[1][0]-P[0][2]*P[2][1]*P[1][0]*P[3][2]-
    P[0][2]*P[3][1]*P[4][2]*P[1][0]-P[0][2]*P[2][1]*P[3][0]*P[4][2]-P[0][2]*P[2][1]*P[1][2]*P[4][0]-P[0][2]*P[2][0]*P[4][1]*P[3][2]+P[0][2]*P[2][0]*P[3][1]*P[4][2]-P[0][2]*P[2][2]*P[1][1]*P[4][0]+P[0][2]*P[1][1]*
    P[4][2]*P[3][0]+P[0][2]*P[1][1]*P[2][2]*P[3][0]-P[0][2]*P[1][1]*P[3][2]*P[4][0]+P[0][2]*P[2][1]*P[4][0]*P[3][2]+P[0][2]*P[4][1]*P[1][0]*P[3][2]+P[0][2]*P[2][1]*P[3][0]*P[1][2]+P[0][2]*P[4][2]*P[2][1]*P[1][0]+
    P[0][2]*P[3][1]*P[1][2]*P[4][0]-P[0][2]*P[4][1]*P[2][2]*P[3][0]-P[0][2]*P[4][2]*P[1][1]*P[2][0]+P[0][2]*P[2][0]*P[1][1]*P[3][2]+P[0][2]*P[4][1]*P[1][2]*P[2][0]-P[0][2]*P[2][0]*P[3][1]*P[1][2]+two*P[0][1]*
    P[3][2]*P[2][2]*P[1][0]+two*P[3][2]*P[4][0]*P[2][2]*P[1][1]-two*P[0][1]*P[3][2]*P[2][2]*P[4][0]+two*P[0][0]*P[3][2]*P[2][2]*P[4][1]-two*P[3][2]*P[1][0]*P[2][2]*P[4][1]-two*P[0][0]
    *P[3][2]*P[2][2]*P[1][1]-two*P[0][1]*P[2][2]*P[4][2]*P[1][0]+two*P[0][1]*P[2][2]*P[4][2]*P[3][0]+two*P[0][1]*P[2][2]*P[1][2]*P[4][0]-two*P[2][2]*P[3][0]*P[4][2]*P[1][1]+two*
    P[2][2]*P[3][0]*P[4][1]*P[1][2]-two*P[2][2]*P[4][0]*P[3][1]*P[1][2]-two*P[0][1]*P[2][2]*P[1][2]*P[3][0]+two*P[0][0]*P[2][2]*P[4][2]*P[1][1]-two*P[0][0]*P[2][2]*P[4][1]*P[1][2]+two
    *P[0][0]*P[2][2]*P[3][1]*P[1][2]-two*P[0][0]*P[2][2]*P[3][1]*P[4][2]+two*P[2][2]*P[1][0]*P[3][1]*P[4][2])*r;
  M[0][3] = -(two*P[0][0]*P[4][2]*P[1][1]-two*P[0][0]*P[4][1]*P[1][2]+two*P[0][1]*P[1][2]*P[4][0]-two*P[0][1]*P[4][2]*P[1][0]-two*P[0][2]*P[1][1]
    *P[4][0]+two*P[0][2]*P[4][1]*P[1][0]+two*P[0][0]*P[3][2]*P[4][1]-two*P[0][1]*P[3][2]*P[4][0]-two*P[0][2]*P[4][1]*P[3][0]+two*P[0][2]*P[3][1]*P[4][0]-two*P[0][0]*P[3][2]*
    P[1][1]+two*P[0][2]*P[1][1]*P[3][0]+two*P[0][1]*P[3][2]*P[1][0]-two*P[0][2]*P[3][1]*P[1][0]-P[2][2]*P[4][1]*P[1][0]+P[2][2]*P[1][1]*P[4][0]-P[4][2]*P[1][1]*P[2][0]+P[4][2]*P[2][1]*P[1][0]
    -two*P[0][1]*P[1][2]*P[3][0]+two*P[0][1]*P[4][2]*P[3][0]-two*P[0][0]*P[3][1]*P[4][2]+two*P[0][0]*P[3][1]*P[1][2]+P[4][1]*P[1][2]*P[2][0]-P[2][1]*P[1][2]*P[4][0]-P[3][1]*P[1][2]*P[4][0]
    +P[4][1]*P[3][0]*P[1][2]+P[2][1]*P[3][0]*P[1][2]-P[2][0]*P[3][1]*P[1][2]+P[4][1]*P[2][2]*P[3][0]-P[1][1]*P[2][2]*P[3][0]-P[3][1]*P[2][2]*P[4][0]+P[3][1]*P[2][2]*P[1][0]+P[1][1]*P[3][2]*P[4][0]-P[2][1]*P[1][0]*
    P[3][2]-P[4][1]*P[1][0]*P[3][2]+P[2][1]*P[4][0]*P[3][2]+P[2][0]*P[1][1]*P[3][2]-P[2][0]*P[4][1]*P[3][2]-P[1][1]*P[4][2]*P[3][0]+P[3][1]*P[4][2]*P[1][0]-P[2][1]*P[3][0]*P[4][2]+P[2][0]*P[3][1]*P[4][2])*r;
  M[1][0] = (-two*P[1][0]*P[3][1]*P[2][2]*P[4][0]-two*P[0][1]*P[1][0]*P[2][2]*P[3][0]-two*P[0][1]*P[1][0]*P[3][2]*P[4][0]-two*P[0][2]*P[1][0]*P[4][1]*
    P[3][0]-two*P[0][2]*P[1][0]*P[2][1]*P[4][0]+two*P[0][2]*P[1][0]*P[2][1]*P[3][0]+two*P[0][2]*P[1][0]*P[3][1]*P[4][0]+two*P[1][0]*P[2][1]*P[4][0]*P[3][2]-two*P[1][0]*P[2][1]*P[3][0]
    *P[4][2]+two*P[1][0]*P[4][1]*P[2][2]*P[3][0]+two*P[0][1]*P[1][0]*P[2][2]*P[4][0]+two*P[0][1]*P[1][0]*P[4][2]*P[3][0]+two*P[0][2]*P[2][0]*P[4][1]*P[1][0]-two*P[0][2]*P[2][0]*
    P[3][1]*P[1][0]+two*P[0][1]*P[2][0]*P[3][2]*P[1][0]-two*P[0][1]*P[2][0]*P[4][2]*P[1][0]-P[0][0]*P[2][1]*P[1][0]*P[3][2]-P[0][0]*P[3][1]*P[4][2]*P[1][0]+P[0][0]*P[2][1]*P[3][0]*P[4][2]+P[0][0]*P[2][1]
    *P[1][2]*P[4][0]-two*P[1][0]*P[2][0]*P[4][1]*P[3][2]+two*P[1][0]*P[2][0]*P[3][1]*P[4][2]+P[0][0]*P[4][1]*P[3][0]*P[1][2]+P[0][0]*P[3][1]*P[2][2]*P[4][0]+P[0][0]*P[3][1]*P[2][2]*P[1][0]+P[0][0]*
    P[2][0]*P[3][1]*P[1][2]+P[0][0]*P[2][0]*P[4][1]*P[3][2]-P[0][0]*P[2][0]*P[3][1]*P[4][2]-P[0][0]*P[2][2]*P[1][1]*P[4][0]-P[0][0]*P[1][1]*P[4][2]*P[3][0]+P[0][0]*P[1][1]*P[2][2]*P[3][0]+P[0][0]*P[1][1]*P[3][2]*
    P[4][0]-P[0][0]*P[2][1]*P[4][0]*P[3][2]+P[0][0]*P[4][1]*P[1][0]*P[3][2]-P[0][0]*P[2][1]*P[3][0]*P[1][2]+P[0][0]*P[4][2]*P[2][1]*P[1][0]-P[0][0]*P[3][1]*P[1][2]*P[4][0]-P[0][0]*P[4][1]*P[2][2]*P[3][0]-P[0][0]*
    P[2][2]*P[4][1]*P[1][0]+P[0][0]*P[4][2]*P[1][1]*P[2][0]-P[0][0]*P[2][0]*P[1][1]*P[3][2]-P[0][0]*P[4][1]*P[1][2]*P[2][0])*r;
  M[1][1] = (P[0][1]*P[2][1]*P[1][0]*P[3][2]+P[0][1]*P[3][1]*P[4][2]*P[1][0]+P[0][1]*P[2][1]*P[3][0]*P[4][2]+P[0][1]*P[2][1]*P[1][2]*P[4][0]-P[0][1]*P[4][1]*P[2][2]*P[3][0]+P[0][1]
    *P[4][1]*P[3][0]*P[1][2]+P[0][1]*P[3][1]*P[2][2]*P[4][0]-P[0][1]*P[3][1]*P[2][2]*P[1][0]+P[0][1]*P[2][2]*P[4][1]*P[1][0]-P[0][1]*P[2][1]*P[4][0]*P[3][2]-P[0][1]*P[4][1]*P[1][0]*P[3][2]-P[0][1]*P[2][1]*P[3][0]*
    P[1][2]-P[0][1]*P[4][2]*P[2][1]*P[1][0]-P[0][1]*P[3][1]*P[1][2]*P[4][0]-P[0][1]*P[2][0]*P[3][1]*P[4][2]+P[0][1]*P[2][2]*P[1][1]*P[4][0]+P[0][1]*P[1][1]*P[4][2]*P[3][0]-P[0][1]*P[1][1]*P[2][2]*P[3][0]-P[0][1]*
    P[1][1]*P[3][2]*P[4][0]+P[0][1]*P[2][0]*P[1][1]*P[3][2]-P[0][1]*P[4][1]*P[1][2]*P[2][0]+P[0][1]*P[2][0]*P[3][1]*P[1][2]+P[0][1]*P[2][0]*P[4][1]*P[3][2]-P[0][1]*P[4][2]*P[1][1]*P[2][0]+two*P[0][0]*P[2][1]*
    P[4][2]*P[1][1]+two*P[2][1]*P[4][0]*P[3][2]*P[1][1]-two*P[2][1]*P[3][0]*P[4][2]*P[1][1]+two*P[0][2]*P[2][1]*P[1][1]*P[3][0]-two*P[0][2]*P[2][1]*P[1][1]*P[4][0]-two*P[0][0]*P[2][1]
    *P[3][2]*P[1][1]+two*P[0][2]*P[1][1]*P[4][1]*P[2][0]-two*P[0][2]*P[1][1]*P[3][1]*P[2][0]-two*P[0][2]*P[1][1]*P[4][1]*P[3][0]+two*P[0][2]*P[1][1]*P[3][1]*P[4][0]-two*P[1][1]*
    P[4][0]*P[2][2]*P[3][1]+two*P[1][1]*P[3][0]*P[2][2]*P[4][1]-two*P[1][1]*P[2][0]*P[4][1]*P[3][2]+two*P[1][1]*P[2][0]*P[3][1]*P[4][2]-two*P[0][0]*P[1][1]*P[3][1]*P[4][2]-two*P[0][0]
    *P[1][1]*P[2][2]*P[4][1]+two*P[0][0]*P[1][1]*P[2][2]*P[3][1]+two*P[0][0]*P[1][1]*P[4][1]*P[3][2])*r;
  M[1][2] = (P[0][2]*P[2][2]*P[4][1]*P[1][0]-P[0][2]*P[4][1]*P[3][0]*P[1][2]+P[0][2]*P[3][1]*P[2][2]*P[4][0]-P[0][2]*P[3][1]*P[2][2]*P[1][0]+P[0][2]*P[2][1]*P[1][0]*P[3][2]+P[0][2]
    *P[3][1]*P[4][2]*P[1][0]+P[0][2]*P[2][1]*P[3][0]*P[4][2]-P[0][2]*P[2][1]*P[1][2]*P[4][0]+P[0][2]*P[2][0]*P[4][1]*P[3][2]-P[0][2]*P[2][0]*P[3][1]*P[4][2]-P[0][2]*P[2][2]*P[1][1]*P[4][0]-P[0][2]*P[1][1]*P[4][2]*
    P[3][0]+P[0][2]*P[1][1]*P[2][2]*P[3][0]+P[0][2]*P[1][1]*P[3][2]*P[4][0]-P[0][2]*P[2][1]*P[4][0]*P[3][2]-P[0][2]*P[4][1]*P[1][0]*P[3][2]+P[0][2]*P[2][1]*P[3][0]*P[1][2]-P[0][2]*P[4][2]*P[2][1]*P[1][0]+P[0][2]*
    P[3][1]*P[1][2]*P[4][0]-P[0][2]*P[4][1]*P[2][2]*P[3][0]+P[0][2]*P[4][2]*P[1][1]*P[2][0]-P[0][2]*P[2][0]*P[1][1]*P[3][2]+P[0][2]*P[4][1]*P[1][2]*P[2][0]-P[0][2]*P[2][0]*P[3][1]*P[1][2]-two*P[0][1]*P[3][2]*
    P[1][2]*P[4][0]+two*P[3][2]*P[4][0]*P[2][1]*P[1][2]+two*P[0][1]*P[3][2]*P[1][2]*P[2][0]-two*P[2][0]*P[3][2]*P[4][1]*P[1][2]-two*P[0][0]*P[3][2]*P[2][1]*P[1][2]+two*P[0][0]*P[3][2]
    *P[4][1]*P[1][2]+two*P[0][1]*P[2][2]*P[1][2]*P[4][0]+two*P[2][2]*P[3][0]*P[4][1]*P[1][2]-two*P[2][2]*P[4][0]*P[3][1]*P[1][2]-two*P[0][1]*P[2][2]*P[1][2]*P[3][0]-two*P[0][0]*
    P[2][2]*P[4][1]*P[1][2]+two*P[0][0]*P[2][2]*P[3][1]*P[1][2]-two*P[0][1]*P[1][2]*P[4][2]*P[2][0]+two*P[0][1]*P[1][2]*P[4][2]*P[3][0]-two*P[0][0]*P[1][2]*P[3][1]*P[4][2]+two*P[0][0]
    *P[1][2]*P[2][1]*P[4][2]+two*P[1][2]*P[2][0]*P[3][1]*P[4][2]-two*P[1][2]*P[3][0]*P[2][1]*P[4][2])*r;
  M[1][3] = (two*P[0][0]*P[2][1]*P[4][2]-two*P[0][0]*P[2][2]*P[4][1]-two*P[0][1]*P[4][2]*P[2][0]+two*P[0][1]*P[2][2]*P[4][0]+two*P[0][0]*P[3][2]*
    P[4][1]+two*P[0][2]*P[4][1]*P[2][0]-two*P[0][1]*P[3][2]*P[4][0]-two*P[0][2]*P[4][1]*P[3][0]-two*P[0][2]*P[2][1]*P[4][0]+two*P[0][2]*P[3][1]*P[4][0]+two*P[0][1]*P[3][2]*
    P[2][0]-two*P[0][0]*P[3][2]*P[2][1]+two*P[0][2]*P[2][1]*P[3][0]-two*P[0][2]*P[2][0]*P[3][1]-two*P[0][1]*P[2][2]*P[3][0]+P[2][2]*P[4][1]*P[1][0]-P[2][2]*P[1][1]*P[4][0]+P[4][2]*P[1][1]*
    P[2][0]-P[4][2]*P[2][1]*P[1][0]+two*P[0][0]*P[3][1]*P[2][2]+two*P[0][1]*P[4][2]*P[3][0]-two*P[0][0]*P[3][1]*P[4][2]-P[4][1]*P[1][2]*P[2][0]+P[2][1]*P[1][2]*P[4][0]-P[3][1]*P[1][2]*P[4][0]+
    P[4][1]*P[3][0]*P[1][2]-P[2][1]*P[3][0]*P[1][2]+P[2][0]*P[3][1]*P[1][2]+P[4][1]*P[2][2]*P[3][0]+P[1][1]*P[2][2]*P[3][0]-P[3][1]*P[2][2]*P[4][0]-P[3][1]*P[2][2]*P[1][0]+P[1][1]*P[3][2]*P[4][0]+P[2][1]*P[1][0]*
    P[3][2]-P[4][1]*P[1][0]*P[3][2]+P[2][1]*P[4][0]*P[3][2]-P[2][0]*P[1][1]*P[3][2]-P[2][0]*P[4][1]*P[3][2]-P[1][1]*P[4][2]*P[3][0]+P[3][1]*P[4][2]*P[1][0]-P[2][1]*P[3][0]*P[4][2]+P[2][0]*P[3][1]*P[4][2])*r;
  M[2][0] = (two*P[0][1]*P[3][0]*P[2][2]*P[4][0]+two*P[3][0]*P[2][1]*P[1][2]*P[4][0]-two*P[0][1]*P[3][0]*P[1][2]*P[4][0]+two*P[0][2]*P[3][0]*P[1][1]*
    P[4][0]-two*P[0][2]*P[3][0]*P[2][1]*P[4][0]-two*P[0][1]*P[1][0]*P[2][2]*P[3][0]-two*P[0][2]*P[1][0]*P[4][1]*P[3][0]+two*P[0][2]*P[1][0]*P[2][1]*P[3][0]-two*P[1][0]*P[2][1]*P[3][0]
    *P[4][2]+two*P[1][0]*P[4][1]*P[2][2]*P[3][0]+two*P[0][1]*P[1][0]*P[4][2]*P[3][0]-two*P[0][2]*P[2][0]*P[1][1]*P[3][0]+two*P[0][2]*P[2][0]*P[4][1]*P[3][0]+two*P[0][1]*P[2][0]*
    P[1][2]*P[3][0]-two*P[0][1]*P[2][0]*P[4][2]*P[3][0]-two*P[2][0]*P[4][1]*P[3][0]*P[1][2]+two*P[2][0]*P[1][1]*P[4][2]*P[3][0]+P[0][0]*P[2][1]*P[1][0]*P[3][2]+P[0][0]*P[3][1]*P[4][2]*P[1][0]+
    three*P[0][0]*P[2][1]*P[3][0]*P[4][2]+P[0][0]*P[2][1]*P[1][2]*P[4][0]+three*P[0][0]*P[4][1]*P[3][0]*P[1][2]+P[0][0]*P[3][1]*P[2][2]*P[4][0]-P[0][0]*P[3][1]*P[2][2]*P[1][0]+P[0][0]*P[2][0]*P[3][1]*P[1][2]
    +P[0][0]*P[2][0]*P[4][1]*P[3][2]-P[0][0]*P[2][0]*P[3][1]*P[4][2]-P[0][0]*P[2][2]*P[1][1]*P[4][0]-three*P[0][0]*P[1][1]*P[4][2]*P[3][0]+three*P[0][0]*P[1][1]*P[2][2]*P[3][0]+P[0][0]*P[1][1]*P[3][2]*
    P[4][0]-P[0][0]*P[2][1]*P[4][0]*P[3][2]-P[0][0]*P[4][1]*P[1][0]*P[3][2]-three*P[0][0]*P[2][1]*P[3][0]*P[1][2]-P[0][0]*P[4][2]*P[2][1]*P[1][0]-P[0][0]*P[3][1]*P[1][2]*P[4][0]-three*P[0][0]*P[4][1]*P[2][2]
    *P[3][0]+P[0][0]*P[2][2]*P[4][1]*P[1][0]+P[0][0]*P[4][2]*P[1][1]*P[2][0]-P[0][0]*P[2][0]*P[1][1]*P[3][2]-P[0][0]*P[4][1]*P[1][2]*P[2][0]-two*P[3][0]*P[2][2]*P[1][1]*P[4][0])*r;
  M[2][1] = (P[0][1]*P[2][1]*P[1][0]*P[3][2]+three*P[0][1]*P[3][1]*P[4][2]*P[1][0]+P[0][1]*P[2][1]*P[3][0]*P[4][2]+P[0][1]*P[2][1]*P[1][2]*P[4][0]-P[0][1]*P[4][1]*P[2][2]*
    P[3][0]+P[0][1]*P[4][1]*P[3][0]*P[1][2]+three*P[0][1]*P[3][1]*P[2][2]*P[4][0]-three*P[0][1]*P[3][1]*P[2][2]*P[1][0]+P[0][1]*P[2][2]*P[4][1]*P[1][0]-P[0][1]*P[2][1]*P[4][0]*P[3][2]-P[0][1]*P[4][1]*P[1][0]
    *P[3][2]-P[0][1]*P[2][1]*P[3][0]*P[1][2]-P[0][1]*P[4][2]*P[2][1]*P[1][0]-three*P[0][1]*P[3][1]*P[1][2]*P[4][0]-three*P[0][1]*P[2][0]*P[3][1]*P[4][2]-P[0][1]*P[2][2]*P[1][1]*P[4][0]-P[0][1]*P[1][1]*
    P[4][2]*P[3][0]+P[0][1]*P[1][1]*P[2][2]*P[3][0]+P[0][1]*P[1][1]*P[3][2]*P[4][0]-P[0][1]*P[2][0]*P[1][1]*P[3][2]-P[0][1]*P[4][1]*P[1][2]*P[2][0]+three*P[0][1]*P[2][0]*P[3][1]*P[1][2]+P[0][1]*P[2][0]*P[4][1]*
    P[3][2]+P[0][1]*P[4][2]*P[1][1]*P[2][0]+two*P[2][1]*P[4][0]*P[3][1]*P[1][2]-two*P[0][2]*P[2][1]*P[3][1]*P[4][0]+two*P[0][2]*P[2][1]*P[3][1]*P[1][0]-two*P[2][1]*P[1][0]*P[3][1]*P[4][2]
    -two*P[0][0]*P[2][1]*P[3][1]*P[1][2]+two*P[0][0]*P[2][1]*P[3][1]*P[4][2]-two*P[0][2]*P[1][1]*P[3][1]*P[2][0]+two*P[0][2]*P[1][1]*P[3][1]*P[4][0]-two*P[1][1]*P[4][0]*P[2][2]*
    P[3][1]+two*P[1][1]*P[2][0]*P[3][1]*P[4][2]-two*P[0][0]*P[1][1]*P[3][1]*P[4][2]+two*P[0][0]*P[1][1]*P[2][2]*P[3][1]+two*P[0][2]*P[3][1]*P[4][1]*P[2][0]-two*P[0][2]*P[3][1]*P[4][1]
    *P[1][0]-two*P[3][1]*P[2][0]*P[4][1]*P[1][2]+two*P[3][1]*P[1][0]*P[2][2]*P[4][1]-two*P[0][0]*P[3][1]*P[2][2]*P[4][1]+two*P[0][0]*P[3][1]*P[4][1]*P[1][2])*r;
  M[2][2] = -(-P[0][2]*P[2][2]*P[4][1]*P[1][0]-P[0][2]*P[4][1]*P[3][0]*P[1][2]-P[0][2]*P[3][1]*P[2][2]*P[4][0]+P[0][2]*P[3][1]*P[2][2]*P[1][0]-three*P[0][2]*P[2][1]*P[1][0]
    *P[3][2]-P[0][2]*P[3][1]*P[4][2]*P[1][0]-P[0][2]*P[2][1]*P[3][0]*P[4][2]-P[0][2]*P[2][1]*P[1][2]*P[4][0]-three*P[0][2]*P[2][0]*P[4][1]*P[3][2]+P[0][2]*P[2][0]*P[3][1]*P[4][2]+P[0][2]*P[2][2]*P[1][1]*P[4][0]
    +P[0][2]*P[1][1]*P[4][2]*P[3][0]-P[0][2]*P[1][1]*P[2][2]*P[3][0]-three*P[0][2]*P[1][1]*P[3][2]*P[4][0]+three*P[0][2]*P[2][1]*P[4][0]*P[3][2]+three*P[0][2]*P[4][1]*P[1][0]*P[3][2]+P[0][2]*P[2][1]*
    P[3][0]*P[1][2]+P[0][2]*P[4][2]*P[2][1]*P[1][0]+P[0][2]*P[3][1]*P[1][2]*P[4][0]+P[0][2]*P[4][1]*P[2][2]*P[3][0]-P[0][2]*P[4][2]*P[1][1]*P[2][0]+three*P[0][2]*P[2][0]*P[1][1]*P[3][2]+P[0][2]*P[4][1]*P[1][2]*
    P[2][0]-P[0][2]*P[2][0]*P[3][1]*P[1][2]+two*P[0][1]*P[3][2]*P[2][2]*P[1][0]+two*P[0][1]*P[3][2]*P[1][2]*P[4][0]-two*P[0][1]*P[3][2]*P[4][2]*P[1][0]-two*P[3][2]*P[4][0]*P[2][1]*P[1][2]+
    two*P[3][2]*P[4][0]*P[2][2]*P[1][1]-two*P[0][1]*P[3][2]*P[1][2]*P[2][0]+two*P[0][1]*P[3][2]*P[4][2]*P[2][0]-two*P[0][1]*P[3][2]*P[2][2]*P[4][0]+two*P[2][0]*P[3][2]*P[4][1]*P[1][2]
    +two*P[3][2]*P[1][0]*P[2][1]*P[4][2]-two*P[2][0]*P[3][2]*P[4][2]*P[1][1]+two*P[0][0]*P[3][2]*P[2][1]*P[1][2]-two*P[0][0]*P[3][2]*P[4][1]*P[1][2]+two*P[0][0]*P[3][2]*P[2][2]*
    P[4][1]-two*P[0][0]*P[3][2]*P[2][1]*P[4][2]-two*P[3][2]*P[1][0]*P[2][2]*P[4][1]-two*P[0][0]*P[3][2]*P[2][2]*P[1][1]+two*P[0][0]*P[3][2]*P[4][2]*P[1][1])*r;
  M[2][3] = -(-two*P[0][1]*P[2][0]*P[1][2]+two*P[0][2]*P[2][0]*P[1][1]-two*P[0][2]*P[1][0]*P[2][1]+two*P[0][1]*P[1][0]*P[2][2]+two*P[0][0]*
    P[2][1]*P[1][2]-two*P[0][0]*P[1][1]*P[2][2]+two*P[0][0]*P[4][2]*P[1][1]-two*P[0][0]*P[2][1]*P[4][2]-two*P[0][0]*P[4][1]*P[1][2]+two*P[0][0]*P[2][2]*P[4][1]+two*P[0][1]*
    P[1][2]*P[4][0]+two*P[0][1]*P[4][2]*P[2][0]-two*P[0][1]*P[2][2]*P[4][0]-two*P[0][1]*P[4][2]*P[1][0]-two*P[0][2]*P[1][1]*P[4][0]+two*P[0][2]*P[4][1]*P[1][0]-two*P[0][2]*
    P[4][1]*P[2][0]+two*P[0][2]*P[2][1]*P[4][0]-three*P[2][2]*P[4][1]*P[1][0]+three*P[2][2]*P[1][1]*P[4][0]-three*P[4][2]*P[1][1]*P[2][0]+three*P[4][2]*P[2][1]*P[1][0]+three*P[4][1]*
    P[1][2]*P[2][0]-three*P[2][1]*P[1][2]*P[4][0]+P[3][1]*P[1][2]*P[4][0]-P[4][1]*P[3][0]*P[1][2]+P[2][1]*P[3][0]*P[1][2]-P[2][0]*P[3][1]*P[1][2]+P[4][1]*P[2][2]*P[3][0]-P[1][1]*P[2][2]*P[3][0]-P[3][1]*P[2][2]*
    P[4][0]+P[3][1]*P[2][2]*P[1][0]-P[1][1]*P[3][2]*P[4][0]-P[2][1]*P[1][0]*P[3][2]+P[4][1]*P[1][0]*P[3][2]+P[2][1]*P[4][0]*P[3][2]+P[2][0]*P[1][1]*P[3][2]-P[2][0]*P[4][1]*P[3][2]+P[1][1]*P[4][2]*P[3][0]-P[3][1]*
    P[4][2]*P[1][0]-P[2][1]*P[3][0]*P[4][2]+P[2][0]*P[3][1]*P[4][2])*r;
  M[3][0] = P[0][0];
  M[3][1] = P[0][1];
  M[3][2] = P[0][2];
  M[3][3] = 1.0;
#else
  /* Agressive common sub-expressions folding via temporary variables. */
  FPvalue t1, t2, t3, t4, t5, t6, t7, t9, t10, t11, t12, t13, t15, t16, t17, t18;
  FPvalue t19, t21, t22, t25, t28, t30, t31, t34, t36, t37, t38, t39, t40, t41;
  FPvalue t42, t43, t44, t45, t46, t47, t48, t50, t51, t53, t54, t55, t56, t58;
  FPvalue t59, t60, t61, t65, t66, t67, t68, t69, t70, t71, t72, t74, t75, t76;
  FPvalue t77, t78, t79, t80, t81, t82, t83, t84, t85, t86, t87, t88, t89, t90;
  FPvalue t93, t94, t95, t96, t97, t98, t99, t100, t101, t102, t103, t104, t105;
  FPvalue t106, t107, t108, t109, t110, t111, t112, t113, t114, t115, t116, t117;
  FPvalue t118, t119, t120, t121, t122, t123, t124, t125, t126, t127, t129, t130;
  FPvalue t131, t134, t135, t136, t137, t138, t139, t140, t141, t142, t143, t144;
  FPvalue t146, t147, t148, t150, t152, t153, t154, t155, t156, t157, t158, t159;
  FPvalue t160, t162, t163, t164, t165, t166, t167, t168, t169, t171, t173, t175;
  FPvalue t179, t181, t182, t185, t186, t187, t189, t190, t195, t196, t197, t198;
  FPvalue t199, t200, t201, t202, t203, t204, t205, t206, t207, t208, t210, t211;
  FPvalue t212, t213, t214, t219, t221, t223, t225, t226, t228, t229, t232, t233;
  FPvalue t235, t236, t237, t238, t239, t240, t241, t242, t243, t245, t249, t255;
  FPvalue t257, t259, t261, t265, t266, t272, t273, t276, t278, t280, t282, t284;
  FPvalue t285, t287, t289, t292, t294, t306, t310, t311, t319, t321, t325, t326;
  FPvalue t327, t329, t330, t331, t336, t337, t342, t344, t345, t348, t349, t353;
  FPvalue t357, t364, t367, t370, t371, t372, t373, t374, t375, t376, t377, t386;
  FPvalue t389, t393, t404, t411, t413, t415, t417, t419, t420, t432, t436, t441;
  FPvalue t457, t474, t482, t498, t504, t513, t530;
  const FPvalue two = (FPvalue) 2.0;
  const FPvalue three = (FPvalue) 3.0;
  
  t1 = P[0][2]*P[2][0];
  t2 = P[1][1]*P[3][0];
  t3 = t1*t2;
  t4 = P[0][1]*P[2][0];
  t5 = P[3][2]*P[1][0];


  t6 = t4*t5;
  t7 = P[1][2]*P[4][0];
  t9 = P[1][2]*P[3][0];
  t10 = t4*t9;
  t11 = P[4][1]*P[3][0];
  t12 = t1*t11;
  t13 = P[3][1]*P[4][0];
  t15 = P[3][1]*P[1][0];
  t16 = t1*t15;
  t17 = P[4][1]*P[1][0];
  t18 = t1*t17;
  t19 = P[1][1]*P[4][0];
  t21 = P[4][1]*P[2][0];
  t22 = t21*t9;
  t25 = P[3][2]*P[4][0];
  t28 = P[4][2]*P[3][0];
  t30 = two*t4*t28;
  t31 = P[3][1]*P[2][0];
  t34 = P[4][2]*P[1][0];
  t36 = two*t4*t34;
  t37 = P[0][0]*P[3][1];
  t38 = P[2][2]*P[1][0];
  t39 = t37*t38;
  t40 = P[0][0]*P[2][0];
  t41 = P[3][1]*P[1][2];
  t42 = t40*t41;
  t43 = P[0][0]*P[4][1];
  t44 = t43*t9;
  t45 = P[2][2]*P[4][0];
  t46 = t37*t45;
  t47 = P[1][0]*P[2][0];
  t48 = P[4][1]*P[3][2];
  t50 = two*t47*t48;
  t51 = P[3][1]*P[4][2];
  t53 = two*t47*t51;
  t54 = P[0][0]*P[2][1];
  t55 = t54*t28;
  t56 = two*t4*t25-t30+two*t31*t7+t36-t39-t42+t44+t46+t50-t53+t55;
  t58 = t54*t7;
  t59 = t37*t34;
  t60 = t54*t5;
  t61 = P[1][1]*P[2][0];
  t65 = two*t61*t28;
  t66 = t54*t25;
  t67 = t43*t5;
  t68 = P[0][0]*P[1][1];
  t69 = P[2][2]*P[3][0];
  t70 = t68*t69;
  t71 = t68*t25;
  t72 = t68*t28;
  t74 = t40*t51;
  t75 = P[0][0]*P[2][2];
  t76 = t75*t19;
  t77 = t40*t48;
  t78 = P[3][2]*P[1][1];
  t79 = t40*t78;
  t80 = P[1][2]*P[2][0];
  t81 = t43*t80;
  t82 = t75*t17;
  t83 = P[0][0]*P[4][2];
  t84 = t83*t61;
  t85 = t37*t7;
  t86 = t43*t69;
  t87 = t54*t9;
  t88 = P[2][1]*P[1][0];
  t89 = t83*t88;
  t90 = t74-t76-t77+t79+t81+t82-t84-t85-t86-t87-t89;
  t93 = P[4][2]*P[1][1];
  t94 = t93*P[2][0];
  t95 = t61*P[3][2];
  t96 = P[4][1]*P[1][2];
  t97 = t96*P[2][0];
  t98 = t31*P[1][2];
  t99 = t21*P[3][2];
  t100 = t31*P[4][2];
  t101 = P[2][2]*P[1][1];
  t102 = t101*P[4][0];
  t103 = t93*P[3][0];
  t104 = t101*P[3][0];
  t105 = t78*P[4][0];
  t106 = P[2][1]*P[4][0];
  t107 = t106*P[3][2];
  t108 = t17*P[3][2];
  t109 = -t94+t95+t97-t98-t99+t100+t102+t103-t104-t105+t107+t108;
  t110 = P[2][1]*P[3][0];
  t111 = t110*P[1][2];
  t112 = P[2][1]*P[4][2];
  t113 = t112*P[1][0];
  t114 = t41*P[4][0];
  t115 = P[2][2]*P[4][1];
  t116 = t115*P[3][0];
  t117 = t115*P[1][0];
  t118 = t11*P[1][2];
  t119 = P[2][2]*P[3][1];
  t120 = t119*P[4][0];
  t121 = t119*P[1][0];
  t122 = t88*P[3][2];
  t123 = t51*P[1][0];
  t124 = t110*P[4][2];
  t125 = P[2][1]*P[1][2];
  t126 = t125*P[4][0];
  t127 = t111+t113+t114+t116-t117-t118-t120+t121-t122-t123-t124-t126;
  t129 = 1/(t109+t127);
  M[0][0] = (-two*t3-two*t6-two*t4*t7+two*t10+two*t12-two*t1*t13+two*t16
    -two*t18+two*t1*t19-two*t22+t56+t58+t59+t60-two*t61*t25+t65-t66-t67+t70+t71-t72
    +t90)*t129;
  t130 = P[0][1]*P[3][1];
  t131 = t130*t34;
  t134 = P[0][1]*P[2][1];
  t135 = t134*t5;
  t136 = t134*t28;
  t137 = P[0][1]*P[4][2];
  t138 = t137*t88;
  t139 = t130*t7;
  t140 = P[0][1]*P[4][1];
  t141 = t140*t5;
  t142 = t134*t9;
  t143 = t130*t45;
  t144 = t130*t38;
  t146 = t140*t9;
  t147 = t134*t7;
  t148 = t140*t69;
  t150 = two*t106*t41;
  t152 = two*t54*t93;
  t153 = t4*t48;
  t154 = t137*t61;
  t155 = t4*t41;
  t156 = P[0][1]*P[1][1];
  t157 = t156*t25;
  t158 = t4*t78;
  t159 = t140*t80;
  t160 = -t146+t147+t148-t150+t152-t153-t154-t155-t157+t158+t159;
  t162 = t156*t28;
  t163 = t156*t69;
  t164 = t4*t51;
  t165 = P[0][1]*P[2][2];
  t166 = t165*t19;
  t167 = t165*t17;
  t168 = t134*t25;
  t169 = P[0][2]*P[2][1];
  t171 = two*t169*t15;
  t173 = two*t169*t19;
  t175 = two*t88*t51;
  t179 = t169*t13;
  t181 = t169*t2;
  t182 = t106*t78;
  t185 = t110*t93;
  t186 = t54*t51;
  t187 = t54*t78;
  t189 = t54*t41;
  t190 = t179-t169*t11+t181+t182+t110*t96+t54*t48-t185-t186-t187-t54*t96+t189;
  M[0][1] = -(-t131-two*t88*t48+t135+t136-t138+t139+t141-t142-t143+t144+
    t160+t162-t163+t164+t166-t167-t168-t171-t173+t175+two*t169*t17+two*t190)*t129;
  t195 = P[0][2]*P[3][1];
  t196 = t195*t38;
  t197 = t195*t34;
  t198 = t169*t5;
  t199 = P[0][2]*P[4][1];
  t200 = t199*t9;
  t201 = t195*t45;
  t202 = P[0][2]*P[2][2];
  t203 = t202*t17;
  t204 = t169*t25;
  t205 = P[0][2]*P[1][1];
  t206 = t205*t25;
  t207 = t205*t28;
  t208 = t205*t69;
  t210 = t202*t19;
  t211 = t1*t48;
  t212 = t1*t51;
  t213 = t169*t28;
  t214 = t169*t7;
  t219 = P[0][0]*P[3][2];
  t221 = two*t219*t101;
  t223 = two*t5*t115;
  t225 = two*t25*t101;
  t226 = P[0][1]*P[3][2];
  t228 = two*t226*t45;
  t229 = -t210-t211+t212-t213-t214-two*t165*t34+two*t165*t28-t221-t223+t225-
    t228;
  t232 = two*t219*t115;
  t233 = t1*t41;
  t235 = two*t226*t38;
  t236 = t199*t5;
  t237 = t169*t9;
  t238 = t1*t78;
  t239 = t199*t80;
  t240 = t199*t69;
  t241 = P[0][2]*P[4][2];
  t242 = t241*t61;
  t243 = t241*t88;
  t245 = t195*t7;
  t249 = two*t75*t41;
  t255 = two*t45*t41;
  t257 = two*t165*t9;
  t259 = two*t75*t96;
  t261 = two*t69*t96;
  t265 = two*t165*t7;
  t266 = t245-two*t75*t51+t249+two*t38*t51+two*t75*t93-t255-t257-t259+t261
    -two*t69*t93+t265;
  M[0][2] = -(-t196-t197-t198-t200+t201+t203+t204-t206+t207+t208+t229+t232-
    t233+t235+t236+t237+t238+t239-t240-t242+t243+t266)*t129;
  t272 = two*t199*P[3][0];
  t273 = -t122+t121-t120+t118-t117+t116-t114+t113+t111-t108-t272;
  t276 = two*t226*P[4][0];
  t278 = two*t219*P[4][1];
  t280 = two*t199*P[1][0];
  t282 = two*t205*P[4][0];
  t284 = two*t137*P[1][0];
  t285 = P[0][1]*P[1][2];
  t287 = two*t285*P[4][0];
  t289 = two*t43*P[1][2];
  t292 = two*t83*P[1][1];
  t294 = two*t137*P[3][0];
  t306 = two*t195*P[4][0];
  t310 = two*t37*P[4][2];
  t311 = t292+t294-two*t285*P[3][0]-two*t195*P[1][0]+two*t226*P[1][0]+two*t205*P[3][0]-t94-two*
    t219*P[1][1]+t306+two*t37*P[1][2]-t310;
  M[0][3] = -(t107+t105-t104-t103+t102+t100-t99-t126-t124+t123+t273-t276+
    t278+t280-t282+t95-t284-t98+t97+t287-t289+t311)*t129;
  t319 = t55+t58-t59-t60-t66+t67+t70+t71-t72-t74-t76;
  t321 = P[0][1]*P[1][0];
  t325 = P[0][2]*P[1][0];
  t326 = t325*t11;
  t327 = t321*t69;
  t329 = t321*t28;
  t330 = t17*t69;
  t331 = t88*t28;
  t336 = t325*t110;
  t337 = -t326-t327-t15*t45+t329+t330-t331+t325*t13+t321*t45+t88*t25-t325*
    t106+t336;
  M[1][0] = (two*t6-two*t16+two*t18-t36+t39+t42+t44+t46-t50+t53+t319+t77-
    t79-t81-t82+t84-t85-t86-t87+t89-two*t321*t25+two*t337)*t129;
  t342 = t205*t31;
  t344 = t205*t13;
  t345 = t19*t119;
  t348 = t61*t51;
  t349 = t68*t51;
  t353 = t68*t119;
  t357 = two*t353+two*t68*t48+t131+t135+t136-t138-t139-t141-t142+t143-t144;
  t364 = t162-t163-t164+t166+t167-t168-t173+two*t181+two*t182-two*t185-two*
    t187;
  M[1][1] = (two*t205*t21-two*t342-two*t205*t11+two*t344-two*t345+two*t2*
    t115-two*t61*t48+two*t348-two*t349-two*t68*t115+t357+t146+t147-t148+t152+t153-
    t154+t155-t157+t158-t159+t364)*t129;
  t367 = P[4][2]*P[2][0];
  t370 = t226*t7;
  t371 = t25*t125;
  t372 = t226*t80;
  t373 = P[2][0]*P[3][2];
  t374 = t373*t96;
  t375 = t219*t125;
  t376 = t219*t96;
  t377 = P[0][0]*P[1][2];
  t386 = two*t80*t51-two*t9*t112-t196+t197+t198-t200+t201+t203-t204+t206-t207
    ;
  t389 = t239-t240+t242-t243+t245+t249-t255-t257-t259+t261+t265;
  M[1][2] = (-two*t285*t367+two*t285*t28-two*t370+two*t371+two*t372-two*
    t374-two*t375+two*t376-two*t377*t51+two*t377*t112+t386+t208-t210+t211-t212+t213
    -t214-t233-t236+t237-t238+t389)*t129;
  t393 = t122-t121-t120+t118+t117+t116-t114-t113-t111-t108-t272;
  t404 = two*t169*P[4][0];
  t411 = two*t165*P[4][0];
  t413 = two*t54*P[4][2];
  t415 = two*t75*P[4][1];
  t417 = two*t137*P[2][0];
  t419 = two*t199*P[2][0];
  t420 = t411+t98-t97+t413-t415-t417+t419+t294+t94+t306-t310;
  M[1][3] = (t107+t105+t104-t103-t102+t100-t99+t126-t124+t123+t393-t276+
    t278-t95+two*t169*P[3][0]-two*t1*P[3][1]-two*t165*P[3][0]+two*t37*P[2][2]-t404+two*t226*P[2][0]-two*t219
    *P[2][1]+t420)*t129;
  t432 = t58+t59+t60+t65-t66-t67+three*t70+t71-three*t72-t74-t76;
  t436 = P[0][2]*P[3][0];
  t441 = P[0][1]*P[3][0];
  t457 = -t85-three*t86-three*t87-t89-two*t69*t19-two*t326-two*t327+two*t329+two*
    t330-two*t331+two*t336;
  M[2][0] = (-two*t3+two*t10+two*t12-two*t22-t30-t39+t42+three*t44+t46+three*
    t55+t432+t77-t79-t81+two*t110*t7+two*t436*t19-two*t436*t106-two*t441*t7+two*
    t441*t45+t82+t84+t457)*t129;
  t474 = -two*t31*t96+two*t15*t115+three*t131+t135+t136-t138-three*t139-t141-t142
    +three*t143-three*t144;
  t482 = -t162+t163-three*t164-t166+t167-t168+t171-t175-two*t179+two*t186-two*
    t189;
  M[2][1] = (-two*t342+two*t344-two*t345+two*t348-two*t349+two*t353-two*t37
    *t115+two*t37*t96+two*t195*t21-two*t195*t17+t474+t146+t147-t148+t150+t153+t154+
    three*t155+t157-t158-t159+t482)*t129;
  t498 = two*t219*t93-two*t226*t34+t196-t197-three*t198-t200-t201-t203+three*t204
    -three*t206+t207;
  t504 = t232-t233+t235+three*t236+t237+three*t238+t239+t240-t242+t243+t245;
  M[2][2] = -(two*t370-two*t371-two*t372+two*t374+two*t375-two*t376+two*
    t226*t367+two*t5*t112-two*t373*t93-two*t219*t112+t498-t208+t210-three*t211+t212-
    t213-t214-t221-t223+t225-t228+t504)*t129;
  t513 = -t122+t121-t120-t118-three*t117+t116+t114+three*t113+t111+t108+t280;
  t530 = -t419+t287-t289+t292-three*t94+two*t54*P[1][2]-two*t68*P[2][2]-two*t4*P[1][2]+two*t1*
    P[1][1]-two*t325*P[2][1]+two*t321*P[2][2];
  M[2][3] = -(t107-t105-t104+t103+three*t102+t100-t99-three*t126-t124-t123+t513
    -t282+t95-t284+t404-t411-t98+three*t97-t413+t415+t417+t530)*t129;
  M[3][0] = P[0][0];
  M[3][1] = P[0][1];
  M[3][2] = P[0][2];
  M[3][3] = 1.0;
#endif
}

/** A few simple double-precision analytical 3D geometry routines. ***/

double
vlengthd(const double *v)
{
    return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}

void
vscaled(double *v, double div)
{
    v[0] *= div;
    v[1] *= div;
    v[2] *= div;
}

void
vnormald(double *v)
{
    vscaled(v,1.0/vlengthd(v));
}

/* GRAPHICS GEMS I, "Useful 3D Geometry", p297 */
void
nearestPointOnPlaneToPoint(double J[4], double P[3], double Q[3])
{
  double JdotP = P[0]*J[0]+P[1]*J[1]+P[2]*J[2];
  double JdotJ = J[0]*J[0]+J[1]*J[1]+J[2]*J[2];
  double k = (J[3] + JdotP)/JdotJ;

  Q[0] = P[0] - k * J[0];
  Q[1] = P[1] - k * J[1];
  Q[2] = P[2] - k * J[2];
}

/* GRAPHICS GEMS I, "Useful 3D Geometry", p299 */
void
pointAtLineAndPlaneIntersection(double J[4], 
                                 double U[3], double V[3],
                                 double P[3])
{
  double t = - (J[3] + (U[0]*J[0]+U[1]*J[1]+U[2]*J[2]))/
               (V[0]*J[0]+V[1]*J[1]+V[2]*J[2]);

  P[0] = U[0] + V[0] * t;
  P[1] = U[1] + V[1] * t;
  P[2] = U[2] + V[2] * t;
}

/* Find the plane equation given 3 points. */
void
findPlane(double plane[4],
  double v0[3], double v1[3], double v2[3])
{
  double vec0[3], vec1[3];

  /* Need 2 vectors to find cross product. */
  vec0[0] = v1[0] - v0[0];
  vec0[1] = v1[1] - v0[1];
  vec0[2] = v1[2] - v0[2];

  vec1[0] = v2[0] - v0[0];
  vec1[1] = v2[1] - v0[1];
  vec1[2] = v2[2] - v0[2];

  /* find cross product to get A, B, and C of plane equation */
  plane[0] = vec0[1] * vec1[2] - vec0[2] * vec1[1];
  plane[1] = -(vec0[0] * vec1[2] - vec0[2] * vec1[0]);

  plane[2] = vec0[0] * vec1[1] - vec0[1] * vec1[0];

  plane[3] = -(plane[0] * v0[0] + plane[1] * v0[1] + plane[2] * v0[2]);
}

/* Compute a 4x4 matrix suitable for glMultMatrixd for projecting
   into a plane. */
void
computeShadowVolume(
  double L[3],  /* IN: light location */
  double O[3],  /* IN: cut-out plane */
  double M[3],  /* IN: cut-out X axis direction */
  double N[3],  /* IN: cut-out Y axis direction */
  double G[4],  /* IN: ground plane */
  double P[5][3])  /* OUT: five mutually non-coplanar points
                      defining the shadow volume projection */
{
  double LO[3];  /* Normalized direction vector from origin to light. */
  double EP[3], LEP[3];

  P[0][0] = O[0];
  P[0][1] = O[1];
  P[0][2] = O[2];

  P[1][0] = O[0] + N[0];
  P[1][1] = O[1] + N[1];
  P[1][2] = O[2] + N[2];

  P[2][0] = O[0] + M[0];
  P[2][1] = O[1] + M[1];
  P[2][2] = O[2] + M[2];

  LO[0] = L[0] - O[0];
  LO[1] = L[1] - O[1];
  LO[2] = L[2] - O[2];
  vnormald(LO);

  pointAtLineAndPlaneIntersection(Pbg, O, LO, P[3]);

  EP[0] = O[0] - N[0] - M[0];
  EP[1] = O[1] - N[1] - M[1];
  EP[2] = O[2] - N[2] - M[2];

  LEP[0] = L[0] - EP[0];
  LEP[1] = L[1] - EP[1];
  LEP[2] = L[2] - EP[2];
  vnormald(LEP);

  pointAtLineAndPlaneIntersection(Pbg, EP, LEP, P[4]);
}

#ifndef GLU_VERSION_1_2
/* nvshadow is written to tesselate the NVIDIA logo from its boundary
   points, but this is done with the GLU 1.2 tessellator.  Most official
   OpenGL implementations provide the GLU 1.2 API, but Mesa 3.1 (and
   earlier) do not.  For this reason, if GLU_VERSION_1_2 is not defined,
   use pre-tesselated version of the NVIDIA logo. */
#define bTriFan     glBegin(GL_TRIANGLE_FAN);
#define bTriStrip   glBegin(GL_TRIANGLE_STRIP);
#define bTris       glBegin(GL_TRIANGLES);
#define bQuadStrip  glBegin(GL_QUAD_STRIP);
#define e           glEnd();
#define v(x,y)      glVertex2f(x,y);
static void preTesselatedCutOut(void) {
bTriFan
v(-1.06485,2.81135) v(-1.06478,2.05894) v(-2.05154,1.85121)
v(-2.83056,1.43573) v(-3.5057,0.9609) v(-5,5) v(5,5) v(3.97311,2.81135) e
bTriStrip
v(-2.39205,0.536978) v(-2.36314,1.18348) v(-2.04704,0.806418)
v(-1.79187,1.49508) v(-1.63302,1.01671) v(-1.06478,1.65831)
v(-1.035,1.16786) e
bTriFan
v(-2.53006,0.260968) v(-2.41505,-0.047902) v(-2.10347,-1.22034)
v(-2.57088,-0.760353) v(-2.83056,-0.315201) v(-2.36314,1.18348)
v(-2.39205,0.536978) e
bTriFan
v(-2.83056,0.84219) v(-2.36314,1.18348) v(-2.83056,-0.315201)
v(-3.09023,0.055759) v(-3.1941,0.471234) e
bTriStrip
v(-2.04704,-0.606495) v(-1.79403,-0.875934) v(-1.58413,-1.53195)
v(-1.47202,-1.09937) v(-1.02549,-1.73179) v(-1.02549,-1.24824) e
bTriFan
v(-2.10347,-1.22034) v(-2.41505,-0.047902) v(-2.25405,-0.337056)
v(-2.04704,-0.606495) v(-1.58413,-1.53195) e
bTriStrip
v(3.97311,2.81135) v(5,5) v(3.97311,-2.81135) v(5,-5) v(-1.02549,-2.81135)
v(-5,-5) v(-2.93443,-1.20551) v(-3.29797,-0.760353) e
bTriFan
v(-1.02549,-2.81135) v(-2.93443,-1.20551) v(-2.41508,-1.63582)
v(-1.8438,-1.96226) v(-1.02549,-2.15912) e
bTriFan
v(-5,-5) v(-5,5) v(-3.97311,0.545426) v(-3.86924,0.204143)
v(-3.60957,-0.315201) v(-3.29797,-0.760353) e
bTriStrip
v(1.29669,-1.07956) v(1.69028,-0.820916) v(1.69028,-1.65308)
v(2.00515,-0.51729) v(2.20194,-1.45066) v(2.4381,-0.179928)
v(2.67425,-1.21451) v(3.34336,-0.719707) v(3.10721,-0.955861) e
bTriStrip
v(-0.474465,-2.15912) v(-0.198952,-1.7093) v(0.155277,-2.10289)
v(0.273356,-1.5856) v(0.706303,-1.97919) v(0.86374,-1.3382)
v(1.17861,-1.833) v(1.29669,-1.07956) v(1.69028,-1.65308) e
bTriFan
v(-1.02549,-2.15912) v(-1.02549,-1.73179) v(-0.198952,-1.7093)
v(-0.474465,-2.15912) e
bTriStrip
v(2.12323,0.494796) v(1.92643,0.269888) v(1.84771,0.83216)
v(1.69028,-0.033738) v(1.53284,1.12454) v(1.3242,0.426719)
v(1.17861,1.39443) v(1.01259,0.753164) v(0.706303,1.7093)
v(0.64905,1.09445) v(0.233576,1.36154) e
bTriStrip
v(-0.441573,-0.87906) v(-0.078033,-0.745515) v(0.115919,-1.1133)
v(0.389378,-0.5081) v(0.588227,-0.899634) v(0.804852,-0.166817)
v(0.942455,-0.652235) v(1.06452,0.189305) v(1.29669,-0.348609)
v(1.3242,0.426719) v(1.69028,-0.033738) e
bTriFan
v(-0.474465,-1.25949) v(-1.02549,-1.24824) v(-1.01353,-0.88907)
v(-0.441573,-0.87906) v(0.115919,-1.1133) e
bTriStrip
v(0.706303,1.7093) v(0.233576,1.36154) v(0.076562,1.96794)
v(-0.233835,1.59895) v(-0.395747,2.0804) v(-1.06478,1.65831)
v(-1.06478,2.05894) e
bTriStrip
v(-1.90903,0.07696) v(-1.86303,0.326685) v(-1.74803,-0.15962)
v(-1.61002,0.530407) v(-1.54102,-0.448774) v(-1.38001,0.681555)
v(-1.31101,-0.685355) v(-1.01285,0.723487) v(-1.01353,-0.88907) e
bTriFan
v(-0.078033,-0.107463) v(-0.233835,0.218981) v(-0.129966,1.00542)
v(0.285509,0.753164) v(0.54518,0.441557) e
bTriFan
v(-0.493508,0.560265) v(-1.01285,0.723487) v(-0.545442,1.1538)
v(-0.129966,1.00542) v(-0.233835,0.218981) e
bTris
v(-0.545442,1.1538) v(-1.01285,0.723487) v(-1.035,1.16786)
v(-3.97311,0.545426) v(-5,5) v(-3.5057,0.9609) e
}
#undef v
#define v(x,y,z)  glVertex3f(x,y,z);
static void preTesselatedCutOutVolume(void) {
bQuadStrip
v(-5,-5,0) v(-5,-5,1) v(-5,5,0) v(-5,5,1) v(5,5,0) v(5,5,1) v(5,-5,0)
v(5,-5,1) v(-5,-5,0) v(-5,-5,1) e
bQuadStrip
v(-3.97311,0.545426,0) v(-3.97311,0.545426,1) v(-3.86924,0.204143,0)
v(-3.86924,0.204143,1) v(-3.60957,-0.315201,0) v(-3.60957,-0.315201,1)
v(-3.29797,-0.760353,0) v(-3.29797,-0.760353,1) v(-2.93443,-1.20551,0)
v(-2.93443,-1.20551,1) v(-2.41508,-1.63582,0) v(-2.41508,-1.63582,1)
v(-1.8438,-1.96226,0) v(-1.8438,-1.96226,1) v(-1.02549,-2.15912,0)
v(-1.02549,-2.15912,1) v(-1.02549,-2.81135,0) v(-1.02549,-2.81135,1)
v(3.97311,-2.81135,0) v(3.97311,-2.81135,1) v(3.97311,2.81135,0)
v(3.97311,2.81135,1) v(-1.06485,2.81135,0) v(-1.06485,2.81135,1)
v(-1.06478,2.05894,0) v(-1.06478,2.05894,1) v(-2.05154,1.85121,0)
v(-2.05154,1.85121,1) v(-2.83056,1.43573,0) v(-2.83056,1.43573,1)
v(-3.5057,0.9609,0) v(-3.5057,0.9609,1) v(-3.97311,0.545426,0)
v(-3.97311,0.545426,1) e
bQuadStrip
v(-2.57088,-0.760353,0) v(-2.57088,-0.760353,1) v(-2.83056,-0.315201,0)
v(-2.83056,-0.315201,1) v(-3.09023,0.055759,0) v(-3.09023,0.055759,1)
v(-3.1941,0.471234,0) v(-3.1941,0.471234,1) v(-2.83056,0.84219,0)
v(-2.83056,0.84219,1) v(-2.36315,1.18348,0) v(-2.36315,1.18348,1)
v(-1.79187,1.49508,0) v(-1.79187,1.49508,1) v(-1.06478,1.65831,0)
v(-1.06478,1.65831,1) v(-1.035,1.16786,0) v(-1.035,1.16786,1)
v(-1.63302,1.01671,0) v(-1.63302,1.01671,1) v(-2.04704,0.806418,0)
v(-2.04704,0.806418,1) v(-2.39205,0.536978,0) v(-2.39205,0.536978,1)
v(-2.53006,0.260968,0) v(-2.53006,0.260968,1) v(-2.41505,-0.047902,0)
v(-2.41505,-0.047902,1) v(-2.25405,-0.337056,0) v(-2.25405,-0.337056,1)
v(-2.04704,-0.606495,0) v(-2.04704,-0.606495,1) v(-1.79403,-0.875934,0)
v(-1.79403,-0.875934,1) v(-1.47202,-1.09937,0) v(-1.47202,-1.09937,1)
v(-1.02549,-1.24824,0) v(-1.02549,-1.24824,1) v(-1.02549,-1.73179,0)
v(-1.02549,-1.73179,1) v(-1.58413,-1.53195,0) v(-1.58413,-1.53195,1)
v(-2.10347,-1.22034,0) v(-2.10347,-1.22034,1) v(-2.57088,-0.760353,0)
v(-2.57088,-0.760353,1) e bQuadStrip v(3.34336,-0.719707,0)
v(3.34336,-0.719707,1) v(3.10721,-0.955861,0) v(3.10721,-0.955861,1)
v(2.67425,-1.21451,0) v(2.67425,-1.21451,1) v(2.20194,-1.45066,0)
v(2.20194,-1.45066,1) v(1.69027,-1.65308,0) v(1.69027,-1.65308,1)
v(1.17861,-1.833,0) v(1.17861,-1.833,1) v(0.706303,-1.97919,0)
v(0.706303,-1.97919,1) v(0.155277,-2.10289,0) v(0.155277,-2.10289,1)
v(-0.474465,-2.15912,0) v(-0.474465,-2.15912,1) v(-1.02549,-2.15912,0)
v(-1.02549,-2.15912,1) v(-1.02549,-1.73179,0) v(-1.02549,-1.73179,1)
v(-0.198952,-1.7093,0) v(-0.198952,-1.7093,1) v(0.273356,-1.5856,0)
v(0.273356,-1.5856,1) v(0.86374,-1.3382,0) v(0.86374,-1.3382,1)
v(1.29669,-1.07956,0) v(1.29669,-1.07956,1) v(1.69027,-0.820916,0)
v(1.69027,-0.820916,1) v(2.00515,-0.51729,0) v(2.00515,-0.51729,1)
v(2.4381,-0.179928,0) v(2.4381,-0.179928,1) v(3.34336,-0.719707,0)
v(3.34336,-0.719707,1) e
bQuadStrip
v(1.53284,1.12454,0) v(1.53284,1.12454,1) v(1.84771,0.83216,0)
v(1.84771,0.83216,1) v(2.12323,0.494796,0) v(2.12323,0.494796,1)
v(1.92643,0.269888,0) v(1.92643,0.269888,1) v(1.69027,-0.033738,0)
v(1.69027,-0.033738,1) v(1.29669,-0.348609,0) v(1.29669,-0.348609,1)
v(0.942455,-0.652235,0) v(0.942455,-0.652235,1) v(0.588227,-0.899634,0)
v(0.588227,-0.899634,1) v(0.115919,-1.1133,0) v(0.115919,-1.1133,1)
v(-0.474465,-1.25949,0) v(-0.474465,-1.25949,1) v(-1.02549,-1.24824,0)
v(-1.02549,-1.24824,1) v(-1.01353,-0.88907,0) v(-1.01353,-0.88907,1)
v(-0.441573,-0.87906,0) v(-0.441573,-0.87906,1) v(-0.078033,-0.745515,0)
v(-0.078033,-0.745515,1) v(0.389378,-0.5081,0) v(0.389378,-0.5081,1)
v(0.804852,-0.166817,0) v(0.804852,-0.166817,1) v(1.06452,0.189305,0)
v(1.06452,0.189305,1) v(1.3242,0.426719,0) v(1.3242,0.426719,1)
v(1.01259,0.753164,0) v(1.01259,0.753164,1) v(0.64905,1.09445,0)
v(0.64905,1.09445,1) v(0.233576,1.36154,0) v(0.233576,1.36154,1)
v(-0.233835,1.59895,0) v(-0.233835,1.59895,1) v(-1.06478,1.65831,0)
v(-1.06478,1.65831,1) v(-1.06478,2.05894,0) v(-1.06478,2.05894,1)
v(-0.395747,2.0804,0) v(-0.395747,2.0804,1) v(0.076562,1.96794,0)
v(0.076562,1.96794,1) v(0.706303,1.7093,0) v(0.706303,1.7093,1)
v(1.17861,1.39443,0) v(1.17861,1.39443,1) v(1.53284,1.12454,0)
v(1.53284,1.12454,1) e
bQuadStrip
v(-1.01285,0.723487,0) v(-1.01285,0.723487,1) v(-1.01353,-0.88907,0)
v(-1.01353,-0.88907,1) v(-1.31101,-0.685355,0) v(-1.31101,-0.685355,1)
v(-1.54102,-0.448774,0) v(-1.54102,-0.448774,1) v(-1.74803,-0.15962,0)
v(-1.74803,-0.15962,1) v(-1.90903,0.07696,0) v(-1.90903,0.07696,1)
v(-1.86303,0.326685,0) v(-1.86303,0.326685,1) v(-1.61002,0.530407,0)
v(-1.61002,0.530407,1) v(-1.38001,0.681555,0) v(-1.38001,0.681555,1)
v(-1.01285,0.723487,0) v(-1.01285,0.723487,1) e
bQuadStrip
v(0.54518,0.441557,0) v(0.54518,0.441557,1) v(-0.078033,-0.107463,0)
v(-0.078033,-0.107463,1) v(-0.233835,0.218981,0) v(-0.233835,0.218981,1)
v(-0.493508,0.560265,0) v(-0.493508,0.560265,1) v(-1.01285,0.723487,0)
v(-1.01285,0.723487,1) v(-1.035,1.16786,0) v(-1.035,1.16786,1)
v(-0.545442,1.1538,0) v(-0.545442,1.1538,1) v(-0.129966,1.00542,0)
v(-0.129966,1.00542,1) v(0.285509,0.753164,0) v(0.285509,0.753164,1)
v(0.54518,0.441557,0) v(0.54518,0.441557,1) e
}

#undef v
#undef bTriStrip
#undef bTriFan
#undef bTris
#undef bQuadStrip
#endif  /* !defined(GLU_VERSION_1_2) */

void
doContour(GLUtesselator *tess, int n, GLdouble o[][3])
{
  int i;

  gluTessBeginContour(tess);
  for (i=0; i< n; i++) {
    gluTessVertex(tess, o[i], o[i]);
  }
  gluTessEndContour(tess);
}

void
doCutOut(void)
{
  gluTessBeginPolygon(tess, NULL);
    doContour(tess, SIZE(base), base);
    doContour(tess, SIZE(nvlogo0), nvlogo0);
    doContour(tess, SIZE(nvlogo1), nvlogo1);
    doContour(tess, SIZE(nvlogo2), nvlogo2);
    doContour(tess, SIZE(nvlogo3), nvlogo3);
  gluTessEndPolygon(tess);
}

static int first;
static GLfloat fv[2];

static void CALLBACK
begin(GLenum type)
{
  glBegin(type);
}

/* ARGSUSED */
static void CALLBACK
bside(GLenum type)
{
  first = 1;
  glBegin(GL_QUAD_STRIP);
}

static void CALLBACK
eside(void)
{
  GLfloat v[3];

  v[0] = fv[0];
  v[1] = fv[1];
  v[2] = 0.0;
  glVertex3fv(v);
  v[2] = 1.0;
  glVertex3fv(v);
  glEnd();
}

static void CALLBACK
vside(void *data)
{
  GLdouble *d = (GLdouble*) data;
  GLfloat v[3];

  if (first) {
    fv[0] = d[0];
    fv[1] = d[1];
    first = 0;
  }

  v[0] = d[0];
  v[1] = d[1];
  v[2] = 0.0;
  glVertex3fv(v);
  v[2] = 1.0;
  glVertex3fv(v);
}

static void CALLBACK
end(void)
{
   glEnd();
}

void
initCutOut(void)
{
#ifdef GLU_VERSION_1_2
  gluTessProperty(tess, GLU_TESS_BOUNDARY_ONLY, GL_FALSE);
  gluTessProperty(tess, GLU_TESS_WINDING_RULE, GLU_TESS_WINDING_ODD);
  gluTessCallback(tess, GLU_TESS_BEGIN, (void (CALLBACK*)()) &begin);
  gluTessCallback(tess, GLU_TESS_VERTEX, (void (CALLBACK*)()) &glVertex3dv);
  gluTessCallback(tess, GLU_TESS_END, (void (CALLBACK*)()) &end);

  glNewList(DL_CUT_OUT, GL_COMPILE);
    doCutOut();
  glEndList();
#else
  /* When GLU 1.2 is not supported. */
  glNewList(DL_CUT_OUT, GL_COMPILE);
    preTesselatedCutOut();
  glEndList();
#endif
}

void
drawCutOut(void)
{
  glNormal3f(0,0,1);
  glCallList(1);
  glNormal3f(0,0,-1);
  glCallList(1);
}

void CALLBACK
combine(GLdouble coords[3], GLdouble *data[4],
        GLfloat weight[4], GLdouble **dataOut )
{
   GLdouble *vertex;
   vertex = (GLdouble *) malloc(3 * sizeof(GLdouble));

   vertex[0] = coords[0];
   vertex[1] = coords[1];
   vertex[2] = coords[2];
   *dataOut = vertex;
}

static void CALLBACK
error(GLenum errno)
{
  fprintf(stderr, "ERROR: %s\n", gluErrorString(errno));
}

GLUtesselator *
initTess(void)
{
  GLUtesselator *tess;

  tess = gluNewTess();
  if (tess == NULL) {
    return NULL;
  }
  gluTessCallback(tess, GLU_TESS_BEGIN, (void (CALLBACK*)()) &begin);
  gluTessCallback(tess, GLU_TESS_VERTEX, (void (CALLBACK*)()) &glVertex3dv);
  gluTessCallback(tess, GLU_TESS_COMBINE, (void (CALLBACK*)()) &combine);
  gluTessCallback(tess, GLU_TESS_END, (void (CALLBACK*)()) &end);
  gluTessCallback(tess, GLU_TESS_ERROR, (void (CALLBACK*)()) &error);
  return tess;
}

void
initCutOutVolume(void)
{
#ifdef GLU_VERSION_1_2
  gluTessProperty(tess, GLU_TESS_BOUNDARY_ONLY, GL_TRUE);
  gluTessProperty(tess, GLU_TESS_WINDING_RULE, GLU_TESS_WINDING_ODD);
  gluTessCallback(tess, GLU_TESS_BEGIN, (void (CALLBACK*)()) &bside);
  gluTessCallback(tess, GLU_TESS_VERTEX, (void (CALLBACK*)()) &vside);
  gluTessCallback(tess, GLU_TESS_END, (void (CALLBACK*)()) &eside);

  glNewList(DL_CUT_OUT_VOLUME, GL_COMPILE);
    doCutOut();
  glEndList();
#else
  /* When GLU 1.2 is not supported. */
  glNewList(DL_CUT_OUT_VOLUME, GL_COMPILE);
    preTesselatedCutOutVolume();
  glEndList();
#endif
}

void
drawCutOutVolume(void)
{
  /* Note that this routine does not cap the ends of the shadow volume.
     The cut-out object itself and a ground plane are assumed to cap
     the shadow volume. */
  glCallList(DL_CUT_OUT_VOLUME);
}

void
drawFloor(void)
{
  float x, y;
  float d = 1.0;

  /* Draw ground. */
  if (useTextures) {
    glEnable(GL_TEXTURE_2D);
  }
  glColor3f(1,1,0);
  glNormal3f(0,1,0);
  /* Tesselate floor so lighting looks reasonable. */
  for (x=-8; x<=7; x+=d) {
    glBegin(GL_QUAD_STRIP);
    for (y=-8; y<=8; y+=d) {
      glTexCoord2f(x+1, y);
      glVertex3f(x+1, -1, y);
      glTexCoord2f(x, y);
      glVertex3f(x, -1, y);
    }
    glEnd();
  }
  if (useTextures) {
    glDisable(GL_TEXTURE_2D);
  }
}

void
initObjects(void)
{
  glNewList(DL_SPHERE, GL_COMPILE);
    glutSolidSphere(0.5, 7, 7);
  glEndList();

  glNewList(DL_TEAPOT, GL_COMPILE);
    glDisable(GL_CULL_FACE);
    fastTeapot(3, 1.6);
    glEnable(GL_CULL_FACE);
  glEndList();
}

void
drawObjects(void)
{
  float h;

  glPushMatrix();
    h = (-(1-time1)*(1-time1) + 1) * 1.5;
    glTranslatef(0,-0.5+h, 0);
    glCallList(DL_SPHERE);
  glPopMatrix();

  glPushMatrix();
    glTranslatef(2,-0.5,2+sin(time2));
    glCallList(DL_SPHERE);
  glPopMatrix();

  glPushMatrix();
    glTranslatef(2+cos(time2),-0.5,-2+0.4*sin(time2));
    glCallList(DL_SPHERE);
  glPopMatrix();

  glPushMatrix();
    glTranslatef(-2,-0.0,2);
    glRotatef(180, 0,1,0);
    glRotatef(teapotAngle, 0,1,0);
    glCallList(DL_TEAPOT);
  glPopMatrix();
}

void
drawRoom(void)
{
  drawFloor();
  drawObjects();
}

void
ambientLight(int enable)
{
  GLfloat on[4] = { 0.92, 0.92, 0.92, 1.0 };
  GLfloat off[4] = { 0.0, 0.0, 0.0, 1.0 };

  if (enable) {
    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, on);
  } else {
    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, off);
  }
}

void
renderWithSoftShadows(int numLightSamples)
{
  GLdouble P[5][3];
  GLdouble matrix[4][4];
  GLfloat Lf[4];
  GLdouble Ld[4];
  GLfloat Lc[4];
  float sampleAngle, sampleDeltaAngle;
  float lightRadius = 0.1;
  int doOnce;
  int i;

  glDepthMask(1);
  glColorMask(1,1,1,1);
  glStencilMask(~0u);
  /* The subsquent shadow passes use additive blending and to prevent the
     possiblity of "double blending" due to "depth testing ties", stencil
     is used to accept the first of any depth testing ties.  Therefore,
     use GL_LESS initially since it too accepts just the first of any
     depth testing ties. */
  glDepthFunc(GL_LESS);

  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);
  
  /* The samples will be in a ring in a constant Z plane. */
  sampleDeltaAngle = 3.14159 * 2.0 / numLightSamples;

  /* Draw light source position as a yellow dot. */
  glDisable(GL_LIGHTING);
  glColor3f(1,1,0);
  glBegin(GL_POINTS);
  for (i=0, sampleAngle=0.0;
       i<numLightSamples;
       i++, sampleAngle += sampleDeltaAngle) {
    glVertex3d(
      L[0] + lightRadius * cos(sampleAngle),
      L[1],
      L[2] + lightRadius * sin(sampleAngle));
  }
  glEnd();

  /* Draw room with no lights sources on, but still have ambient lighting. */
  glEnable(GL_LIGHTING);
  glDisable(GL_LIGHT0);
  ambientLight(1);
  /* Use back face culling for better performance. */
  /* XXX Too bad that crappy glutSolidTeapot model is not a true closed
     surface.  You may notice minor cracks in the teapot when backface
     culling is enabled. */
  glEnable(GL_CULL_FACE);
  glCullFace(GL_BACK);
  if (forceStencilHack) {
    glEnable(GL_STENCIL_TEST);
    glStencilFunc(GL_ALWAYS, 0, 0xffffffff);
    glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
  } else {
    /* XXX Note that the non-hack case does not enable or
       disable stencil.  The hack case enables stencil
       testing but sets up the stencil modes so that stencil
       testing is effectively disabled.  If you wanted
       stencil testing on during the renderSceneFunc, you won't
       need to have the hack enabled though! */
  }
  drawRoom();

  /* Set the light source samples to be 1/nth of the strength of the true
     area light source. */
  Lc[0] = 1.0 / numLightSamples;
  Lc[1] = 1.0 / numLightSamples;
  Lc[2] = 1.0 / numLightSamples;
  Lc[3] = 1.0 / numLightSamples;
  glLightfv(GL_LIGHT0, GL_DIFFUSE, &Lc[0]);
  glLightfv(GL_LIGHT0, GL_SPECULAR, &Lc[0]);
  glEnable(GL_LIGHT0);

  doOnce = 1;

  for (i=0, sampleAngle=0.0;
       i<numLightSamples;
       i++, sampleAngle += sampleDeltaAngle) {

    /* Compute each light sample position. */
    Lf[0] = Ld[0] = L[0] + lightRadius * cos(sampleAngle);
    Lf[1] = Ld[1] = L[1] + lightRadius * sin(sampleAngle);
    Lf[2] = Ld[2] = L[2];
    Lf[3] = Ld[3] = L[3];

    /* Tell OpenGL the light sample's position. */
    glLightfv(GL_LIGHT0, GL_POSITION, &Lf[0]);

    /* Given: the light source position (L), the cut-out origin (O), the
       cut-out Y direction (N), the cut-out X direction (M), the ground
       plane equation, compute five points in world space to project
       the cut-out space volume into. */
    computeShadowVolume(Ld, O, M, N, Pbg, P);
  
    /* Generate a cut-out shadow volume projection matrix from the five
       points in world space. */
    constructShadowVolumeMatrixd(P, matrix);
  
    glPushMatrix();
      /* Concatenate the previously computed 4x4 projective matrix
         onto the current modelview matrix.  With this matrix in place,
         the "drawCutOut" routine will draw the cut-out in the proper
         orientation as described by the parameters to the call to
         computeShadowVolume and the "drawCutOutVolume" routine will
         draw the cut-out's shadow volume extended to the ground plane as
         described by the parameters to the computeShadowVolume call. */
      glMultMatrixd(&matrix[0][0]);

      glDepthFunc(GL_ALWAYS);  /* XXX hack */
      drawCutOut();

      if (doOnce) {
        /* Disable the global ambient light for subsequent passes. */
        ambientLight(0);
        /* Enable additive blending to add in each light source sample's
           lighitng contribution. */
        glEnable(GL_BLEND);
        glBlendFunc(GL_ONE, GL_ONE);

        glDepthMask(0);
        doOnce = 0;
      }

      /* Draw the shadow volume. */
      glEnable(GL_STENCIL_TEST);
      glDepthFunc(GL_LESS);
      glColorMask(0,0,0,0);
      glStencilOp(GL_KEEP, GL_KEEP, GL_INVERT);
      glStencilFunc(GL_ALWAYS, 0, 1 << i);
      glStencilMask(1 << i);
      glDisable(GL_CULL_FACE);
      drawCutOutVolume();
      glEnable(GL_CULL_FACE); 
    glPopMatrix();

    glColorMask(1,1,1,1);
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
    glStencilFunc(GL_NOTEQUAL, 1 << i, 1 << i);
    glDepthFunc(GL_EQUAL);
    /* Draw the room with the ith light sample enabled. */
    drawRoom();

    if (forceStencilHack) {
      glEnable(GL_STENCIL_TEST);
      glStencilFunc(GL_ALWAYS, 0, 0xffffffff);
      glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
    } else {
      glDisable(GL_STENCIL_TEST);
    }
  }

  glDisable(GL_BLEND);
}

void
reshape(int w, int h)
{
  windowWidth = w;
  windowHeight = h;
  glViewport(0,0,w,h);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(60.0, (GLdouble)w/h, 1.0, 20.0);
  glPushMatrix();
  glMatrixMode(GL_MODELVIEW);  
}

void
render(void)
{
  GLdouble P[5][3];
  GLdouble matrix[4][4];
  GLfloat shadowMat[4][4];
  GLfloat Lf[4];
  GLfloat Lc[4];

  /* Use a white light. */
  Lc[0] = 1.0;
  Lc[1] = 1.0;
  Lc[2] = 1.0;
  Lc[3] = 1.0;
  glLightfv(GL_LIGHT0, GL_DIFFUSE, &Lc[0]);
  glLightfv(GL_LIGHT0, GL_SPECULAR, &Lc[0]);

  glColorMask(1,1,1,1);
  glDepthMask(1);
  glStencilMask(~0u);
  /* GL_LESS means that depth test "ties" go to the FIRST tieing fragment.
     This is OpenGL's default depth test, but when depth testing with
     shadow volumes, we later use GL_EQUAL to match depth values.  However,
     GL_EQUAL means that depth test "ties" go to the LAST tieing fragment.
     Therefore, we use GL_LEQUAL initially so we can match the later
     "depth test tie goes to the last fragment" passes. */
  glDepthFunc(GL_LEQUAL);
  if (forceStencilHack) {
    glEnable(GL_STENCIL_TEST);
    glStencilFunc(GL_ALWAYS, 0, 0xffffffff);
    glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
  } else {
    glDisable(GL_STENCIL_TEST);
  }
  glEnable(GL_LIGHT0);

  ambientLight(1);

  Lf[0] = L[0];
  Lf[1] = L[1];
  Lf[2] = L[2];
  Lf[3] = L[3];
  glLightfv(GL_LIGHT0, GL_POSITION, &Lf[0]);
  
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);
  
  /* Given: the light source position (L), the cut-out origin (O), the
     cut-out Y direction (N), the cut-out X direction (M), the ground
     plane equation, compute five points in world space to project the
     cut-out space volume into. */
  computeShadowVolume(L, O, M, N, Pbg, P);
  
  /* Generate a volume projection matrix from the five points in world
     space. */
  constructShadowVolumeMatrixd(P, matrix);

  shadowMatrix(shadowMat, Pg, Lf);
   
  glDisable(GL_LIGHTING);
  /* Draw light source position as a yellow dot. */
  glColor3f(1,1,0);
  glBegin(GL_POINTS);
  glVertex3dv(L);
  glEnd();
  glEnable(GL_LIGHTING);

  glEnable(GL_CULL_FACE);
  glCullFace(GL_BACK);
  drawObjects();

  glPushMatrix();
    glMultMatrixd(&matrix[0][0]);
    drawCutOut();
  glPopMatrix();

  if (displayMode == M_PLANAR_SHADOWS || displayMode == M_HYBRID_SHADOWS) {
    glEnable(GL_STENCIL_TEST);
    glStencilFunc(GL_ALWAYS, 0x2, 0);
    glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
  }
  drawFloor();

  if (displayMode != M_NO_SHADOWS) {
    glPushMatrix();
    glMultMatrixd(&matrix[0][0]);
    if (displayMode == M_VOLUME_SHADOWS || displayMode == M_HYBRID_SHADOWS) {
      
      /* Draw the shadow volume. */
      glColorMask(0,0,0,0);
      glDepthMask(0);
      glEnable(GL_STENCIL_TEST);
      glStencilOp(GL_KEEP, GL_KEEP, GL_INVERT);
      glStencilFunc(GL_ALWAYS, 0, 0x1);
      glStencilMask(0x1);
      glDisable(GL_CULL_FACE);
      drawCutOutVolume();
      glEnable(GL_CULL_FACE); 
    }
    if (displayMode == M_SHADOW_VOLUME) {
      glDisable(GL_LIGHTING);
      glEnable(GL_BLEND);
      glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
      glDepthMask(0);

      glColor4f(1,0,1,0.8);
      glColorMask(1,0,0,0);
      glCullFace(GL_FRONT);
      drawCutOutVolume();

      glColorMask(0,0,1,0);
      glColor4f(1,0,1,0.8);
      glCullFace(GL_BACK);
      drawCutOutVolume();

      glDepthMask(1);
      glDisable(GL_BLEND);
      glEnable(GL_LIGHTING);
      glColorMask(1,1,1,1);
    }
    glPopMatrix();
  }

  if (displayMode == M_VOLUME_SHADOWS || displayMode == M_HYBRID_SHADOWS) {
    glColorMask(1,1,1,1);
    glStencilOp(GL_KEEP, GL_KEEP, GL_ZERO);
    glStencilFunc(GL_EQUAL, 0x1, 0x1);
    glDepthFunc(GL_EQUAL);
    glDisable(GL_LIGHT0);
    
    drawRoom();
    
    glEnable(GL_LIGHT0);
  }

  if (displayMode == M_PLANAR_SHADOWS || displayMode == M_HYBRID_SHADOWS) {
    glDisable(GL_DEPTH_TEST);

    glPushMatrix();
      glMultMatrixf((GLfloat *) shadowMat);
      glStencilFunc(GL_EQUAL, 0x3, 0x2);
      glStencilOp(GL_KEEP, GL_KEEP, GL_REPLACE);
      drawObjects();

      if (displayMode == M_PLANAR_SHADOWS) {
        glPushMatrix();
          glMultMatrixd(&matrix[0][0]);
          drawCutOut();
        glPopMatrix();
      }
    glPopMatrix();

    glStencilFunc(GL_EQUAL, 0x1, 0x1);
    glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
    glColorMask(1,1,1,1);
    glDisable(GL_LIGHT0);
    drawFloor();
    glEnable(GL_DEPTH_TEST); 
  }

  if (forceStencilHack) {
    glEnable(GL_STENCIL_TEST);
    glStencilFunc(GL_ALWAYS, 0, 0xffffffff);
    glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);
  } else {
    glDisable(GL_STENCIL_TEST);
  }
}

void
display(void)
{
  if (displayMode == M_SOFT_SHADOWS) {
    /* Handle with a special routine. */
    renderWithSoftShadows(numLightSamples);
  } else {
    /* Handle all the non-soft shadow cases. */
    render();
  }
  if (doubleBuffer) {
    glutSwapBuffers();
  } else {
    glFlush();
  }
  if (oneFrame) {
    exit(0);
  }
}

void
updateWindowTitle(void)
{
  char title[100];

  switch (displayMode) {
  case M_VOLUME_SHADOWS:
    glutSetWindowTitle("nvshadow: Hard Volmetric Shadows");
    break;
  case M_SOFT_SHADOWS:
    sprintf(title, "nvshadow: Soft Shadows (%d samples)",
      numLightSamples);
    glutSetWindowTitle(title);
    break;
  case M_NO_SHADOWS:
    glutSetWindowTitle("nvshadow: No Shadows");
    break;
  case M_SHADOW_VOLUME:
    glutSetWindowTitle("nvshadow: Visualize Shadow Volume");
    break;
  case M_PLANAR_SHADOWS:
    glutSetWindowTitle("nvshadow: Planar Projected Shadows");
    break;
  case M_HYBRID_SHADOWS:
    glutSetWindowTitle("nvshadow: Hybrid Shadows");
    break;
  }
}

void
updateAxis(void)
{
  M[0] = 0.3 * cos(logoAngle/180.0*M_PI);
  M[1] = 0.3 * 0.0;
  M[2] = 0.3 * -sin(logoAngle/180.0*M_PI);

  N[0] = 0.3 * sin(logoAngle/180.0*M_PI);
  N[1] = 0.3 * 0.0;
  N[2] = 0.3 * cos(logoAngle/180.0*M_PI);
}

void
benchmark(void)
{
  float start, end;
  int frames = 25;
  int i;

  L[0] = 3;
  L[1] = 4;
  L[2] = 2;
  L[3] = 1;

  O[0] = 3;
  O[1] = 2;
  O[2] = 2;

  updateAxis();

  start = glutGet(GLUT_ELAPSED_TIME);
  for (i=0; i<frames; i++) {
    L[0] = L[0] + 0.15;
    display();
  }
  end = glutGet(GLUT_ELAPSED_TIME);
  printf("\nseconds = %g\n", (end - start)/1000.0);
  printf("seconds/frame = %.3g\n", ((end - start)/1000.0)/frames);
  printf("frames/second = %.2g\n", frames/((end - start)/1000.0));
  printf("size = %dx%d\n",
    glutGet(GLUT_WINDOW_WIDTH), glutGet(GLUT_WINDOW_HEIGHT));
}

void
doAnimate(void)
{
  if (animateMode != 2) {
    teapotAngle = (teapotAngle + 5) % 360;
    time1 += 0.1;
    if (time1 > 2.0) {
      time1 = time1 - 2.0;
    }
    time2 += 0.05;
  }
  if (animateMode != 0) {
    logoAngle = (logoAngle + 3) % 360;
    updateAxis();
  }
  glutPostRedisplay();
}

static void
visible(int vis)
{
  if (vis == GLUT_VISIBLE) {
    if (animate) 
      glutIdleFunc(doAnimate);
  } else {
    if (!animate)
      glutIdleFunc(NULL);
  }
}

static void
keyboard(unsigned char c, int x, int y)
{
  switch (c) {
  case 27:
    exit(0);
    return;
  case 'h':
  case 'H':
    O[0] -= 0.25;
    break;
  case 'l':
  case 'L':
    O[0] += 0.25;
    break;
  case 'j':
  case 'J':
    O[2] -= 0.25;
    break;
  case 'k':
  case 'K':
    O[2] += 0.25;
    break;
  case 'a':
  case 'A':
    O[1] += 0.25;
    break;
  case 'z':
  case 'Z':
    O[1] -= 0.25;
    break;
  case '2':
  case '3':
  case '4':
  case '5':
  case '6':
  case '7':
  case '8':
    numLightSamples = c - '0';
    updateWindowTitle();
    break;
  case 'b':
  case 'B':
    benchmark();
    break;
  case 'm':
    animateMode = (animateMode+1) % 3;
    break;
  case ' ':
    animate = !animate;
    if (animate) {
      glutIdleFunc(doAnimate);
    } else {
      glutIdleFunc(NULL);
    }
    break;
  default:
    return;
  }
  glutPostRedisplay();
}

void
special(int k, int x, int y)
{
  switch (k) {
  case GLUT_KEY_F1:
    displayMode = M_VOLUME_SHADOWS;
    updateWindowTitle();
    break;
  case GLUT_KEY_F2:
    displayMode = M_NO_SHADOWS;
    updateWindowTitle();
    break;
  case GLUT_KEY_F3:
    displayMode = M_SHADOW_VOLUME;
    updateWindowTitle();
    break;
  case GLUT_KEY_F4:
    displayMode = M_SOFT_SHADOWS;
    updateWindowTitle();
    break;
  case GLUT_KEY_F5:
    displayMode = M_PLANAR_SHADOWS;
    updateWindowTitle();
    break;
  case GLUT_KEY_F6:
    displayMode = M_HYBRID_SHADOWS;
    updateWindowTitle();
    break;
  case GLUT_KEY_UP:
    L[2] -= 0.25;
    break;
  case GLUT_KEY_DOWN:
    L[2] += 0.25;
    break;
  case GLUT_KEY_RIGHT:
    L[0] += 0.25;
    break;
  case GLUT_KEY_LEFT:
    L[0] -= 0.25;
    break;
  case GLUT_KEY_PAGE_UP:
    L[1] += 0.25;
    break;
  case GLUT_KEY_PAGE_DOWN:
    L[1] -= 0.25;
    break;
#if 0  /* If you wanted keys to raise and lower the virtual ground plane. */
  case GLUT_KEY_HOME:
    Pbg[3] += 0.25;
    break;
  case GLUT_KEY_END:
    Pbg[3] -= 0.25;
    break;
#endif
  default:
    /* Avoid the glutPostRedisplay below. */
    return;
  }
  glutPostRedisplay();
}

static void
mouse(int button, int state, int x, int y)
{
  if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
    moving = 1;
    beginx = x;
    beginy = y;
  }
  if (button == GLUT_LEFT_BUTTON && state == GLUT_UP) {
    moving = 0;
  }
}

static void
motion(int x, int y)
{
  if (moving) {
    logoAngle = logoAngle + (beginx - x);

    updateAxis();

    beginx = x;
    beginy = y;
    glutPostRedisplay();
  }
}

static void
menu(int item)

{
  switch (item) {
  case M_ANIMATE:
    keyboard(' ', 0, 0);
    break;
  case M_BENCHMARK:
    keyboard('b', 0, 0);
    break;
  case M_VOLUME_SHADOWS:
  case M_SOFT_SHADOWS:
  case M_NO_SHADOWS:
  case M_SHADOW_VOLUME:
  case M_PLANAR_SHADOWS:
  case M_HYBRID_SHADOWS:
    displayMode = item;
    updateWindowTitle();
    break;
  }
  glutPostRedisplay();
}

static void
numLightSamplesMenu(int samples)
{
  numLightSamples = samples;
  updateWindowTitle();
  glutPostRedisplay();
}

/* XXX RIVA 128 board vendors may change their GL_VENDOR
   and GL_RENDERER strings. */
int
needsStencilRenderingInvariantHack(void)
{
  const char *renderer;
  GLint bits;

  renderer = glGetString(GL_RENDERER);
  /* Stencil rendering on RIVA 128 and RIVA 128 ZX 
     is not invariant with stencil-disabled rendering
     in 16-bit hardware accelerated mode. */
  if (!strncmp("RIVA 128", renderer, 8)) {
    glGetIntegerv(GL_INDEX_BITS, &bits);
    return bits == 16;
  }
  /* Stencil rendering on RIVA 128 and RIVA 128 ZX 
     is not invariant with stencil-disabled rendering
     in 16-bit hardware accelerated mode.  32-bit mode
     is invariant (and hardware accelerated though!).  */
  if (!strncmp("RIVA TNT", renderer, 8)) {
    glGetIntegerv(GL_INDEX_BITS, &bits);
    return bits == 16;
  }
  return 1;
}

int
main(int argc, char **argv)
{
  int i;
  int sub1;

  glutInit(&argc, argv);

  for (i=1; i<argc; i++) {
    if(!strcmp("-notex", argv[i])) {
      useTextures = 0;
    }
    else
    if(!strcmp("-stencilhack", argv[i])) {
      forceStencilHack = 1;
    }
    else
    if(!strcmp("-fullscreen", argv[i])) {
      fullscreen = 1;
    }
    else
    if(!strcmp("-400x300", argv[i])) {
      fullscreen = 1;
      smallFullscreen = 1;
    }
    else
    if(!strcmp("-soft", argv[i])) {
      displayMode = M_SOFT_SHADOWS;
    }
    else
    if(!strcmp("-volume", argv[i])) {
      displayMode = M_SHADOW_VOLUME;
    }
    else
    if(!strcmp("-planar", argv[i])) {
      displayMode = M_PLANAR_SHADOWS;
    }
    else
    if(!strcmp("-noshadow", argv[i])) {
      displayMode = M_NO_SHADOWS;
    }
    else
    if(!strcmp("-nomipmaps", argv[i])) {
      useMipmaps = 0;
    }
    else
    if(!strcmp("-oneframe", argv[i])) {
      oneFrame = 1;
    }
    else
    if(!strcmp("-cool", argv[i])) {
      time1 = 0.5;
      time2 = 2.42;
      teapotAngle = 65;
      logoAngle = 78;
      updateAxis();
    }
    else
    if(!strcmp("-sb", argv[i])) {
      doubleBuffer = 0;
    }
  }

  if (doubleBuffer) {
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH | GLUT_STENCIL);
  } else {
    glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE | GLUT_DEPTH | GLUT_STENCIL);
  }
  if (fullscreen) {
    if (smallFullscreen) {
      glutGameModeString("400x300:32");
    } else {
      glutGameModeString("640x480:32");
    }
    glutEnterGameMode();
    glutSetCursor(GLUT_CURSOR_NONE);
  } else {
    glutCreateWindow("nvshadow");
    updateWindowTitle();
  }
  glutReshapeFunc(reshape);
  glutDisplayFunc(display);
  glutVisibilityFunc(visible);

  glutSpecialFunc(special);
  glutKeyboardFunc(keyboard);
  glutMouseFunc(mouse);
  glutMotionFunc(motion);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective(60.0, 1.0, 1.0, 20.0);

  glMatrixMode(GL_MODELVIEW);
  glLoadIdentity();
  gluLookAt(
    0,5,10,          /* eye */
    0,0,0,           /* center */
    0,1,0);          /* up */

  glEnable(GL_DEPTH_TEST);
  glPointSize(4.0);

  /* Determine the plane equation of the ground plane from three points. */
  findPlane(Pbg, Sbg, Tbg, Rbg);
  findPlane(Pg, Sg, Tg, Rg);

  if (useTextures) {
    makeFloorTexture();
  }

  glEnable(GL_LIGHTING);
  glEnable(GL_LIGHT0);

  tess = initTess();

  initCutOut();
  initCutOutVolume();
  initObjects();

  if (!fullscreen) {
    sub1 = glutCreateMenu(numLightSamplesMenu);
    glutAddMenuEntry("2", 2);
    glutAddMenuEntry("3", 3);
    glutAddMenuEntry("4", 4);
    glutAddMenuEntry("5", 5);
    glutAddMenuEntry("6", 6);
    glutAddMenuEntry("7", 7);
    glutAddMenuEntry("8", 8);

    glutCreateMenu(menu);
    glutAddMenuEntry("Toggle animation", M_ANIMATE);
    glutAddMenuEntry("Volumetric shadows", M_VOLUME_SHADOWS);
    glutAddMenuEntry("Soft shadows", M_SOFT_SHADOWS);
    glutAddMenuEntry("No shadows", M_NO_SHADOWS);
    glutAddMenuEntry("Shadow volumes", M_SHADOW_VOLUME);
    glutAddMenuEntry("Planar shadows", M_PLANAR_SHADOWS);
    glutAddMenuEntry("Hybrid shadows", M_HYBRID_SHADOWS);
    glutAddSubMenu("Light samples", sub1);
    glutAddMenuEntry("Do benchmark", M_BENCHMARK);
    glutAttachMenu(GLUT_RIGHT_BUTTON);
  }

  if (needsStencilRenderingInvariantHack()) {
    forceStencilHack = 1;
  }
  if (forceStencilHack) {
    printf("nvshadow: using stencil hack\n");
  }

  /* Scale the orientations for the cut-out. */
  M[0] = 0.3;
  M[1] = 0.0;
  M[2] = 0.0;

  N[0] = 0.0;
  N[1] = 0.0;
  N[2] = 0.3;

  glutMainLoop();
  return 0;
}

/* Rim, body, lid, and bottom data must be reflected in x and
   y; handle and spout data across the y axis only.  */

static int patchdata[][16] =
{
    /* rim */
  {102, 103, 104, 105, 4, 5, 6, 7, 8, 9, 10, 11,
    12, 13, 14, 15},
    /* body */
  {12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
    24, 25, 26, 27},
  {24, 25, 26, 27, 29, 30, 31, 32, 33, 34, 35, 36,
    37, 38, 39, 40},
    /* lid */
  {96, 96, 96, 96, 97, 98, 99, 100, 101, 101, 101,
    101, 0, 1, 2, 3,},
  {0, 1, 2, 3, 106, 107, 108, 109, 110, 111, 112,
    113, 114, 115, 116, 117},
    /* bottom */
  {118, 118, 118, 118, 124, 122, 119, 121, 123, 126,
    125, 120, 40, 39, 38, 37},
    /* handle */
  {41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52,
    53, 54, 55, 56},
  {53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64,
    28, 65, 66, 67},
    /* spout */
  {68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
    80, 81, 82, 83},
  {80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91,
    92, 93, 94, 95}
};
/* *INDENT-OFF* */

static float cpdata[][3] =
{
    {0.2, 0, 2.7}, {0.2, -0.112, 2.7}, {0.112, -0.2, 2.7}, {0,
    -0.2, 2.7}, {1.3375, 0, 2.53125}, {1.3375, -0.749, 2.53125},
    {0.749, -1.3375, 2.53125}, {0, -1.3375, 2.53125}, {1.4375,
    0, 2.53125}, {1.4375, -0.805, 2.53125}, {0.805, -1.4375,

    2.53125}, {0, -1.4375, 2.53125}, {1.5, 0, 2.4}, {1.5, -0.84,
    2.4}, {0.84, -1.5, 2.4}, {0, -1.5, 2.4}, {1.75, 0, 1.875},
    {1.75, -0.98, 1.875}, {0.98, -1.75, 1.875}, {0, -1.75,
    1.875}, {2, 0, 1.35}, {2, -1.12, 1.35}, {1.12, -2, 1.35},
    {0, -2, 1.35}, {2, 0, 0.9}, {2, -1.12, 0.9}, {1.12, -2,
    0.9}, {0, -2, 0.9}, {-2, 0, 0.9}, {2, 0, 0.45}, {2, -1.12,
    0.45}, {1.12, -2, 0.45}, {0, -2, 0.45}, {1.5, 0, 0.225},
    {1.5, -0.84, 0.225}, {0.84, -1.5, 0.225}, {0, -1.5, 0.225},
    {1.5, 0, 0.15}, {1.5, -0.84, 0.15}, {0.84, -1.5, 0.15}, {0,
    -1.5, 0.15}, {-1.6, 0, 2.025}, {-1.6, -0.3, 2.025}, {-1.5,
    -0.3, 2.25}, {-1.5, 0, 2.25}, {-2.3, 0, 2.025}, {-2.3, -0.3,
    2.025}, {-2.5, -0.3, 2.25}, {-2.5, 0, 2.25}, {-2.7, 0,
    2.025}, {-2.7, -0.3, 2.025}, {-3, -0.3, 2.25}, {-3, 0,
    2.25}, {-2.7, 0, 1.8}, {-2.7, -0.3, 1.8}, {-3, -0.3, 1.8},
    {-3, 0, 1.8}, {-2.7, 0, 1.575}, {-2.7, -0.3, 1.575}, {-3,
    -0.3, 1.35}, {-3, 0, 1.35}, {-2.5, 0, 1.125}, {-2.5, -0.3,
    1.125}, {-2.65, -0.3, 0.9375}, {-2.65, 0, 0.9375}, {-2,
    -0.3, 0.9}, {-1.9, -0.3, 0.6}, {-1.9, 0, 0.6}, {1.7, 0,
    1.425}, {1.7, -0.66, 1.425}, {1.7, -0.66, 0.6}, {1.7, 0,
    0.6}, {2.6, 0, 1.425}, {2.6, -0.66, 1.425}, {3.1, -0.66,
    0.825}, {3.1, 0, 0.825}, {2.3, 0, 2.1}, {2.3, -0.25, 2.1},
    {2.4, -0.25, 2.025}, {2.4, 0, 2.025}, {2.7, 0, 2.4}, {2.7,
    -0.25, 2.4}, {3.3, -0.25, 2.4}, {3.3, 0, 2.4}, {2.8, 0,
    2.475}, {2.8, -0.25, 2.475}, {3.525, -0.25, 2.49375},
    {3.525, 0, 2.49375}, {2.9, 0, 2.475}, {2.9, -0.15, 2.475},
    {3.45, -0.15, 2.5125}, {3.45, 0, 2.5125}, {2.8, 0, 2.4},
    {2.8, -0.15, 2.4}, {3.2, -0.15, 2.4}, {3.2, 0, 2.4}, {0, 0,
    3.15}, {0.8, 0, 3.15}, {0.8, -0.45, 3.15}, {0.45, -0.8,
    3.15}, {0, -0.8, 3.15}, {0, 0, 2.85}, {1.4, 0, 2.4}, {1.4,
    -0.784, 2.4}, {0.784, -1.4, 2.4}, {0, -1.4, 2.4}, {0.4, 0,
    2.55}, {0.4, -0.224, 2.55}, {0.224, -0.4, 2.55}, {0, -0.4,
    2.55}, {1.3, 0, 2.55}, {1.3, -0.728, 2.55}, {0.728, -1.3,
    2.55}, {0, -1.3, 2.55}, {1.3, 0, 2.4}, {1.3, -0.728, 2.4},
    {0.728, -1.3, 2.4}, {0, -1.3, 2.4}, {0, 0, 0}, {1.425,
    -0.798, 0}, {1.5, 0, 0.075}, {1.425, 0, 0}, {0.798, -1.425,
    0}, {0, -1.5, 0.075}, {0, -1.425, 0}, {1.5, -0.84, 0.075},
    {0.84, -1.5, 0.075}
};

/* *INDENT-ON* */

static void
fastTeapot(GLint grid, GLdouble scale)
{
  float p[4][4][3], q[4][4][3], r[4][4][3], s[4][4][3];
  long i, j, k, l;

  glEnable(GL_AUTO_NORMAL);
  glEnable(GL_MAP2_VERTEX_3);
  glRotatef(270.0, 1.0, 0.0, 0.0);
  glScalef(0.5 * scale, 0.5 * scale, 0.5 * scale);
  glTranslatef(0.0, 0.0, -1.5);
  for (i = 0; i < 10; i++) {
    for (j = 0; j < 4; j++) {
      for (k = 0; k < 4; k++) {
        for (l = 0; l < 3; l++) {
          p[j][k][l] = cpdata[patchdata[i][j * 4 + k]][l];
          q[j][k][l] = cpdata[patchdata[i][j * 4 + (3 - k)]][l];
          if (l == 1)
            q[j][k][l] *= -1.0;
          if (i < 6) {
            r[j][k][l] =
              cpdata[patchdata[i][j * 4 + (3 - k)]][l];
            if (l == 0)
              r[j][k][l] *= -1.0;
            s[j][k][l] = cpdata[patchdata[i][j * 4 + k]][l];
            if (l == 0)
              s[j][k][l] *= -1.0;
            if (l == 1)
              s[j][k][l] *= -1.0;
          }
        }
      }
    }
    glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4, 0, 1, 12, 4,
      &p[0][0][0]);
    glMapGrid2f(grid, 0.0, 1.0, grid, 0.0, 1.0);
    glEvalMesh2(GL_FILL, 0, grid, 0, grid);
    glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4, 0, 1, 12, 4,
      &q[0][0][0]);
    glEvalMesh2(GL_FILL, 0, grid, 0, grid);
    if (i < 6) {
      glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4, 0, 1, 12, 4,
        &r[0][0][0]);
      glEvalMesh2(GL_FILL, 0, grid, 0, grid);
      glMap2f(GL_MAP2_VERTEX_3, 0, 1, 3, 4, 0, 1, 12, 4,
        &s[0][0][0]);
      glEvalMesh2(GL_FILL, 0, grid, 0, grid);
    }
  }
}