#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
#	arm.c
#	chull.c
#	cube.c
#	dt.c
#	extreme.c
#	graham.c
#	i.poly.18
#	inpoly.c
#	int.c
#	macros.h
#	quadratic.c
#	sphere.c
#	stack.c
#	tri.c
#	vector.c
# This archive created: Wed May 22 13:11:23 1996
export PATH; PATH=/bin:$PATH
echo shar: extracting "'arm.c'" '(7257 characters)'
if test -f 'arm.c'
then
	echo shar: will not over-write existing file "'arm.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'arm.c'
	X/*
	X	arm.c
	X	Prints out one arm configuration to reach given target.
	X	Assumes number of links >= 3.
	X*/
	X
	X#include 	<stdio.h>
	X#include	<math.h>
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	X#define	NLINKS	20
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X
	X#define	FUZZ	1.0e-10
	X
	Xvoid	ReadTarget( tPointi target );
	Xint	ReadLinks( void );
	Xvoid	PrintLinks( void );
	Xvoid	PrintPointi( tPointi p );
	Xvoid	PrintPointd( tPointd p );
	Xbool	Solve3( int l1, int l2, int l3, tPointi target );
	Xint	TwoCircles( tPointi c1, int r1, tPointi c2, int r2, tPointd p );
	Xint	TwoCircles0a( int r1, tPointi c2, int r2, tPointd p );
	Xint	TwoCircles0b( int r1, tPointi c2, int r2, tPointd p );
	Xvoid	TwoCircles00( int r1, double a2, int r2, tPointd p );
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c );
	Xbool	EqPointi( tPointi a, tPointi b );
	Xint	Length2( tPointi v );
	Xbool	Solve2( int l1, int l2, tPointi target, tPointd J );
	Xbool	Nested( int l1, int l2, int l3, tPointi target );
	Xbool	Solven( int nlinks );
	Xdouble	RadDeg( double x );
	X
	Xint	linklen[NLINKS];	/* link lengths */
	Xint	nlinks;			/* number of links */
	XtPointi	target;			/* target point */
	X
	Xmain(argc, argv)
	Xint	argc;
	Xchar	*argv[];
	X{
	X
	X	ReadTarget( target );
	X	nlinks = ReadLinks( );
	X
	X	if ( !Solven( nlinks ) ) 
	X		printf("Solven: no solutions!\n");
	X}
	X
	X
	Xbool    Solven( int nlinks )
	X{
	X	int	i;
	X        int     m;              /* index of median link */
	X        int     l1, l2, l3;     /* length of links between kinks */
	X        int     totlength;      /* total length of all links */
	X        int     halflength;     /* floor of half of total */
	X
	X	/* Compute total & half length. */
	X        totlength = 0;
	X        for ( i = 0; i < nlinks; i++ ) 
	X                totlength += linklen[i];
	X        halflength = totlength / 2;
	X
	X	/* Find median link. */
	X        l1 = 0;
	X        for ( m = 0; m < nlinks; m++ ) {
	X                if ( (l1 + linklen[m]) > halflength)
	X                        break;
	X                l1 += linklen[m];
	X        }
	X
	X        l2 = linklen[m];
	X        l3 = totlength - l1 - l2;
	X        if ( Solve3( l1, l2, l3, target ) )
	X                return  TRUE;
	X        else    return  FALSE;
	X}
	X
	Xbool    Solve3( int l1, int l2, int l3, tPointi target )
	X{
	X	tPointd Jk;	/* coords of kinked joint returned by Solve2 */
	X	tPointi J1;	/* Joint1 on x-axis */
	X        tPointi Ttarget;/* translated target */
	X
	X	printf("==>Solve3: links = %d,%d,%d\n", l1, l2, l3);
	X	if ( Solve2( l1 + l2, l3, target, Jk ) ) {
	X		printf("Solve3: link1=%d, link2=%d, joint=",
	X			l1 + l2, l3);
	X		PrintPointd( Jk );
	X		return	TRUE;
	X	}
	X	else if ( Solve2( l1, l2 + l3, target, Jk ) ) {
	X		printf("Solve3: link1=%d, link2=%d, joint=",
	X			l1, l2 + l3);
	X		PrintPointd( Jk );
	X		return	TRUE;
	X	}
	X	else {	/* pin J0 to 0. */
	X		/* Shift so J1 is origin. */
	X		J1[X] = l1;	J1[Y] = 0;
	X		SubVec( target, J1, Ttarget );
	X		if ( Solve2( l2, l3, Ttarget, Jk ) ) {
	X			/* Shift solution back to origin. */
	X			Jk[X] += l1;
	X			printf("Solve3: link1=%d, link2=%d, link3=%d,\
	X				joints=", l1, l2, l3);
	X			PrintPointi( J1 );
	X			PrintPointd( Jk );
	X			return	TRUE;
	X		}
	X		else	return	FALSE;
	X	}
	X}
	Xbool    Solve2( int l1, int l2, tPointi target, tPointd J )
	X{
	X        tPointi c1 = {0,0};     /* center of circle 1 */
	X        int     nsoln;          /* # of solns: 0,1,2,3(infinite) */
	X
	X        nsoln = TwoCircles( c1, l1, target, l2, J );
	X	printf("<==Solve2: l1=%d, l2=%d; nsoln=%d\n", l1, l2, nsoln); 
	X        return	nsoln != 0;
	X}
	X
	X
	X/*
	X        TwoCircles finds the intersection points between two circles.
	X        This is the general routine, no assumptions.
	X	One intersection point is placed in p.
	X*/
	Xint     TwoCircles( tPointi c1, int r1, tPointi c2, int r2, tPointd p)
	X{
	X        tPointi c;
	X
	X	SubVec( c2, c1, c );
	X	return TwoCircles0a( r1, c, r2, p );
	X}
	X
	X/*
	X        TwoCircles0a assumes that the first circle is centered 
	X	on the origin.  It handles all the special cases, returning
	X        the number of intersection points: 0, 1, 2, 3 (inf).
	X*/
	Xint     TwoCircles0a( int r1, tPointi c2, int r2, tPointd p )
	X{
	X	int	dc2;	/* dist to center 2 squared */
	X	int	rplus2, rminus2;	/* (r1 +/- r2)^2 */
	X	double	f;	/* fraction along c2 for nsoln=1 */
	X
	X	/* Handle special cases. */
	X	dc2 = Length2( c2 );
	X	rplus2  = (r1 + r2) * (r1 + r2);
	X	rminus2 = (r1 - r2) * (r1 - r2);
	X
	X	/* No solution if c2 out of reach + or -. */
	X	if ( ( dc2 > rplus2 ) || ( dc2 < rminus2 ) )
	X		return	0;
	X
	X	/* One solution if c2 just reached. */
	X	/* Then solution is r1-of-the-way (f) to c2. */
	X	if ( dc2 == rplus2 ) {
	X		f = r1 / (double)(r1 + r2);
	X		p[X] = f * c2[X];	p[Y] = f * c2[Y];
	X		return	1;
	X	}
	X	if ( dc2 == rminus2 ) {
	X		if ( rminus2 == 0 ) {	/* Circles coincide. */
	X			p[X] = r1;	p[Y] = 0;
	X			return	3;
	X		}
	X		f = r1 / (double)(r1 - r2);
	X		p[X] = f * c2[X];	p[Y] = f * c2[Y];
	X		return	1;
	X	}
	X
	X        /* Two intersections. */
	X	return	TwoCircles0b( r1, c2, r2, p );
	X
	X}
	X/*
	X        TwoCircles0b also assumes that the first circle is centered       
	X        on the origin.  It rotates c2 to lie on the x-axis.
	X*/
	Xint     TwoCircles0b( int r1, tPointi c2, int r2, tPointd p )
	X{
	X        double  a2;     /* center of 2nd circle when rotated to x-axis */
	X        tPointd q;      /* one solution when c2 on x-axis */
	X        double  cost, sint;     /* sine and cosine of angle of c2 */
	X
	X	/* Rotate c2 to a2 on x-axis. */
	X        a2 = sqrt( (double) Length2( c2 ) );
	X        cost = c2[X] / a2;
	X        sint = c2[Y] / a2;
	X
	X        TwoCircles00( r1, a2, r2, q );
	X
	X        /* Rotate back */
	X        p[X] =  cost * q[X] + -sint * q[Y];
	X        p[Y] =  sint * q[X] +  cost * q[Y];
	X        
	X        return	2;
	X}
	X
	X/*
	X        TwoCircles00 assume circle one is origin-centered, and
	X        that the center of the second circle is at (a2,0).
	X*/
	Xvoid     TwoCircles00( int r1, double a2, int r2, tPointd p )
	X{
	X        /* Two intersections (only one returned in p). */
	X        p[X] = ( a2 + ( r1*r1 - r2*r2 ) / a2 ) / 2;
	X        p[Y] = sqrt( r1*r1 - p[X]*p[X] );
	X}
	Xvoid	ReadTarget( tPointi target )
	X{
	X	int	i;
	X
	X	for( i = 0; i < DIM; i++)
	X		scanf("%d", &target[i] );
	X	printf("Target (target) point = ");
	X	PrintPointi( target );
	X	putchar('\n');
	X}
	X
	Xvoid	PrintPointi( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X
	Xvoid	PrintPointd( tPointd p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%g", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X	putchar('\n');
	X}
	X
	Xint	ReadLinks( void )
	X{
	X	int	length;
	X
	X	nlinks = 0;
	X	while( scanf("%d", &length) != EOF ) {
	X		linklen[nlinks] = length;
	X		nlinks++;
	X	}
	X	PrintLinks( );
	X	return	nlinks;
	X}
	X		
	Xvoid	PrintLinks( )
	X{
	X	int	i;
	X
	X	for ( i=0; i < nlinks; i++)
	X		printf("%d:%d\t", i, linklen[i]);
	X	putchar('\n');
	X}
	Xbool	EqPointi( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i=0; i < DIM; i++ )
	X		if ( a[i] != b[i] )
	X			return FALSE;
	X	return TRUE;
	X}
	X
	Xint	Length2( tPointi v )
	X{
	X	int	i;
	X	int	ss;
	X
	X	ss = 0;
	X	for ( i=0; i < DIM; i++ )
	X		ss += v[i] * v[i];
	X	return ss;
	X}
	X
	X
	X
	Xdouble	RadDeg( double x )
	X{
	X	return 180 * x / M_PI;
	X}
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X	int	i;
	X
	X	for ( i=0; i < DIM; i++ ) {
	X		c[i] = a[i] - b[i];
	X	}
	X}
SHAR_EOF
if test 7257 -ne "`wc -c < 'arm.c'`"
then
	echo shar: error transmitting "'arm.c'" '(should have been 7257 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'chull.c'" '(28085 characters)'
if test -f 'chull.c'
then
	echo shar: will not over-write existing file "'chull.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'chull.c'
	X/*
	X	chull.c
	X	Written by Joseph O'Rourke
	X		and John Kutcher, Catherine Schevon, Susan Weller.
	X	Last modified: 14 August 1995
	X*/
	X#include <stdio.h>
	X#include <math.h>
	X
	X/* Define Boolean type */
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X/* Define vertex indices. */
	X#define X   0
	X#define Y   1
	X#define Z   2
	X
	X/* Define structures for vertices, edges and faces */
	Xtypedef struct tVertexStructure tsVertex;
	Xtypedef tsVertex *tVertex;
	X
	Xtypedef struct tEdgeStructure tsEdge;
	Xtypedef tsEdge *tEdge;
	X
	Xtypedef struct tFaceStructure tsFace;
	Xtypedef tsFace *tFace;
	X
	Xstruct tVertexStructure {
	X	int	v[3];
	X	int	vnum;
	X	tEdge	duplicate;	/* pointer to incident cone edge (or NULL) */
	X	bool 	onhull;		/* T iff point on hull. */
	X	bool	mark;		/* T iff point already processed. */
	X	tVertex next, prev;
	X};
	X
	Xstruct tEdgeStructure {
	X	tFace 	adjface[2];
	X	tVertex endpts[2];
	X	tFace	newface;	/* pointer to incident cone face. */
	X	bool 	delete;		/* T iff edge should be delete. */
	X	tEdge 	next, prev;
	X};
	X
	Xstruct tFaceStructure {
	X	tEdge 	edge[3];
	X	tVertex vertex[3];
	X	bool	visible;	/* T iff face visible from new point. */
	X	tFace 	next, prev;
	X};
	X
	X/* Define flags */
	X#define ONHULL   	TRUE
	X#define REMOVED  	TRUE
	X#define VISIBLE  	TRUE
	X#define PROCESSED	TRUE
	X
	X/* Global variable definitions */
	XtVertex vertices	= NULL;
	XtEdge edges    		= NULL;
	XtFace faces    		= NULL;
	Xbool debug = FALSE;
	X
	XtVertex vertices;
	XtEdge 	edges;
	XtFace 	faces;
	Xbool 	debug;
	X
	X/* Function declarations */
	XtVertex MakeVertex( void );
	Xvoid    ReadVertices( void );
	Xvoid    Print( void );
	Xvoid    Tetrahedron( void );
	Xvoid    ConstructHull( void );
	Xbool	AddOne( tVertex p );
	Xint     Volume6(tFace f, tVertex p);
	XtFace	MakeStructs( tEdge e, tVertex p );
	Xvoid    MakeCcw( tFace f, tEdge e, tVertex p );
	XtEdge   MakeEdge( void );
	XtFace   MakeFace( void );
	Xvoid    CleanUp( void );
	Xvoid    CleanEdges( void );
	Xvoid    CleanFaces( void );
	Xvoid    CleanVertices( void );
	Xbool	Collinear( tVertex a, tVertex b, tVertex c );
	Xvoid    CheckEuler(int V, int E, int F );
	Xdouble 	Volumed( tFace f, tVertex p );
	Xvoid	PrintPoint( tVertex p );
	Xvoid    Checks( void );
	Xvoid	Consistency( void );
	Xvoid	Convexity( void );
	Xvoid	PrintOut( tVertex v );
	Xvoid	PrintVertices( void );
	Xvoid	PrintEdges( void );
	Xvoid	PrintFaces( void );
	X
	X#include "macros.h"
	X
	X/*--------------------------------------------------------------------------*/
	Xmain( int argc, char *argv[] )
	X{
	X	char *s;
	X
	X	if ( argc > 1 && argv[1][0] == '-' )
	X             for ( s = argv[1] + 1; *s; ++s )
	X                   switch ( *s ) {
	X                          case 'd': debug = TRUE;
	X                                    fprintf( stderr, "Debug mode\n");
	X                                    break;
	X                          default : debug = FALSE;
	X                          }
	X
	X        else if ( argc > 1 && argv[1][0] != '-' ) {
	X             	printf ("Usage:  %s -d[ebug] \n", *argv );
	X		printf ("x y z coords of vertices from stdin\n");
	X             	exit(1);
	X	}
	X
	X        ReadVertices();
	X        Tetrahedron();
	X        ConstructHull();
	X        Print();
	X}
	X/*------------------------------------------------------------------
	XMakeVertex: Makes a vertex, nulls out fields.
	X--------------------------------------------------------------------*/ 
	XtVertex	MakeVertex( void )
	X{
	X	tVertex	v;
	X
	X	NEW( v, tsVertex );
	X	v->duplicate = NULL;
	X	v->onhull = !ONHULL;
	X	v->mark = !PROCESSED;
	X	ADD( vertices, v );
	X
	X	return v;
	X}
	X/*------------------------------------------------------------------
	XReadVertices: Reads in the vertices, and links them into a circular
	Xlist with MakeVertex.  There is no need for the # of vertices to be
	Xthe first line: the function looks for EOF instead.
	X--------------------------------------------------------------------*/ 
	Xvoid	ReadVertices( void )
	X{
	X	tVertex v;
	X	int 	x, y, z;
	X	int	vnum = 0;
	X
	X	while ( scanf ("%d %d %d", &x, &y, &z ) != EOF )  {
	X		v = MakeVertex();
	X		v->v[X] = x;
	X		v->v[Y] = y;
	X		v->v[Z] = z;
	X		v->vnum = vnum++;
	X		if ( abs(x) > 512 ) printf("coord %d might be too large\n", x);
	X		if ( abs(y) > 512 ) printf("coord %d might be too large\n", y);
	X		if ( abs(z) > 512 ) printf("coord %d might be too large\n", z);
	X	}
	X}
	X
	X/*------------------------------------------------------------------
	XPrint: Prints out the vertices and the faces.  Uses the vnum indices 
	Xcorresponding to the order in which the vertices were input.
	X--------------------------------------------------------------------*/
	Xvoid	Print( void )
	X{
	X	/* Pointers to vertices, edges, faces. */
	X	tVertex	v;
	X	tEdge   e;
	X	tFace   f;
	X	/* Counters for Euler's formula. */
	X	int 	V = 0, E = 0 , F = 0;
	X	/* Note: lowercase==pointer, uppercase==counter. */
	X
	X	/* Vertices. */
	X	printf("\n");
	X	v = vertices;
	X	do {                                 
	X		V++;           
	X		v = v->next;
	X	} while ( v != vertices );
	X	printf("Vertices:\tV = %d\n", V );
	X
	X	printf("index:\tx\ty\tz\n");
	X	do {                                 
	X		printf("%5d:\t%d\t%d\t%d\n", 
	X			v->vnum, v->v[X], v->v[Y], v->v[Z] );
	X		v = v->next;
	X	} while ( v != vertices );
	X	
	X	/* Faces. */
	X	f = faces;
	X	do {
	X		++F;                              
	X		f  = f ->next;
	X	} while ( f  != faces );
	X	printf("Faces:\t\tF = %d\n", F );
	X
	X	printf("\tv0\tv1\tv2\t(vertex indices)\n");
	X	do {                           
	X		printf("\t%d\t%d\t%d\n", 
	X			f->vertex[0]->vnum,
	X			f->vertex[1]->vnum, 
	X			f->vertex[2]->vnum );
	X		f = f->next;
	X	} while ( f != faces );
	X
	X	/* Edges. */	
	X	e = edges;
	X        do {
	X                E++;
	X                e = e->next;
	X        } while ( e != edges );
	X        printf("Edges:\tE = %d\n", E );
	X	/* Edges not printed out (but easily added). */
	X	debug = TRUE;
	X	CheckEuler( V, E, F );
	X
	X}
	X
	X/*-----------------------------------------------------------------------
	X Tetrahedron builds the initial tetrahedron.  It first finds 3 noncollinear
	X points and makes a face out of them, and then finds a fourth point that
	X is not coplanar with that face.  The vertices are stored in the face
	X structure in counterclockwise order so that the volume between the face
	X and the point is negative. Lastly, the 3 newfaces to the fourth point
	X are constructed and the data structures are cleaned up. 
	X -----------------------------------------------------------------------*/
	Xvoid	Tetrahedron( void )
	X{
	X	tVertex v1, v4, t;
	X	tFace 	f;
	X	tEdge 	e1, e2, e3, s;
	X	int 	vol;
	X	
	X
	X	/* Find 3 non-Collinear points. */
	X	v1 = vertices;
	X	while ( Collinear( v1, v1->next, v1->next->next ) ) 
	X      	      	if ( ( v1 = v1->next ) == vertices ) {
	X         	     printf("All points are Collinear!\n");
	X         	     exit(0);
	X        	} 
	X	
	X	/* Mark the vertices as processed. */
	X	v1->mark = PROCESSED;  
	X	v1->next->mark = PROCESSED; 
	X	v1->next->next->mark = PROCESSED;              
	X                                               
	X	/* Create edges of the initial triangle. */
	X	e1 = MakeEdge(); 
	X	e2 = MakeEdge();
	X	e3 = MakeEdge();
	X	e1->endpts[0] = v1;              e1->endpts[1] = v1->next;
	X	e2->endpts[0] = v1->next;        e2->endpts[1] = v1->next->next;
	X	e3->endpts[0] = v1->next->next;  e3->endpts[1] = v1;
	X
	X	/* Create face for triangle. */
	X	f = MakeFace();   
	X	f->edge[0] = e1;   f->edge[1] = e2;   f->edge[2] = e3;
	X	f->vertex[0] = v1;  f->vertex[1] = v1->next;
	X	f->vertex[2] = v1->next->next;
	X
	X	/* Link edges to face. */
	X	e1->adjface[0] = e2->adjface[0] = e3->adjface[0] = f; 
	X	
	X	
	X	/* Find a fourth, non-coplanar point to form tetrahedron. */
	X	v4 = v1->next->next->next;
	X	vol = Volume6( f, v4 );
	X	while ( !vol )   {
	X      		if ( ( v4 = v4->next ) == v1 ) {
	X         	       printf("All points are coplanar!\n");
	X         	       exit(0);
	X        	} 
	X       		vol = Volume6( f, v4 );
	X       	}
	X	v4->mark = PROCESSED;
	X
	X	/* Store vertices in ccw order. */
	X	if( vol < 0 ) {
	X            SWAP( t, f->vertex[1], f->vertex[2] );
	X            SWAP( s, f->edge[1], f->edge[2] );
	X  	}
	X
	X	
	X	/* Construct the faces and edges between the original
	X   	   triangle and the fourth point. */
	X	e1->adjface[1] = MakeStructs( e1, v4 );
	X	e2->adjface[1] = MakeStructs( e2, v4 );
	X	e3->adjface[1] = MakeStructs( e3, v4 );
	X	
	X	CleanUp();
	X}
	X
	X/*-------------------------------------------------------------------------
	XConstructHull adds the vertices to the hull one at a time.  The hull
	Xvertices are those in the list marked as onhull.
	X  -------------------------------------------------------------------------*/
	Xvoid	ConstructHull( void )
	X{
	X	tVertex	v;
	X	tFace 	f;
	X	int 	vol;
	X	bool	changed;	/* T if addition changes hull; not used. */
	X
	X	v = vertices;
	X	f = faces;
	X	do {
	X		if ( !v->mark ) {
	X			v->mark = PROCESSED;
	X			changed = AddOne( v );
	X			CleanUp();
	X			if ( debug ) {
	X				fprintf(stderr,"after cleanup:\n");
	X				PrintOut( v );
	X				Checks();
	X			}
	X		}
	X		v = v->next;
	X	} while ( v != vertices );
	X}
	X/*-------------------------------------------------------------------------
	XAddOne is passed a vertex.  It first determines all faces visible from 
	Xthat point.  If none are visible then the point is marked as not onhull.
	XNext is a loop over edges.  If both faces adjacent to an edge are
	Xvisible, then the edge is marked for deletion.  If just one of the adjacent 
	Xfaces is visible then a new face is constructed.
	X--------------------------------------------------------------------------*/
	Xbool 	AddOne( tVertex p )
	X{
	X	tFace 	f; 
	X	tEdge 	e;
	X        int 	vol;
	X	bool	vis;
	X
	X	/* Mark faces visible from p. */
	X        f = faces;
	X	do {
	X		vol = Volume6( f, p );
	X/*
	X		if (debug) fprintf(stderr,"faddr: %6x   paddr: %6x   Vol = %d\n",f,p,vol);
	X*/
	X		if ( vol < 0 ) {
	X			f->visible = VISIBLE;  
	X			vis = TRUE;                      
	X		}
	X		f = f->next;
	X	} while ( f != faces );
	X
	X	/* If no faces are visible from p, then p is inside the hull. */
	X        if ( !vis ) {
	X		p->onhull = !ONHULL;  
	X		return FALSE; 
	X	}
	X
	X	/* Mark edges in interior of visible region for deletion.
	X	   Erect a newface based on each border edge. */
	X        e = edges;
	X	do {
	X		tEdge temp;
	X		temp = e->next;
	X		if ( e->adjface[0]->visible && e->adjface[1]->visible )
	X			/* e interior: mark for deletion. */
	X			e->delete = REMOVED;
	X		else if ( e->adjface[0]->visible || e->adjface[1]->visible ) 
	X			/* e border: make a new face. */
	X			e->newface = MakeStructs( e, p );
	X		e = temp;
	X	} while ( e != edges );
	X	return TRUE;
	X}
	X
	X/*-------------------------------------------------------------------------
	XVolume6 returns six times the volume of the tetrahedron determined by f 
	Xand p.  Volume6 is positive iff p is on the negative side of f,
	Xwhere the positive side is determined by the rh-rule.  So the volume 
	Xis positive if the ccw normal to f points outside the tetrahedron.
	X--------------------------------------------------------------------------*/
	Xint 	Volume6( tFace f, tVertex p )
	X{
	X	int 	vol;
	X	int 	ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz;
	X        int	bxdx, bydy, bzdz, cxdx, cydy, czdz;
	X	double	vold;
	X	int	i;
	X
	X	ax = f->vertex[0]->v[X];
	X	ay = f->vertex[0]->v[Y];
	X	az = f->vertex[0]->v[Z];
	X	bx = f->vertex[1]->v[X];
	X	by = f->vertex[1]->v[Y];
	X	bz = f->vertex[1]->v[Z];
	X	cx = f->vertex[2]->v[X];
	X	cy = f->vertex[2]->v[Y];
	X	cz = f->vertex[2]->v[Z];
	X	dx = p->v[X];
	X	dy = p->v[Y];
	X	dz = p->v[Z];
	X
	X	/* This is the expression used in the text.  Now replaced.
	X	vol = 	 -az * by * cx + ay * bz * cx + az * bx * cy - ax * bz * cy 
	X		- ay * bx * cz + ax * by * cz + az * by * dx - ay * bz * dx 
	X		- az * cy * dx + bz * cy * dx + ay * cz * dx - by * cz * dx 
	X		- az * bx * dy + ax * bz * dy + az * cx * dy - bz * cx * dy 
	X		- ax * cz * dy + bx * cz * dy + ay * bx * dz - ax * by * dz 
	X		- ay * cx * dz + by * cx * dz + ax * cy * dz - bx * cy * dz;
	X	*/
	X	/* This expression is algebraically equivalent to the above, but
	X	uses fewer multiplications.
	X	vol =	 -(az-dz) * (by-dy) * (cx-dx) 
	X		+ (ay-dy) * (bz-dz) * (cx-dx) 
	X		+ (az-dz) * (bx-dx) * (cy-dy) 
	X		- (ax-dx) * (bz-dz) * (cy-dy) 
	X		- (ay-dy) * (bx-dx) * (cz-dz) 
	X		+ (ax-dx) * (by-dy) * (cz-dz); */
	X	/* And this one has even fewer arithmetic operations:
	X	   (thanks to Robert Fraczkiewicz): */
	X        bxdx=bx-dx;
	X        bydy=by-dy;
	X        bzdz=bz-dz;
	X        cxdx=cx-dx;
	X        cydy=cy-dy;
	X        czdz=cz-dz;
	X        vol =    (az-dz)*(bxdx*cydy-bydy*cxdx)
	X               + (ay-dy)*(bzdz*cxdx-bxdx*czdz)
	X               + (ax-dx)*(bydy*czdz-bzdz*cydy);
	X
	X	/* Compare integer volume with double volume for saftey. */
	X	vold = Volumed( f, p );
	X	if (debug)
	X	{
	X		fprintf(stderr,"Face = %6x\tVertex = %d, vol(int) = %d\tvol(double) = %lf\n",
	X			f,p->vnum,vol,vold);
	X	}
	X	if ( fabs( vol - vold ) >= 1.0 ) {
	X		printf("Likely integer overflow in volume calculation\n");
	X		printf("because coordinates are too large:\n");
	X		printf("\tvol(int) = %d;\tvol(double) = %lf\n");
	X		printf("\tfour points:\n");
	X		for ( i=0; i < 3; i++ )
	X			PrintPoint( f->vertex[i] );
	X		PrintPoint( p );
	X		printf("Basing decision based on double volume.\n");
	X		/* Return based on double volume. */
	X		if ( vold >= 1.0 )
	X			return	1;
	X		else if ( vold <= -1.0 )
	X			return	-1;
	X		else	return	0;
	X	}
	X	return vol;
	X}		
	X/*-------------------------------------------------------------------------
	XVolumed is the same as Volume6 but computed with doubles.  For protection
	Xagainst overflow.
	X--------------------------------------------------------------------------*/
	Xdouble 	Volumed( tFace f, tVertex p )
	X{
	X	double	vol;
	X	double 	ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz;
	X        double	bxdx, bydy, bzdz, cxdx, cydy, czdz;
	X
	X	ax = f->vertex[0]->v[X];
	X	ay = f->vertex[0]->v[Y];
	X	az = f->vertex[0]->v[Z];
	X	bx = f->vertex[1]->v[X];
	X	by = f->vertex[1]->v[Y];
	X	bz = f->vertex[1]->v[Z];
	X	cx = f->vertex[2]->v[X];
	X	cy = f->vertex[2]->v[Y];
	X	cz = f->vertex[2]->v[Z];
	X	dx = p->v[X];
	X	dy = p->v[Y];
	X	dz = p->v[Z];
	X
	X        bxdx=bx-dx;
	X        bydy=by-dy;
	X        bzdz=bz-dz;
	X        cxdx=cx-dx;
	X        cydy=cy-dy;
	X        czdz=cz-dz;
	X        vol =    (az-dz)*(bxdx*cydy-bydy*cxdx)
	X               + (ay-dy)*(bzdz*cxdx-bxdx*czdz)
	X               + (ax-dx)*(bydy*czdz-bzdz*cydy);
	X
	X
	X	return vol;
	X}
	X
	X/*-------------------------------------------------------------------------*/
	Xvoid	PrintPoint( tVertex p )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < 3; i++ )
	X		printf("\t%d", p->v[i]);
	X	putchar('\n');
	X}
	X/*----------------------------------------------------------------------
	XMakeStructs makes a new face and two new edges between the 
	Xedge and the point that are passed to it. It returns a pointer to
	Xthe new face.
	X----------------------------------------------------------------------*/
	XtFace	MakeStructs( tEdge e, tVertex p )
	X{
	X	tEdge 	new_edge[2];
	X        tFace 	new_face;
	X        int 	i, j;
	X
	X	/* Make two new edges (if don't already exist). */
	X        for ( i=0; i < 2; ++i ) 
	X		/* If the edge exists, copy it into new_edge. */
	X		if ( !( new_edge[i] = e->endpts[i]->duplicate) ) {
	X			/* Otherwise (duplicate is NULL), MakeEdge. */
	X			new_edge[i] = MakeEdge();
	X			new_edge[i]->endpts[0] = e->endpts[i];
	X			new_edge[i]->endpts[1] = p;
	X			e->endpts[i]->duplicate = new_edge[i];
	X		}
	X
	X        /* Make the new face. */
	X        new_face = MakeFace();   
	X        new_face->edge[0] = e;
	X        new_face->edge[1] = new_edge[0];
	X        new_face->edge[2] = new_edge[1];
	X        MakeCcw( new_face, e, p ); 
	X        
	X        /* Set the adjacent face pointers. */
	X        for ( i=0; i < 2; ++i )
	X	for ( j=0; j < 2; ++j )  
	X		/* Only one NULL link should be set to new_face. */
	X		if ( !new_edge[i]->adjface[j] ) {
	X			new_edge[i]->adjface[j] = new_face;
	X			break;
	X		}
	X        
	X        return new_face;
	X}
	X
	X/*------------------------------------------------------------------------
	XMakeCcw puts the vertices in the face structure in counterclockwise order.  
	XIf there is no adjacent face[1] then we know that
	Xwe are working with the first face of the initial tetrahedron.  In this
	Xcase we want to store the vertices in the opposite order from the 
	Xinitial face.  Otherwise, we want to store the vertices in the same order
	Xas in the visible face.  The third vertex is always p.
	X------------------------------------------------------------------------*/
	Xvoid	MakeCcw( tFace f, tEdge e, tVertex p )
	X{
	X        int 	i;	/* Index */
	X        tFace 	fi;	/* The invisible face adjacent to e */
	X	tEdge	s;	/* Temporary, for swapping */
	X      
	X	/* If this is the initial tetrahedron, then e has only one
	X	   adjacent face, and use that for fi.  Otherwise, use the
	X	   invisible face. */ 
	X	if ( !e->adjface[1] )
	X		fi = e->adjface[0];
	X	else {
	X		if  ( !e->adjface[0]->visible )      
	X			fi = e->adjface[0];
	X		else	fi = e->adjface[1];
	X	}
	X       
	X	/* Set vertex[0] & [1] of f to have the opposite orientation
	X	   as do the corresponding vertices of fi. */ 
	X	/* Find the index i of e->endpoint[1] in fi. */
	X	for ( i=0; fi->vertex[i] != e->endpts[1]; ++i )
	X		;
	X	/* Orient f opposite that of fi. */
	X	if ( fi->vertex[ (i+1) % 3 ] != e->endpts[0] ) {
	X		f->vertex[0] = e->endpts[1];  
	X		f->vertex[1] = e->endpts[0];    
	X	}
	X	else {                               
	X		f->vertex[0] = e->endpts[0];   
	X		f->vertex[1] = e->endpts[1];      
	X		SWAP( s, f->edge[1], f->edge[2] );
	X	}
	X 
	X        f->vertex[2] = p;
	X}
	X 
	X/*---------------------------------------------------------------------
	XMakeEdge creates a new cell and initializes all pointers to NULL
	Xand sets all flags to off.  It returns a pointer to the empty cell.
	X---------------------------------------------------------------------*/
	XtEdge 	MakeEdge( void )
	X{
	X	tEdge 	e;
	X
	X        NEW( e, tsEdge );
	X        e->adjface[0] = e->adjface[1] = e->newface = NULL;
	X        e->endpts[0] = e->endpts[1] = NULL;
	X        e->delete = !REMOVED;
	X        ADD( edges, e );
	X        return e;
	X}
	X
	X/*---------------------------------------------------------------------
	XMakeFace creates a new face structure and initializes all of its
	Xflags to NULL and sets all the flags to off.  It returns a pointer
	Xto the empty cell.
	X----------------------------------------------------------------------*/
	XtFace 	MakeFace( void )
	X{
	X	tFace 	f;
	X	int 	i;
	X
	X	NEW( f, tsFace);
	X	for ( i=0; i < 3; ++i ) {
	X    	      f->edge[i] = NULL;
	X    	      f->vertex[i] = NULL;
	X	}
	X	f->visible = !VISIBLE;
	X	ADD( faces, f );
	X	return f;
	X}
	X
	X/*-------------------------------------------------------------------------
	XCleanUp goes through each data structure list and clears all
	Xflags and NULLs out some pointers.  The order of processing
	X(edges, faces, vertices) is important.
	X------------------------------------------------------------------------*/
	Xvoid	CleanUp( void )
	X{
	X	CleanEdges();
	X	CleanFaces();
	X	CleanVertices();
	X}
	X
	X/*------------------------------------------------------------------------
	XCleanEdges runs through the edge list and cleans up the structure.
	XIf there is a newface then it will put that face in place of the 
	Xvisible face and NULL out newface. It also deletes so marked edges.
	X-----------------------------------------------------------------------*/
	Xvoid	CleanEdges( void )
	X{
	X	tEdge 	e;	/* Primary index into edge list. */
	X	tEdge 	t;	/* Temporary edge pointer. */
	X		
	X	/* Integrate the newface's into the data structure. */
	X	/* Check every edge. */
	X	e = edges;
	X	do {
	X		if ( e->newface ) { 
	X			if ( e->adjface[0]->visible )
	X				e->adjface[0] = e->newface; 
	X			else	e->adjface[1] = e->newface;
	X			e->newface = NULL;
	X		}
	X		e = e->next;
	X	} while ( e != edges );
	X
	X	/* Delete any edges marked for deletion. */
	X	while ( edges && edges->delete ) { 
	X		e = edges;
	X		DELETE( edges, e );
	X	}
	X	e = edges->next;
	X	do {
	X		if ( e->delete ) {
	X			t = e;
	X			e = e->next;
	X			DELETE( edges, t );
	X		}
	X		else e = e->next;
	X	} while ( e != edges );
	X}
	X
	X/*------------------------------------------------------------------------
	XCleanFaces runs through the face list and deletes any face marked visible.
	X-----------------------------------------------------------------------*/
	Xvoid	CleanFaces( void )
	X{
	X	tFace 	f;	/* Primary pointer into face list. */
	X	tFace 	t;	/* Temporary pointer, for deleting. */
	X	
	X
	X	while ( faces && faces->visible ) { 
	X		f = faces;
	X		DELETE( faces, f );
	X	}
	X	f = faces->next;
	X	do {
	X		if ( f->visible ) {
	X			t = f;
	X			f = f->next;
	X			DELETE( faces, t );
	X		}
	X		else f = f->next;
	X	} while ( f != faces );
	X}
	X
	X/*-------------------------------------------------------------------------
	XCleanVertices runs through the vertex list and deletes the 
	Xvertices that are marked as processed but are not incident to any undeleted
	Xedges. 
	X-------------------------------------------------------------------------*/
	Xvoid	CleanVertices( void )
	X{
	X	tEdge 	e;
	X	tVertex	v, t;
	X	
	X	/* Mark all vertices incident to some undeleted edge as on the hull. */
	X	e = edges;
	X	do {
	X	         e->endpts[0]->onhull = e->endpts[1]->onhull = ONHULL;
	X	         e = e->next;
	X	} while (e != edges);
	X	
	X	/* Delete all vertices that have been processed but
	X	   are not on the hull. */
	X	while ( vertices && vertices->mark && !vertices->onhull ) { 
	X		v = vertices;
	X		DELETE( vertices, v );
	X	}
	X	v = vertices->next;
	X	do {
	X		if ( v->mark && !v->onhull ) {    
	X			t = v; 
	X			v = v->next;
	X			DELETE( vertices, t )
	X		}
	X		else v = v->next;
	X	} while ( v != vertices );
	X	
	X	/* Reset flags. */
	X	v = vertices;
	X	do {
	X		v->duplicate = NULL; 
	X		v->onhull = !ONHULL; 
	X		v = v->next;
	X	} while ( v != vertices );
	X}
	X/*------------------------------------------------------------------
	XCollinear checks to see if the three points given are collinear,
	Xby checking to see if each element of the cross product is zero.
	X---------------------------------------------------------------------*/
	Xbool	Collinear( tVertex a, tVertex b, tVertex c )
	X{
	X	return ( c->v[Z] - a->v[Z] ) * ( b->v[Y] - a->v[Y] ) -
	X               ( b->v[Z] - a->v[Z] ) * ( c->v[Y] - a->v[Y] ) == 0
	X            && ( b->v[Z] - a->v[Z] ) * ( c->v[X] - a->v[X] ) -
	X               ( b->v[X] - a->v[X] ) * ( c->v[Z] - a->v[Z] ) == 0
	X            && ( b->v[X] - a->v[X] ) * ( c->v[Y] - a->v[Y] ) -
	X               ( b->v[Y] - a->v[Y] ) * ( c->v[X] - a->v[X] ) == 0  ;
	X}
	X
	X
	X/*------------------------------------------------------------------------
	XConsistency runs through the edge list and checks that all
	Xadjacent faces have their endpoints in opposite order.  This verifies
	Xthat the vertices are in counterclockwise order.
	X-----------------------------------------------------------------------*/
	Xvoid	Consistency( void )
	X{
	X	register tEdge e;
	X	register int i, j;
	X
	X	e = edges;
	X
	X	do {
	X    	   /* find index of endpoint[0] in adjacent face[0] */
	X           for ( i = 0; e->adjface[0]->vertex[i] != e->endpts[0]; ++i )
	X                 ;
	X   
	X   	   /* find index of endpoint[0] in adjacent face[1] */
	X   	   for ( j = 0; e->adjface[1]->vertex[j] != e->endpts[0]; ++j )
	X         	 ;
	X
	X           /* check if the endpoints occur in opposite order */
	X           if ( !( e->adjface[0]->vertex[ (i+1) % 3 ] ==
	X                   e->adjface[1]->vertex[ (j+2) % 3 ] ||
	X                   e->adjface[0]->vertex[ (i+2) % 3 ] ==
	X                   e->adjface[1]->vertex[ (j+1) % 3 ] )  )
	X               break;
	X           e = e->next;
	X
	X           } while ( e != edges );
	X
	X        if ( e != edges )
	X             fprintf( stderr, "Checks: edges are NOT consistent.\n");
	X        else
	X             fprintf( stderr, "Checks: edges consistent.\n");
	X
	X}
	X
	X/*----------------------------------------------------------------------
	XConvexity checks that the volume between every face and every
	Xpoint is negative.  This shows that each point is inside every face
	Xand therefore the hull is convex.
	X---------------------------------------------------------------------*/
	Xvoid	Convexity( void )
	X{
	X	register tFace f;
	X        register tVertex v;
	X        int vol;
	X
	X        f = faces;
	X
	X        do {
	X              v = vertices;
	X              do {
	X		   if ( v->mark ) {
	X                      vol = Volume6( f, v );
	X                      if ( vol < 0 )
	X			break;
	X		   }
	X                   v = v->next;
	X                   } while ( v != vertices );
	X
	X	      f = f->next;
	X
	X            } while ( f != faces );
	X
	X        if ( f != faces )
	X           fprintf( stderr, "Checks: NOT convex.\n");
	X        else if ( debug ) 
	X           fprintf( stderr, "Checks: convex.\n");
	X}
	X
	X/*----------------------------------------------------------------------
	XCheckEuler checks Euler's relation, as well as its implications when
	Xall faces are known to be triangles.  Only prints positive information
	Xwhen debug is true, but always prints negative information.
	X  ---------------------------------------------------------------------*/
	Xvoid	CheckEuler( int V, int E, int F )
	X{
	X	if ( debug )
	X             fprintf( stderr, "Checks: V, E, F = %d %d %d:\t", V, E, F);
	X
	X	if ( (V - E + F) != 2 )
	X             fprintf( stderr, "Checks: V-E+F != 2\n");
	X        else if ( debug )
	X             fprintf( stderr, "V-E+F = 2\t");
	X
	X
	X        if ( F != (2 * V - 4) )
	X             fprintf( stderr, "Checks: F=%d != 2V-4=%d; V=%d\n",
	X                F, 2*V-4, V);
	X        else if ( debug ) 
	X             fprintf( stderr, "F = 2V-4\t");
	X
	X        if ( (2 * E) != (3 * F) )
	X             fprintf( stderr, "Checks: 2E=%d != 3F=%d; E=%d, F=%d\n",
	X                2*E, 3*F, E, F );
	X        else if ( debug ) 
	X             fprintf( stderr, "2E = 3F\n");
	X}
	X/*-----------------------------------------------------------------------*/
	Xvoid	Checks( void )
	X{
	X        tVertex v;
	X        tEdge   e;
	X        tFace   f;
	X        int 	V = 0, E = 0 , F = 0;
	X
	X        Consistency();
	X        Convexity();
	X        if ( v = vertices )
	X             do {
	X                if (v->mark) V++;
	X                v = v->next;
	X                } while ( v != vertices );
	X        if ( e = edges )
	X             do {
	X                 E++;
	X                e = e->next;
	X             } while ( e != edges );
	X        if ( f = faces )
	X             do {
	X                F++;
	X                f  = f ->next;
	X                } while ( f  != faces );
	X	CheckEuler( V, E, F );
	X}
	X/*================================================================
	XThese functions are used whenever the debug flag is set.
	XThey print out the entire contents of each data structure.  
	XPrinting is to standard error.  To grab the output in a file in the csh, 
	Xuse this:
	X	chull < i.file >&! o.file
	X==================================================================*/
	X/*-------------------------------------------------------------------*/
	Xvoid	PrintOut( tVertex v )
	X{
	X	fprintf( stderr, "\nAdding vertex %6x :\n", v );
	X	PrintVertices();
	X	PrintEdges();
	X	PrintFaces();
	X}
	X
	X/*-------------------------------------------------------------------------*/
	Xvoid	PrintVertices( void )
	X{
	X	tVertex temp;
	X
	X	temp = vertices;
	X	fprintf (stderr, "Vertex List\n");
	X	if (vertices) do {
	X            fprintf(stderr,"  addr %6x\t", vertices );
	X            fprintf(stderr,"  vnum %4d", vertices->vnum );
	X            fprintf(stderr,"   (%6d,%6d,%6d)",vertices->v[X],
	X                    vertices->v[Y], vertices->v[Z] );
	X            fprintf(stderr,"   active:%3d", vertices->onhull );
	X            fprintf(stderr,"   dup:%5x", vertices->duplicate );
	X            fprintf(stderr,"   mark:%2d\n", vertices->mark );
	X            vertices = vertices->next;
	X            } while ( vertices != temp );
	X
	X}
	X/*-------------------------------------------------------------------------*/
	Xvoid	PrintEdges( void )
	X{
	X	tEdge 	temp;
	X	int 	i;
	X	
	X	temp = edges;
	X	fprintf (stderr, "Edge List\n");
	X	if (edges) do {
	X            fprintf( stderr, "  addr: %6x\t", edges );
	X            fprintf( stderr, "adj: ");
	X            for (i=0; i<2; ++i) 
	X                 fprintf( stderr, "%6x", edges->adjface[i] );
	X            fprintf( stderr, "  endpts:");
	X            for (i=0; i<2; ++i) 
	X                 fprintf( stderr, "%4d", edges->endpts[i]->vnum);
	X            fprintf( stderr, "  del:%3d\n", edges->delete );
	X            edges = edges->next; 
	X        } while (edges != temp );
	X
	X}
	X/*-------------------------------------------------------------------------*/
	Xvoid	PrintFaces( void )
	X{
	X	int 	i;
	X	tFace 	temp;
	X
	X	temp = faces;
	X	fprintf (stderr, "Face List\n");
	X	if (faces) do {
	X            fprintf(stderr, "  addr: %6x\t", faces );
	X            fprintf(stderr, "  edges:");
	X            for( i=0; i<3; ++i )
	X                 fprintf(stderr, "%6x", faces->edge[i] );
	X            fprintf(stderr, "  vert:");
	X            for ( i=0; i<3; ++i)
	X                  fprintf(stderr, "%4d", faces->vertex[i]->vnum );
	X            fprintf(stderr, "  vis: %d\n", faces->visible );
	X            faces= faces->next;
	X        } while ( faces != temp );
	X
	X}
SHAR_EOF
if test 28085 -ne "`wc -c < 'chull.c'`"
then
	echo shar: error transmitting "'chull.c'" '(should have been 28085 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'cube.c'" '(1065 characters)'
if test -f 'cube.c'
then
	echo shar: will not over-write existing file "'cube.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'cube.c'
	X/*-----------------------------------------------------------------------
	X      Cube                                     Susan L. Weller
	X  -----------------------------------------------------------------------
	X           This program will generate a given number of points that
	X   lie inside a cube.  The number of points is given on command line.
	X  -----------------------------------------------------------------------*/
	X
	X#include <stdio.h>
	X#include <math.h>
	X
	X#define MAX_INT	   2147483647
	X
	Xmain( argc, argv )
	Xint argc;
	Xchar *argv[];
	X{
	Xint n;
	Xdouble x, y, z;
	X
	Xsrandom( (int) time( (long *) 0 ) );
	Xif ( argc != 2 ) {
	X     printf("Usage: cube <number of random points>\n");
	X     exit(1);
	X     }
	X
	Xn = atoi( argv[1] );
	X
	Xwhile ( n-- ) {
	X              x = 2.0 * (double) random() / MAX_INT - 1.0;
	X              y = 2.0 * (double) random() / MAX_INT - 1.0;
	X              z = 2.0 * (double) random() / MAX_INT - 1.0;
	X              printf ("%6d  %6d  %6d\n", (int) (100*x ),  (int) ( 100*y ),
	X                        (int) (100 * z) );
	X              }
	X
	X}
	X
	X
	X            
SHAR_EOF
if test 1065 -ne "`wc -c < 'cube.c'`"
then
	echo shar: error transmitting "'cube.c'" '(should have been 1065 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'dt.c'" '(1147 characters)'
if test -f 'dt.c'
then
	echo shar: will not over-write existing file "'dt.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'dt.c'
	X#include	<stdio.h>
	X#define	NMAX	100
	Xmain()
	X{
	X	int	x[NMAX],y[NMAX],z[NMAX];/* input points xy,z=x^2+y^2 */
	X	int	n;			/* number of input points */
	X	int	i, j, k, m;		/* indices of four points */
	X	int	xn, yn, zn;		/* outward vector normal to (i,j,k) */
	X	int	flag;			/* t if m above of (i,j,k) */
	X
	X	/* Input points and compute z = x^2 + y^2. */
	X	scanf("%d", &n);
	X	for ( i = 0; i < n; i++ ) {
	X		scanf("%d %d", &x[i], &y[i]);
	X		z[i] = x[i] * x[i] + y[i] * y[i];
	X	}
	X
	X	/* For each triple (i,j,k) */
	X	for ( i = 0; i < n - 2; i++ )
	X	for ( j = i + 1; j < n; j++ )
	X	for ( k = i + 1; k < n; k++ ) 
	X	if ( j != k ) {
	X		/* Compute normal to triangle (i,j,k). */
	X		xn = (y[j]-y[i])*(z[k]-z[i]) - (y[k]-y[i])*(z[j]-z[i]);
	X		yn = (x[k]-x[i])*(z[j]-z[i]) - (x[j]-x[i])*(z[k]-z[i]);
	X		zn = (x[j]-x[i])*(y[k]-y[i]) - (x[k]-x[i])*(y[j]-y[i]);
	X
	X		/* Only examine faces on bottom of paraboloid: zn < 0. */
	X		if ( flag = (zn < 0) )
	X			/* For each other point m */
	X			for (m = 0; m < n; m++)
	X				/* Check if m above (i,j,k). */
	X				flag = flag && 
	X					((x[m]-x[i])*xn +
	X					 (y[m]-y[i])*yn +
	X					 (z[m]-z[i])*zn <= 0);
	X		if (flag)
	X			printf("%d\t%d\t%d\n", i, j, k);
	X	}
	X}
SHAR_EOF
if test 1147 -ne "`wc -c < 'dt.c'`"
then
	echo shar: error transmitting "'dt.c'" '(should have been 1147 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'extreme.c'" '(5406 characters)'
if test -f 'extreme.c'
then
	echo shar: will not over-write existing file "'extreme.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'extreme.c'
	X/*
	X	extreme.c
	X	Written by Joseph O'Rourke
	X	orourke@cs.smith.edu
	X*/
	X#include	<stdio.h>
	X#include	<math.h>
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X#define PMAX    1000            /* Max # of pts in polygon */
	X   
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c );
	Xint     Dot( tPointi a, tPointi b );
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] );
	Xvoid	PrintPoint( tPointi p );
	Xint	Midway( int a, int b, int n );
	Xint	Extreme( tPointi u, tPolygoni P, int n );
	X
	Xmain()
	X{
	X	int	n, c;
	X	tPolygoni	P;
	X	tPointi	u;
	X
	X	n = ReadPoly( P );
	X	u[0] = 0; u[1] = 1;
	X	c = Extreme( u, P, n );
	X	printf("Extreme returns %d\n", c);
	X	u[0] = 0; u[1] = -1;
	X	c = Extreme( u, P, n );
	X	printf("Extreme returns %d\n", c);
	X	u[0] = 1; u[1] = 0;
	X	c = Extreme( u, P, n );
	X	printf("Extreme returns %d\n", c);
	X	u[0] = -1; u[1] = 0;
	X	c = Extreme( u, P, n );
	X	printf("Extreme returns %d\n", c);
	X}
	X
	X/*
	X	Returns index midway ccw from a to b, mod n.
	X*/
	Xint	Midway( int a, int b, int n )
	X{
	X	if (a < b)
	X		return	( a + b ) / 2;
	X	else
	X		return	( ( a + b + n ) / 2 ) % n;
	X}
	X
	X/*
	X	Returns index of a point extreme in direction u.
	X*/
	Xint	Extreme( tPointi u, tPolygoni P, int n )
	X{
	X	int	a,a1, b;/* [a,b] includes extreme; a1=a+1. */
	X	int	c, c1;	/* index midway; c1 is c +- 1. */
	X	tPointi	A, C, B;/* edge vectors after a, after c, before c. */
	X	int	Adot, Cdot, Bdot;	/* dots with u. */
	X	int	y;	/* height difference: P[a] - P[c] in dir. u. */
	X	tPointi	ur;	/* u rotated cw by pi/2 */
	X
	X	printf("\n==>Extreme: u = "); PrintPoint(u); putchar('\n');
	X	ur[0] = u[1];	ur[1] = -u[0];
	X	a = 0; b = 0;
	X	do {
	X		c = Midway( a, b, n );
	X		printf("Extreme: a,b,c=%d\t%d\t%d\n", a, b, c);
	X
	X		/* Compute basic vectors and dots. */
	X		a1 = ( a + 1 ) % n;
	X		SubVec( P[a1], P[a], A );
	X		c1 = ( c + 1 ) % n;
	X		SubVec( P[c1], P[c], C );
	X		c1 = ( c + ( n-1 ) ) % n;
	X		SubVec( P[c], P[c1], B );
	X		printf("Extreme: A,C,B: ");
	X		PrintPoint( A ); putchar('\t');
	X		PrintPoint( C ); putchar('\t');
	X		PrintPoint( B ); putchar('\n');
	X		Adot = Dot(u,A);
	X		Cdot = Dot(u,C);
	X		Bdot = Dot(u,B);
	X		y = Dot(u,P[a]) - Dot(u,P[c]);
	X		printf("Extreme:\tAdot=%d\t", Adot);
	X		printf("Cdot=%d\t", Cdot);
	X		printf("Bdot=%d\t", Bdot);
	X		printf("y = %d\n", y);
	X
	X		/* Termination conditions */
	X		/* If either A or C points left of u, then at extreme. */
	X		if ( (Adot == 0) && (Dot(ur,A) < 0) ) {
	X			printf("Extreme: A points left, so return a=%d\n",a);
	X			return	a;
	X		}
	X		if ( (Cdot == 0) && (Dot(ur,C) < 0) ) {
	X			printf("Extreme: C points left, so return c=%d\n",c);
	X			return	c;
	X		}
	X		/* From here on, can assume that zero dots put point on bot */
	X
	X		if ( (Cdot < 0) && (Bdot > 0) ) {
	X			printf("Extreme: unique extreme c=%d\n", c);
	X			return	c;
	X		}
	X		if (a == c) {
	X			printf("Extreme: a=c=%d: choosing a or b=a+1\n", a);
	X			if (Adot > 0)
	X				return	b;
	X			else
	X				return	a;
	X		}
	X
	X		/* Halving interval */
	X		if ( (Adot >= 0) && (Cdot <= 0) )
	X			b = c;		/* new: (a,c) */
	X		else if ( (Adot <= 0) && (Cdot >= 0) )
	X			a = c;		/* new: (c,b) */
	X		else if ( (Adot > 0) && (Cdot > 0) ) {
	X			if ( y > 0 )
	X				b = c;	/* new: (a,c) */
	X			else
	X				a = c;	/* new: (c,b) */
	X		}
	X		else if ( (Adot < 0) && (Cdot < 0) ) {
	X			if ( y < 0 )
	X				b = c;	/* new: (a,c) */
	X			else
	X				a = c;	/* new: (c,b) */
	X		}
	X		else {
	X			printf("Error in Extreme: shouldn't reach here\n");
	X			exit(1);
	X		}
	X	}
	X	while	( TRUE );
	X}
	X/*
	X        Puts a - b into c.
	X*/
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X        int     i;
	X        
	X        for( i = 0; i < DIM; i++ )
	X                c[i] = a[i] - b[i];
	X}
	X
	X/*
	X        Returns the dot product of the two input vectors.
	X*/
	Xint     Dot( tPointi a, tPointi b )
	X{
	X        int     i;
	X        int     sum;
	X	sum = 0;
	X/*
	X	printf("Dot: a, b =\n");
	X	PrintPoint(a); putchar('\n');
	X	PrintPoint(b); putchar('\n');
	X*/
	X
	X        for( i = 0; i < DIM; i++ ) {
	X		/*printf("before: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", 
	X			i, sum, a[i], b[i]);*/
	X                sum += a[i] * b[i];
	X		/*printf("after: i=%d, sum=%d, a[i]=%d, b[i]=%d\n", 
	X			i, sum, a[i], b[i]);*/
	X	}
	X	/*printf("Dot: returning %d\n", sum);*/
	X        
	X        return  sum;
	X}
	X
	X/*
	X	Implements a = b, assignment of points/vectors.
	X	Assignment between arrays is not possible in C.
	X*/
	Xvoid	PointAssign( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < DIM; i++ )
	X		a[i] = b[i];
	X}
	X
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X/*
	X	Reads in the coordinates of the vertices of a polygon from stdin,
	X	puts them into P, and returns n, the number of vertices.
	X	Formatting conventions: etc.
	X*/
	Xint	ReadPoly( tPolygoni P )
	X{
	X	int	n = 0;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   l   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]);
	X}
	X
SHAR_EOF
if test 5406 -ne "`wc -c < 'extreme.c'`"
then
	echo shar: error transmitting "'extreme.c'" '(should have been 5406 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'graham.c'" '(3937 characters)'
if test -f 'graham.c'
then
	echo shar: will not over-write existing file "'graham.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'graham.c'
	X/* 
	X   Written by Joseph O'Rourke.
	X   Graham's algorithm for hull of 2-dim points.
	X   See "Computational Geometry in C."
	X*/
	X#include	<stdio.h>
	X#include	<math.h>
	X#include	<stdlib.h>
	X/* #include	<malloc.h> */
	X
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X#define PMAX    1000            /* Max # of pts in polygon */
	X   
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xtypedef	struct tCell *tStack;
	Xstruct tCell {
	X	int	i;
	X	tStack	next;
	X};
	X
	XtStack	Pop( tStack p );
	XtStack	Push( tStack t, tStack p );
	XtStack	GetCell( void );
	Xvoid	PrintStack( tStack t );
	XtStack	GetPush( int i, tStack t );
	XtStack	Graham( tPolygoni P );
	X
	X/*----------from tri.c-------------*/
	Xint     Area2( tPointi a, tPointi b, tPointi c );
	Xbool 	xor( bool x, bool y );
	Xbool    Left( tPointi a, tPointi b, tPointi c );
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c );
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c );
	Xint     Dot( tPointi a, tPointi b );
	Xint     Length2( tPointi p );
	Xvoid	PointAssign( tPointi a, tPointi b );
	Xint	ReadPoly( tPolygoni P );
	Xvoid	PrintPoly( int n, tPolygoni P );
	Xvoid	PrintPoint( tPointi p );
	X/*----------from tri.c-------------*/
	Xint	Compare( tPointi *p1, tPointi *p2 );
	Xvoid	FindLowest( int n, tPolygoni P );
	X
	Xstatic	tPolygoni	P;	/* global so compare can access it. */
	Xint	n;
	Xmain()
	X{
	X	tStack	top;
	X
	X	n = ReadPoly( P );
	X	FindLowest( n, P );
	X	PrintPoly( n, P );
	X	qsort(
	X		(void *) P[1],		/* base: pointer to 1st elem */
	X		n-1,			/* count: # of elems */
	X		sizeof( tPointi ),	/* size of each elem */
	X		Compare			/* -1,0,+1 compare fnc */
	X	);
	X	PrintPoly( n, P );
	X
	X	top = NULL;
	X	top = Graham( P );
	X	printf("Hull:\n");
	X	PrintStack( top );
	X}
	X
	X/*
	X	FindLowest finds the rightmost lowest point and swaps with 0-th.
	X	The lowest point has the min y-coord, and amongst those, the
	X	max x-coord: so it is rightmost among the lowest.
	X*/
	Xvoid	FindLowest( int n, tPolygoni P )
	X{
	X	int 	i;
	X	int	m;	/* Index of lowest so far. */
	X	tPointi	low;	/* To hold point when swapping. */
	X
	X	m = 0;
	X	for ( i = 1; i < n; i++ )
	X		if ( (P[i][Y] < P[m][Y]) ||
	X		    ((P[i][Y] == P[m][Y]) && (P[i][X] > P[m][X])) ) 
	X			m = i;
	X	/* swap */
	X	PointAssign( low, P[m] );
	X	PointAssign( P[m], P[0] );
	X	PointAssign( P[0], low  );
	X}
	X/*
	X	Compare: returns -1,0,+1 if p1 < p2, =, or > respectively;
	X	here "<" means smaller angle.  Follows the conventions of qsort.
	X*/
	Xint	Compare( tPointi *p1, tPointi *p2 )
	X{
	X	int 	a;		/* area */
	X	tPointi	r1, r2;		/* ri = pi - p0 */
	X	int	len1, len2;	/* length of r1 & r2 */
	X
	X	a = Area2( P[0], *p1, *p2 );
	X	if (a > 0)
	X		return -1;
	X	else if (a < 0)
	X		return 1;
	X	else {
	X		SubVec( *p1, P[0], r1 );
	X		SubVec( *p2, P[0], r2 );
	X		len1 = Length2( r1 );
	X		len2 = Length2( r2 );
	X		printf("compare:\n");
	X		PrintPoint(r1);PrintPoint(r2);
	X		printf("\nlen1=%d, len2=%d\n", len1, len2);
	X		return (len1 < len2) ? -1 : (len1 > len2) ? 1 : 0;
	X	}
	X}
	X
	X#include	"stack.c"
	X/*
	X	Get a cell, fill it with i, and push it onto the stack.
	X	Return pointer to stack top.
	X*/
	XtStack	GetPush( int i, tStack top )
	X{
	X	tStack	p;
	X
	X	p = GetCell();
	X	p->i = i;
	X	return Push( top, p );
	X}
	X
	X/*
	X	Performs the Graham scan on an array of angularly sorted points P.
	X*/
	XtStack	Graham( tPolygoni P )
	X{
	X	int	i;
	X	tStack	top;
	X
	X	/* Initialize stack to (P[n-1], P[0]). */
	X	top = NULL;
	X	top = GetPush ( n-1, top );
	X	top = GetPush ( 0, top );
	X
	X	/* Bottom two elements will never be removed. */
	X	i = 1;
	X
	X	while ( i < n ) {
	X		printf("Stack at top of while loop, i=%d:\n", i);
	X		PrintStack( top );
	X		if ( Left( P[ (top->next)->i ], P[ top->i ], P[ i ] ) ) {
	X			top = GetPush ( i, top );
	X			i++;
	X		}
	X		else 	top = Pop( top );
	X		printf("Stack at bot of while loop, i=%d:\n", i);
	X		PrintStack( top );
	X		putchar('\n');
	X	}
	X
	X	/* P[n-1] pushed twice, so pop it off. */
	X	return	Pop(top);
	X
	X}
	X
	X#include	"vector.c"
SHAR_EOF
if test 3937 -ne "`wc -c < 'graham.c'`"
then
	echo shar: error transmitting "'graham.c'" '(should have been 3937 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'i.poly.18'" '(94 characters)'
if test -f 'i.poly.18'
then
	echo shar: will not over-write existing file "'i.poly.18'"
else
sed 's/^	X//' << \SHAR_EOF > 'i.poly.18'
	X0	0
	X10	7
	X12	3
	X20	8
	X13	17
	X10	12
	X12	14
	X13	11
	X7	11
	X6	14
	X10	15
	X6	18
	X-1	15
	X1	13
	X4	14
	X5	10
	X-2	9
	X5	5
SHAR_EOF
if test 94 -ne "`wc -c < 'i.poly.18'`"
then
	echo shar: error transmitting "'i.poly.18'" '(should have been 94 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'inpoly.c'" '(2838 characters)'
if test -f 'inpoly.c'
then
	echo shar: will not over-write existing file "'inpoly.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'inpoly.c'
	X/*
	X	inpoly.c
	X	Written by Joseph O'Rourke
	X	orourke@sophia.smith.edu
	X*/
	X#include	<stdio.h>
	X#include	<math.h>
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X#define PMAX    1000            /* Max # of pts in polygon */
	X   
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xvoid	PrintPoly( int n, tPolygoni P );
	Xvoid	PrintPoint( tPointi p );
	Xbool	InPoly( tPointi q, tPolygoni P, int n );
	X
	Xmain()
	X{
	X	int		n;
	X	tPolygoni	P;
	X	tPointi		q;
	X
	X	n = ReadPoly( P );
	X	/* Last point is the point q */
	X	q[0] = P[n-1][0];
	X	q[1] = P[n-1][1];
	X	n = n - 1;
	X	PrintPoly( n, P );
	X	printf("In = %d\n", InPoly( q, P, n ));
	X}
	X
	X/*
	X	Returns true if q is inside polygon P.
	X*/
	Xbool	InPoly( tPointi q, tPolygoni P, int n )
	X{
	X	int	i, i1;		/* point index; i1 = i-1 mod n */
	X	int	d;		/* dimension index */
	X	double	x;		/* x intersection of e with ray */
	X	int	crossings = 0;	/* number of edge/ray crossings */
	X
	X	printf("\n==>In: u = "); PrintPoint(q); putchar('\n');
	X
	X	/* Shift so that q is the origin. */
	X	for( i = 0; i < n; i++ ) {
	X		for( d = 0; d < DIM; d++ )
	X			P[i][d] = P[i][d] - q[d];
	X	}
	X
	X	/* For each edge e=(i-1,i), see if crosses ray. */
	X	for( i = 0; i < n; i++ ) {
	X		i1 = ( i + n - 1 ) % n;
	X		printf("e=(%d,%d)\t", i1, i);
	X		/* if e straddles the x-axis... */
	X		if( ( ( P[i] [Y] > 0 ) && ( P[i1][Y] <= 0 ) ) ||
	X		    ( ( P[i1][Y] > 0 ) && ( P[i] [Y] <= 0 ) ) ) {
	X			/* e straddles ray, so compute intersection with ray. */
	X			x = (P[i][X] * P[i1][Y] - P[i1][X] * P[i][Y])
	X				/ (double)(P[i1][Y] - P[i][Y]);
	X			printf("straddles: x = %g\t", x);
	X			/* crosses ray if strictly positive intersection. */
	X			if (x > 0) crossings++;
	X		}
	X		printf("crossings=%d\n", crossings);
	X	}
	X	/* q inside if an odd number of crossings. */
	X	if( (crossings % 2) == 1 )
	X		return	TRUE;
	X	else	return	FALSE;
	X}
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X/*
	X	Reads in the coordinates of the vertices of a polygon from stdin,
	X	puts them into P, and returns n, the number of vertices.
	X	Formatting conventions: etc.
	X*/
	Xint	ReadPoly( tPolygoni P )
	X{
	X	int	n = 0;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d\n", i, P[i][0], P[i][1]);
	X}
	X
SHAR_EOF
if test 2838 -ne "`wc -c < 'inpoly.c'`"
then
	echo shar: error transmitting "'inpoly.c'" '(should have been 2838 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'int.c'" '(8242 characters)'
if test -f 'int.c'
then
	echo shar: will not over-write existing file "'int.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'int.c'
	X#include        <stdio.h>
	X#include        <math.h>
	X#define ERROR   1
	X#define X       0
	X#define Y       1
	Xtypedef enum { FALSE, TRUE }    bool;
	Xtypedef enum { Pin, Qin, Unknown }       tInFlag;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X#define PMAX    1000            /* Max # of pts in polygon */
	X
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xint     Area2( tPointi a, tPointi b, tPointi c );
	Xvoid    SubVec( tPointi a, tPointi b, tPointi c );
	Xint     Dot( tPointi a, tPointi b );
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c );
	Xbool    Left( tPointi a, tPointi b, tPointi c );
	Xvoid    PrintPoly( int n, tPolygoni P );
	Xvoid    PrintPoint( tPointi p );
	Xvoid    PrintPointd( tPointd p );
	Xvoid    ConvexIntersect( tPolygoni P, tPolygoni Q, int n, int m );
	XtInFlag InOut( tPointd p, tInFlag inflag, bool aHB, bool bHA );
	Xint     Advance( int a, int *aa, int n, bool inside, tPointi v );
	Xbool    Intersection( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p );
	Xbool    EqPoint( tPointi a, tPointi b );
	X
	X	int     	n, m;
	X	tPolygoni	P, Q;
	X
	Xmain()
	X{
	X
	X        n = ReadPoly( P );
	X        m = ReadPoly( Q );
	X        ConvexIntersect( P, Q, n, m );
	X}
	X
	Xvoid    ConvexIntersect( tPolygoni P, tPolygoni Q, int n, int m )
	X                                /* P has n vertices, Q has m vertices. */
	X{
	X        int     a, b;           /* indices on P and Q (resp.) */
	X        int     a1, b1;         /* a-1, b-1 (resp.) */
	X        tPointi A, B;           /* directed edges on P and Q (resp.) */
	X        int     cross;          /* A x B */
	X        bool    bHA, aHB;       /* b in H(A); a in H(b). */
	X        tPointi Origin = {0,0}; /* (0,0) */
	X        tPointd p;              /* double point of intersection */
	X        tInFlag inflag;         /* {Pin, Qin, Unknown}:
	X                                   which polygon is known to be inside */
	X        int     i;              /* loop counter */
	X	int	aa, ba;		/* # advances on a & b indices 
	X				   (from first intersection) */
	X	bool    Reset = FALSE;  /* have the advance counters ever been reset */
	X
	X        /* Initialize variables. */
	X        a = 0; b = 0; aa = 0; ba = 0;
	X        i = 0;
	X        inflag = Unknown;
	X
	X        do {
	X                /* Computations of key variables. */
	X                a1 = (a + n - 1) % n;
	X                b1 = (b + m - 1) % m;
	X
	X                SubVec( P[a], P[a1], A );
	X                SubVec( Q[b], Q[b1], B );
	X
	X                cross = Area2( Origin, A, B );
	X                bHA = Left( P[a1], P[a], Q[b] );
	X                aHB = Left( Q[b1], Q[b], P[a] );
	X
	X                /* If A & B intersect, update inflag. */
	X                if ( Intersection( P[a1], P[a], Q[b1], Q[b], p ) ) {
	X			if ( inflag == Unknown  && !Reset ) {
	X				aa = ba = 0; 
	X				Reset = TRUE;
	X			}
	X                        inflag = InOut( p, inflag, aHB, bHA );
	X		}
	X                /* Advance rules. */
	X		if ( (cross == 0) && !bHA && !aHB ) {
	X			if ( inflag == Pin )
	X                                b = Advance( b, &ba, m, inflag == Qin, Q[b] );
	X			else
	X                                a = Advance( a, &aa, n, inflag == Pin, P[a] );
	X
	X                }
	X                else if ( cross >= 0 )
	X                        if ( bHA )
	X                                a = Advance( a, &aa, n, inflag == Pin, P[a] );
	X                        else {
	X                                b = Advance( b, &ba, m, inflag == Qin, Q[b] );
	X                }
	X                else /* if ( cross < 0 ) */{
	X                        if ( aHB )
	X                                b = Advance( b, &ba, m, inflag == Qin, Q[b] );
	X                        else
	X                                a = Advance( a, &aa, n, inflag == Pin, P[a] );
	X                }
	X
	X	/* Quit when both adv. indices have cycled, or one has cycled twice. */
	X        } while ( ((aa < n) || (ba < m)) && (aa < 2*n) && (ba < 2*m) );
	X
	X        /* Deal with special cases: not implemented here. */
	X	if ( inflag == Unknown) 
	X		printf("The boundaries of P and Q do not cross:\n");
	X}
	X
	X/*	
	X	Prints out the double point of intersection, and toggles
	X	in/out flag.
	X*/
	XtInFlag InOut( tPointd p, tInFlag inflag, bool aHB, bool bHA )
	X{
	X        PrintPointd( p ); putchar('\n');
	X
	X        /* Update inflag. */
	X        if ( aHB )
	X                return   Pin;
	X        else if ( bHA )
	X                return   Qin;
	X	else 	return	inflag;
	X}
	X/*
	X	Advances and prints out an inside vertex if appropriate.
	X*/
	Xint     Advance( int a, int *aa, int n, bool inside, tPointi v )
	X{
	X        if ( inside ) {
	X                PrintPoint( v ); putchar('\n');
	X        }
	X	(*aa)++;
	X        return  (a+1) % n;
	X}
	Xbool    Intersection( tPointi a, tPointi b, tPointi c, tPointi d, tPointd p )
	X{
	X        double  s, t;   /* The two parameters of the parametric eqns. */
	X        double  denom;  /* Denoninator of solutions. */
	X
	X        denom =
	X                a[0] * ( d[1] - c[1] ) +
	X                b[0] * ( c[1] - d[1] ) +
	X                d[0] * ( b[1] - a[1] ) +
	X                c[0] * ( a[1] - b[1] );
	X
	X        /* If denom is zero, then the line segments are parallel. */
	X        /* In this case, return false even though the segments might overlap. */
	X        if (denom == 0.0)
	X                return  FALSE;
	X
	X        s = (
	X                a[0] * ( d[1] - c[1] ) +
	X                c[0] * ( a[1] - d[1] ) +
	X                d[0] * ( c[1] - a[1] )
	X             ) / denom;
	X        t = -(
	X                a[0] * ( c[1] - b[1] ) +
	X                b[0] * ( a[1] - c[1] ) +
	X                c[0] * ( b[1] - a[1] )
	X             ) / denom;
	X
	X        p[0] = a[0] + s * ( b[0] - a[0] );
	X        p[1] = a[1] + s * ( b[1] - a[1] );
	X
	X        if (    (0.0 <= s) && (s <= 1.0) &&
	X                (0.0 <= t) && (t <= 1.0)
	X           )
	X                return  TRUE;
	X	else	return	FALSE;
	X}
	X
	X/*
	X	Implements a = b, assignment of points/vectors.
	X	Assignment between arrays is not possible in C.
	X*/
	Xvoid	PointAssign( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < DIM; i++ )
	X		a[i] = b[i];
	X}
	X
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	Xvoid    PrintPointd( tPointd p )
	X{
	X        int     i;
	X
	X        putchar('(');
	X        for ( i = 0; i < DIM; i++ ) {
	X                printf("%5.2f", p[i]);
	X                if ( i != DIM-1 ) putchar(',');
	X        }
	X        putchar(')');
	X}
	X
	X/*
	X	Reads in the coordinates of the vertices of a polygon from stdin,
	X	puts them into P, and returns n, the number of vertices.
	X	Formatting conventions: etc.
	X*/
	Xint	ReadPoly( tPolygoni P )
	X{
	X	int	n = 0;
	X	int	nin;
	X
	X	scanf("%d", &nin);
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < nin) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   l   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d%4d\n", i, P[i][0], P[i][1]);
	X}
	X
	X/*
	X	Returns true iff c is strictly to the left of the directed
	X	line through a to b.
	X*/
	Xbool    Left( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) > 0;
	X}
	X
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) >= 0;
	X}
	X
	Xbool    Collinear( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) == 0;
	X}
	X/*
	X        Puts a - b into c.
	X*/
	Xvoid    SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X        int     i;
	X
	X        for( i = 0; i < DIM; i++ )
	X                c[i] = a[i] - b[i];
	X}
	X
	X/*
	X        Returns twice the signed area of the triangle determined by a,b,c,
	X        positive if a,b,c are oriented ccw, and negative if cw.
	X*/
	Xint     Area2( tPointi a, tPointi b, tPointi c )
	X{
	X        return
	X                a[0] * b[1] - a[1] * b[0] +
	X                a[1] * c[0] - a[0] * c[1] +
	X                b[0] * c[1] - c[0] * b[1];
	X}
	X
	Xbool    EqPoint( tPointi a, tPointi b )
	X{
	X        int     i;
	X        for ( i = 0; i < DIM; i++ )
	X                if ( a[i] != b[i])	
	X			return	FALSE;
	X	return	TRUE;
	X}
SHAR_EOF
if test 8242 -ne "`wc -c < 'int.c'`"
then
	echo shar: error transmitting "'int.c'" '(should have been 8242 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'macros.h'" '(956 characters)'
if test -f 'macros.h'
then
	echo shar: will not over-write existing file "'macros.h'"
else
sed 's/^	X//' << \SHAR_EOF > 'macros.h'
	X/*====================================================================
	X    macros.h
	X 
	X 	macros used to access data structures and perform quick tests.
	X
	X  ====================================================================*/
	X
	X/* general-purpose macros */
	X#define SWAP(t,x,y)	{ t = x; x = y; y = t; }
	X
	Xchar *malloc();
	X
	X#define NEW(p,type)	if ((p=(type *) malloc (sizeof(type))) == NULL) {\
	X				printf ("Out of Memory!\n");\
	X				exit(0);\
	X			}
	X
	X#define FREE(p)		if (p) { free ((char *) p); p = NULL; }
	X
	X
	X#define ADD( head, p )  if ( head )  { \
	X				p->next = head->next; \
	X				p->prev = head; \
	X				head->next = p; \
	X				p->next->prev = p; \
	X			} \
	X			else { \
	X				head = p; \
	X				head->next = head->prev = p; \
	X			}
	X
	X#define DELETE( head, p )   if ( head )  { \
	X				if ( head == head->next ) \
	X					head = NULL;  \
	X				else if ( p == head ) \
	X					head = head->next; \
	X				p->next->prev = p->prev;  \
	X				p->prev->next = p->next;  \
	X				FREE( p ); \
	X			} 
	X
SHAR_EOF
if test 956 -ne "`wc -c < 'macros.h'`"
then
	echo shar: error transmitting "'macros.h'" '(should have been 956 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'quadratic.c'" '(10343 characters)'
if test -f 'quadratic.c'
then
	echo shar: will not over-write existing file "'quadratic.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'quadratic.c'
	X/* 
	X   Written by Joseph O'Rourke.
	X   Last modified February 17, 1991.
	X   Code to triangulate a simple polygon.
	X   Quadratic algorithm.
	X   See "Computational Geometry in C."
	X*/
	X#include	<stdio.h>
	X#include	<math.h>
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X#define PMAX    1000            /* Max # of pts in polygon */
	X   
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xint     Area2( tPointi a, tPointi b, tPointi c );
	Xint     AreaPoly2( int n, tPolygoni P );
	Xbool 	xor( bool x, bool y );
	Xbool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d );
	Xbool    Left( tPointi a, tPointi b, tPointi c );
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c );
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c );
	Xint     Dot( tPointi a, tPointi b );
	Xint     Length2( tPointi p );
	Xbool    Between( tPointi a, tPointi b, tPointi c );
	Xbool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d );
	Xbool	Diagonalie( int i, int j, int n, tPolygoni P );
	Xbool	Diagonal( int i, int j, int n, tPolygoni P );
	Xvoid	PointAssign( tPointi a, tPointi b );
	Xvoid	Triangulate2( int n, tPolygoni P, bool Ear[], int labels[]);
	Xvoid	Triangulate( int n, tPolygoni P );
	Xvoid	TriRecurse( int n, tPolygoni P, int labels[] );
	Xbool	InCone( int i, int j, int n, tPolygoni P );
	Xvoid	ClipEar1( int i, int n, tPolygoni P, int labels[] );
	Xvoid	ClipEar2( int i, int n, tPolygoni P, bool Ear[], int labels[] );
	Xvoid	ClipEar( int i, int n, tPolygoni P );
	Xint	ReadPoly( tPolygoni P );
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] );
	Xvoid	PrintPoint( tPointi p );
	Xvoid    UpdateArray( bool Ear[PMAX], int j, int n);
	X
	Xmain()
	X{
	X	int	n;
	X	tPolygoni	P;
	X
	X	n = ReadPoly( P );
	X	Triangulate( n, P );
	X}
	X
	X/* 
	X        Returns twice the signed area of the triangle determined by a,b,c.
	X        The area is positive if a,b,c are oriented ccw, negative if cw,
	X	and zero if the points are collinear.
	X*/
	Xint     Area2( tPointi a, tPointi b, tPointi c )
	X{
	X        return
	X                a[0] * b[1] - a[1] * b[0] +
	X                a[1] * c[0] - a[0] * c[1] +
	X                b[0] * c[1] - c[0] * b[1];
	X}
	X
	X/*      
	X        Returns twice the area of polygon P.
	X*/
	Xint     AreaPoly2( int n, tPolygoni P )
	X{
	X        int     i;
	X        int     sum = 0;        /* Partial area sum */
	X
	X        for (i = 1; i < n-1; i++)
	X                sum += Area2( P[0], P[i], P[i+1] );
	X        return  sum;
	X}
	X
	X/*
	X	Exclusive or: true iff exactly one argument is true.
	X	The arguments are negated to ensure that they are 0/1
	X	values.  Then the bitwise xor operator may apply.
	X	(This idea is due to Michael Baldwin.)
	X*/
	Xbool 	xor( bool x, bool y )
	X{
	X	return	!x ^ !y;
	X}
	X/*
	X	Returns true iff ab properly intersects cd: they share
	X	a point interior to both segments.  The properness of the
	X	intersection is ensured by using strict leftness.
	X*/
	Xbool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d )
	X{
	X        return
	X                xor( Left(a,b,c), Left(a,b,d) )
	X        &&      xor( Left(c,d,a), Left(c,d,b) );
	X}
	X
	X/*
	X	Returns true iff c is strictly to the left of the directed
	X	line through a to b.
	X*/
	Xbool    Left( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) > 0;
	X}
	X
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) >= 0;
	X}
	X
	Xbool    Collinear( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) == 0;
	X}
	X
	X/*
	X        Puts a - b into c.
	X*/
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X        int     i;
	X        
	X        for( i = 0; i < DIM; i++ )
	X                c[i] = a[i] - b[i];
	X}
	X
	X/*
	X        Returns the dot product of the two input vectors.
	X*/
	Xint     Dot( tPointi a, tPointi b )
	X{
	X        int     i;
	X        int     sum = 0;
	X
	X        for( i = 0; i < DIM; i++ )
	X                sum += a[i] * b[i];
	X        
	X        return  sum;
	X}
	X
	X/*
	X        Returns the square of the length of the vector p.
	X*/
	Xint     Length2( tPointi p )
	X{
	X        return  Dot( p, p );
	X}
	X/*
	X        Returns T iff point c lies on the closed segement ab.
	X        First checks that c is collinear with a and b.
	X*/
	Xbool    Between( tPointi a, tPointi b, tPointi c )
	X{
	X        tPointi ba, ca;
	X        if ( ! Collinear( a, b, c ) )
	X                return  FALSE;
	X        SubVec( b, a, ba );
	X        SubVec( c, a, ca );
	X        return
	X                Dot( ca, ba ) >= 0
	X        &&      Length2( ca ) <= Length2( ba );
	X}
	X
	X/*
	X        Returns true iff segments ab and cd intersect,
	X        properly or improperly.
	X*/
	Xbool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d )
	X{
	X        if ( IntersectProp( a, b, c, d ) )
	X                return  TRUE;
	X        else if (       Between( a, b, c )
	X                     || Between( a, b, d )
	X                     || Between( c, d, a )
	X                     || Between( c, d, b )
	X                )
	X                return  TRUE;
	X        else    return  FALSE;
	X}
	X
	X/*
	X	Returns T iff (v_i, v_j) is a proper internal *or* external
	X	diagonal of P, *ignoring edges incident to v_i and v_j*.
	X*/
	Xbool	Diagonalie( int i, int j, int n, tPolygoni P )
	X{
	X	int	k;
	X	int	k1;
	X
	X	/* For each edge (k,k+1) of P */
	X	for ( k = 0; k < n; k++ ) {
	X		k1 = (k+1) % n;
	X		/* Skip edges incident to i or j */
	X		if ( ! (
	X				( k == i ) || ( k1 == i )
	X			||	( k == j ) || ( k1 == j )
	X		       )
	X		   )
	X			if ( Intersect( P[i], P[j], P[k], P[k1] ) )
	X				return	FALSE;
	X	}
	X	return	TRUE;
	X}
	X
	X/*
	X	Implements a = b, assignment of points/vectors.
	X	Assignment between arrays is not possible in C.
	X*/
	Xvoid	PointAssign( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < DIM; i++ )
	X		a[i] = b[i];
	X}
	X
	X
	Xvoid	TriRecurse( int n, tPolygoni P, int labels[] )
	X{
	X	int	i, i1, i2;
	X
	X	printf("TriRecurse: n = %d\n", n);
	X	/*PrintPoly( n, P, labels );*/
	X	if ( n > 3 )
	X		for ( i = 0; i < n; i++ ) {
	X			i1 = (i+1) % n;
	X			i2 = (i+2) % n;
	X			if ( Diagonal( i, i2, n, P ) ) {
	X				printf("%d %d\n", labels[i], labels[i2]);
	X				ClipEar1( i1, n, P, labels );
	X				TriRecurse( n-1, P, labels );
	X				break;
	X			}
	X		}
	X}
	X
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X
	X/*
	X	Prints out n-3 diagonals (as pairs of integer indices)
	X	which form a triangulation of P.
	X	This function initializes the data structures, and calls
	X	Triangulate2 to clip off the ears one by one.
	X*/
	Xvoid	Triangulate( int n, tPolygoni P )
	X{
	X	tPolygoni	Pt;
	X        bool    Ear[PMAX];	/* Ear[i] true iff (i,i+1,i+2) is an ear */
	X	int	labels[PMAX];
	X	int	i, i1, i2;
	X
	X	/* Copy points into temporary array & initialize labels. */
	X	for ( i = 0; i < n; i++ ){
	X		PointAssign( Pt[i], P[i] );
	X		labels[i] = i;
	X	}
	X
	X       	/* Ear[i] is true iff (i,i+1,i+2) is an ear. */
	X	/* Initialize Ear[] for all i. */
	X	if ( n > 3 )
	X		for ( i = 0; i < n; i++ ) {
	X			i1 = (i+1) % n;
	X			i2 = (i+2) % n;
	X			Ear[i] = Diagonal( i, i2, n, P ); 
	X                        }
	X
	X	Triangulate2( n, Pt, Ear, labels );
	X}
	X
	X
	X/*
	X	Triangulate2 is an O(n^2) triangulation function
	X	(or it would be if the array squashing were replaced by
	X	pointer movements).
	X	Assumes all arrays (including Ear) have been initialized.
	X*/
	X
	Xvoid	Triangulate2( int n, tPolygoni P , bool Ear[], int labels[] )
	X{
	X	int	im1, i, i1, i2, i3;	/* i-1,i,i+1,i+2,i+3 */
	X
	X	/* Each step of outer loop removes one ear. */
	X        while ( n > 3 ) {     
	X		/* Inner loop searches for an ear. */
	X                for (i = 0; i < n; i++)
	X                        if (Ear[i]) {
	X				/* Ear found. */
	X				im1 = (i-1+n) % n;
	X				i1 = (i+1) % n;  
	X				i2 = (i+2) % n;  
	X				i3 = (i+3) % n;  
	X
	X				/* (i,i+2) is a diagonal */
	X				printf("%d %d\n", labels[i], labels[i2]);
	X				
	X				/* Replace the two entries for segments
	X				   incident to i+1:
	X				   (i-1,i+1) ==> (i-1,i+2)
	X				   (i+1,i+3) ==> (i,i+3)
	X				*/
	X                                Ear[im1] = Diagonal( im1, i2, n, P);
	X                                Ear[i]   = Diagonal( i,   i3, n, P);
	X
	X				/* Squash out the i1 entry in all arrays */
	X	        		ClipEar2( i1, n, P, Ear, labels );
	X                                n = n - 1;
	X
	X                                break;	/* out of inner loop */
	X                                }
	X                }
	X}
	X
	Xvoid	ClipEar2( int i, int n, tPolygoni P, bool Ear[], int labels[] )
	X{
	X	int	k;
	X	
	X	for ( k = i; k < n-1; k++ ) {
	X		PointAssign( P[k], P[k+1] );
	X		Ear[k] = Ear[k+1];
	X		labels[k] = labels[k+1];
	X	}
	X}
	X
	X
	X/*
	X	Returns true iff the diagonal (i,j) is striclty internal to the 
	X	polygon P in the neighborhood of the i endpoint.
	X*/
	Xbool	InCone( int i, int j, int n, tPolygoni P )
	X{
	X	int	i1;	/* i + 1 */
	X	int	in1;	/* i - 1 */
	X
	X	i1  = (i + 1)     % n;
	X	in1 = (i + n - 1) % n;
	X
	X	/* If P[i] is a convex vertex [ i+1 left of (i-1,i) ]. */
	X	if ( Left( P[in1], P[i], P[i1] ) )
	X		return	Left( P[i], P[j], P[in1] )
	X		     &&	Left( P[j], P[i], P[i1]  );
	X
	X	/* Assume (i-1,i,i+1) not collinear. */
	X	/* else P[i] is reflex. */
	X	else
	X		return	!(	LeftOn( P[i], P[j], P[i1]  )
	X		            &&	LeftOn( P[j], P[i], P[in1] ) );
	X}
	X
	X/*
	X	Returns T iff (v_i, v_j) is a proper internal
	X	diagonal of P.
	X*/
	Xbool	Diagonal( int i, int j, int n, tPolygoni P )
	X{
	X	return	InCone( i, j, n, P ) && Diagonalie( i, j, n, P );
	X}
	X
	X
	X/*
	X	Removes P[i] by copying P[i+1]...P[n-1] left one index.
	X*/
	Xvoid	ClipEar1( int i, int n, tPolygoni P, int labels[] )
	X{
	X	int	k;
	X	
	X	for ( k = i; k < n-1; k++ ) {
	X		PointAssign( P[k], P[k+1] );
	X		labels[k] = labels[k+1];
	X	}
	X}
	X
	Xvoid	ClipEar( int i, int n, tPolygoni P )
	X{
	X	int	k;
	X	
	X	for ( k = i; k < n-1; k++ )
	X		PointAssign( P[k], P[k+1] );
	X}
	X
	X/*
	X	Reads in the coordinates of the vertices of a polygon from stdin,
	X	puts them into P, and returns n, the number of vertices.
	X	Formatting conventions: etc.
	X*/
	Xint	ReadPoly( tPolygoni P )
	X{
	X	int	n = 0;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   l   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d%4d\n", i, labels[i], P[i][0], P[i][1]);
	X}
SHAR_EOF
if test 10343 -ne "`wc -c < 'quadratic.c'`"
then
	echo shar: error transmitting "'quadratic.c'" '(should have been 10343 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'sphere.c'" '(1365 characters)'
if test -f 'sphere.c'
then
	echo shar: will not over-write existing file "'sphere.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'sphere.c'
	X/*----------------------------------------------------------------------
	X     Sphere                                   Susan L. Weller           
	X  ----------------------------------------------------------------------
	X         This program will generate a given number of points that lie
	X    on the unit sphere. The number of points is given on  the command line.
	X  ---------------------------------------------------------------------- */
	X#include <stdio.h>
	X#include <math.h>
	X
	X#define MAX_INT    	2147483647     /* Max int is 2^31 - 1 */
	X#define FUZZ 		0.000001
	X
	Xmain( argc, argv )
	Xint argc;
	Xchar *argv[];
	X{
	Xint n;
	Xdouble x, y, z, r;
	X
	X	srandom( (int) time((long *) 0 ) );
	X	if ( argc != 2 )
	X   	printf( "Usage: sphere <number of random points> \n" ),
	X   	exit(1);
	X
	X	n = atoi( argv[1] );            /* n is number of points */
	X
	X	while (n--) {
	X    	      do {
	X            	    /* the do-while will generate points until one
	X               	       is inside the unit sphere.               */
	X	    	    x = 2.0 * (double) random() / MAX_INT  - 1.0;
	X	    	    y = 2.0 * (double) random() / MAX_INT  - 1.0;
	X	    	    z = 2.0 * (double) random() / MAX_INT  - 1.0;
	X	    	    r = sqrt( x*x + y*y + z*z );
	X	      
	X    	      } while ( r <= FUZZ || r > 1.);
	X    	      printf ("%6d  %6d  %6d\n", (int) (100 * x/r), 
	X				      (int) (100 * y/r), (int) (100 * z/r) );
	X    	}
	X}
SHAR_EOF
if test 1365 -ne "`wc -c < 'sphere.c'`"
then
	echo shar: error transmitting "'sphere.c'" '(should have been 1365 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'stack.c'" '(861 characters)'
if test -f 'stack.c'
then
	echo shar: will not over-write existing file "'stack.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'stack.c'
	X/*----------------------Stack routines------------------*/
	X/*
	X	Pushes cell p on top of stack t, and returns new top.
	X*/
	XtStack	Push( tStack t, tStack p )
	X{
	X	p->next = t;
	X	return	p;
	X}
	X/*
	X	Pops off top elment of stack p, frees up the cell, and
	X	returns new top.
	X*/
	XtStack	Pop( tStack p )
	X{
	X	tStack	top;
	X
	X	top = p->next;
	X	free( (void *) p );
	X	return 	top;
	X}
	X/*
	X	Prints the stack, both point index and point coordinates.
	X*/
	Xvoid	PrintStack( tStack t )
	X{
	X	if (t) { 
	X		printf("i=%d:", t->i); PrintPoint( P[t->i] ); putchar('\n');
	X		PrintStack( t->next );
	X	}
	X}
	X/* 
	X	GetCell returns a pointer to a newly allocated piece of storage
	X	of type tCell, or exits if no space is available.
	X*/
	XtStack	GetCell( void )
	X{
	X	tStack	p;
	X
	X	p = (tStack) malloc( sizeof( struct tCell ) );
	X
	X	if (p == NULL) {
	X		printf("Error in GetCell: Out of memory!\n");
	X		exit(1);
	X	}
	X	else	return p;
	X}
SHAR_EOF
if test 861 -ne "`wc -c < 'stack.c'`"
then
	echo shar: error transmitting "'stack.c'" '(should have been 861 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'tri.c'" '(9289 characters)'
if test -f 'tri.c'
then
	echo shar: will not over-write existing file "'tri.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'tri.c'
	X/*
	X	tri.c
	X	Written by Joseph O'Rourke
	X	orourke@cs.smith.edu
	X	Triangulation code for "Computational Geometry in C."
	X        Assumes polygon vertices are entered in ccw order. 
	X*/
	X#include	<stdio.h>
	X#include	<math.h>
	X#define	ERROR	1
	X#define	X	0
	X#define	Y	1
	Xtypedef	enum { FALSE, TRUE }	bool;
	X
	X#define DIM     2               /* Dimension of points */
	Xtypedef int     tPointi[DIM];   /* type integer point */
	Xtypedef double  tPointd[DIM];   /* type double point */
	X
	X#define PMAX    1000            /* Max # of pts in polygon */
	Xtypedef tPointi tPolygoni[PMAX];/* type integer polygon */
	X
	Xint     Area2( tPointi a, tPointi b, tPointi c );
	Xint     AreaPoly2( int n, tPolygoni P );
	Xbool 	Xor( bool x, bool y );
	Xbool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d );
	Xbool    Left( tPointi a, tPointi b, tPointi c );
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c );
	Xbool    Collinear( tPointi a, tPointi b, tPointi c );
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c );
	Xint     Dot( tPointi a, tPointi b );
	Xint     Length2( tPointi p );
	Xbool    Between( tPointi a, tPointi b, tPointi c );
	Xbool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d );
	Xbool	Diagonalie( int i, int j, int n, tPolygoni P );
	Xbool	Diagonal( int i, int j, int n, tPolygoni P );
	Xvoid	PointAssign( tPointi a, tPointi b );
	Xvoid	Triangulate1( int n, tPolygoni P );
	Xvoid	Triangulate( int n, tPolygoni P );
	Xvoid	TriRecurse( int n, tPolygoni P, int labels[] );
	Xbool	InCone( int i, int j, int n, tPolygoni P );
	Xvoid	ClipEar1( int i, int n, tPolygoni P, int labels[] );
	Xvoid	ClipEar( int i, int n, tPolygoni P );
	Xint	ReadPoints( tPolygoni P );
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] );
	Xvoid	PrintPoint( tPointi p );
	X
	Xmain()
	X{
	X	int	n;
	X	tPolygoni	P;
	X
	X	n = ReadPoints( P );
	X	Triangulate1( n, P );
	X	Triangulate( n, P );
	X}
	X
	X/* 
	X        Returns twice the signed area of the triangle determined by a,b,c,
	X        positive if a,b,c are oriented ccw, and negative if cw.
	X*/
	Xint     Area2( tPointi a, tPointi b, tPointi c )
	X{
	X/*	The text has this:
	X        return
	X                a[0] * b[1] - a[1] * b[0] +
	X                a[1] * c[0] - a[0] * c[1] +
	X                b[0] * c[1] - c[0] * b[1];
	X
	X	The following computation is algebraically equivalent but
	X	uses four fewer multiplications.  It is obtained by shifting the
	X	coordinate system so that point a is the origin.
	X*/
	X
	X	return
	X		(b[0] - a[0]) * (c[1] - a[1]) -
	X		(c[0] - a[0]) * (b[1] - a[1]);
	X}
	X
	X/*      
	X        Returns twice the area of polygon P.
	X*/
	Xint     AreaPoly2( int n, tPolygoni P )
	X{
	X        int     i;
	X        int     sum = 0;        /* Partial area sum */
	X
	X        for (i = 1; i < n-1; i++)
	X                sum += Area2( P[0], P[i], P[i+1] );
	X        return  sum;
	X}
	X
	X/*
	X	Exclusive or: true iff exactly one argument is true.
	X	The arguments are negated to ensure that they are 0/1
	X	values.  Then the bitwise Xor operator may apply.
	X	(This idea is due to Michael Baldwin.)
	X*/
	Xbool 	Xor( bool x, bool y )
	X{
	X	return	!x ^ !y;
	X}
	X/*
	X	Returns true iff ab properly intersects cd: they share
	X	a point interior to both segments.  The properness of the
	X	intersection is ensured by using strict leftness.
	X*/
	Xbool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d )
	X{
	X	/* Eliminate improper cases. */
	X	if (
	X		Collinear(a,b,c) ||
	X		Collinear(a,b,d) ||
	X		Collinear(c,d,a) ||
	X		Collinear(c,d,b)
	X	   )
	X		return FALSE;
	X
	X        return
	X                Xor( Left(a,b,c), Left(a,b,d) )
	X        &&      Xor( Left(c,d,a), Left(c,d,b) );
	X}
	X
	X/*
	X	Returns true iff c is strictly to the left of the directed
	X	line through a to b.
	X*/
	Xbool    Left( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) > 0;
	X}
	X
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) >= 0;
	X}
	X
	Xbool    Collinear( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) == 0;
	X}
	X
	X/*
	X        Returns T iff (a,b,c) are collinear and point c lies 
	X	on the closed segement ab.
	X*/
	Xbool    Between( tPointi a, tPointi b, tPointi c )
	X{
	X        tPointi ba, ca;
	X        if ( ! Collinear( a, b, c ) )
	X                return  FALSE;
	X
	X	/* If ab not vertical, check betweenness on x; else on y. */
	X	if ( a[X] != b[X] ) 
	X		return	((a[X] <= c[X]) && (c[X] <= b[X])) ||
	X			((a[X] >= c[X]) && (c[X] >= b[X]));
	X	else
	X		return	((a[Y] <= c[Y]) && (c[Y] <= b[Y])) ||
	X			((a[Y] >= c[Y]) && (c[Y] >= b[Y]));
	X}
	X
	X/*
	X        Returns true iff segments ab and cd intersect, properly or improperly.
	X*/
	Xbool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d )
	X{
	X        if ( IntersectProp( a, b, c, d ) )
	X                return  TRUE;
	X        else if (       Between( a, b, c )
	X                     || Between( a, b, d )
	X                     || Between( c, d, a )
	X                     || Between( c, d, b )
	X                )
	X                return  TRUE;
	X        else    return  FALSE;
	X}
	X
	X/*
	X	Returns T iff (v_i, v_j) is a proper internal *or* external
	X	diagonal of P, *ignoring edges incident to v_i and v_j*.
	X*/
	Xbool	Diagonalie( int i, int j, int n, tPolygoni P )
	X{
	X	int	k;
	X	int	k1;
	X
	X	/* For each edge (k,k+1) of P */
	X	for ( k = 0; k < n; k++ ) {
	X		k1 = (k+1) % n;
	X		/* Skip edges incident to i or j */
	X		if ( ! (
	X				( k == i ) || ( k1 == i )
	X			||	( k == j ) || ( k1 == j )
	X		       )
	X		   )
	X			if ( Intersect( P[i], P[j], P[k], P[k1] ) )
	X				return	FALSE;
	X	}
	X	return	TRUE;
	X}
	X
	X/*
	X	Implements a = b, assignment of points/vectors.
	X	Assignment between arrays is not possible in C.
	X*/
	Xvoid	PointAssign( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < DIM; i++ )
	X		a[i] = b[i];
	X}
	X
	X/*
	X	Prints out n-3 diagonals (as pairs of integer indices)
	X	which form a triangulation of P.
	X*/
	Xvoid	Triangulate1( int n, tPolygoni P )
	X{
	X	tPolygoni	Pt;
	X	int	labels[PMAX];
	X	int	i;
	X
	X	for ( i = 0; i < n; i++ ){
	X		PointAssign( Pt[i], P[i] );
	X		labels[i] = i;
	X	}
	X
	X	TriRecurse( n, Pt, labels );
	X}
	X
	Xvoid	TriRecurse( int n, tPolygoni P, int labels[] )
	X{
	X	int	i, i1, i2;
	X
	X	printf("TriRecurse: n = %d\n", n);
	X	/*PrintPoly( n, P, labels );*/
	X	if ( n > 3 )
	X		for ( i = 0; i < n; i++ ) {
	X			i1 = (i+1) % n;
	X			i2 = (i+2) % n;
	X			if ( Diagonal( i, i2, n, P ) ) {
	X				printf("%d %d\n", labels[i], labels[i2]);
	X				ClipEar1( i1, n, P, labels );
	X				TriRecurse( n-1, P, labels );
	X				break;
	X			}
	X		}
	X}
	X
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X
	Xvoid	Triangulate( int n, tPolygoni P )
	X{
	X	int	i, i1, i2;
	X
	X	if ( n > 3 )
	X		for ( i = 0; i < n; i++ ) {
	X			i1 = (i+1) % n;
	X			i2 = (i+2) % n;
	X			if ( Diagonal( i, i2, n, P ) ) {
	X				PrintPoint( P[i]  ); putchar('\t');
	X				PrintPoint( P[i2] ); putchar('\n');
	X				ClipEar( i1, n, P );
	X				Triangulate( n-1, P );
	X				break;
	X			}
	X		}
	X}
	X
	X/*
	X	Returns true iff the diagonal (i,j) is striclty internal to the 
	X	polygon P in the neighborhood of the i endpoint.
	X*/
	Xbool	InCone( int i, int j, int n, tPolygoni P )
	X{
	X	int	i1;	/* i + 1 */
	X	int	in1;	/* i - 1 */
	X
	X	i1  = (i + 1)     % n;
	X	in1 = (i + n - 1) % n;
	X
	X	/* If P[i] is a convex vertex [ i+1 left or on (i-1,i) ]. */
	X	if ( LeftOn( P[in1], P[i], P[i1] ) )
	X		return	Left( P[i], P[j], P[in1] )
	X		     &&	Left( P[j], P[i], P[i1]  );
	X
	X	/* Assume (i-1,i,i+1) not collinear. */
	X	/* else P[i] is reflex. */
	X	else
	X		return	!(	LeftOn( P[i], P[j], P[i1]  )
	X		            &&	LeftOn( P[j], P[i], P[in1] ) );
	X}
	X
	X/*
	X	Returns T iff (v_i, v_j) is a proper internal
	X	diagonal of P.
	X*/
	Xbool	Diagonal( int i, int j, int n, tPolygoni P )
	X{
	X	return	InCone( i, j, n, P ) && Diagonalie( i, j, n, P );
	X}
	X
	X
	X/*
	X	Removes P[i] by copying P[i+1]...P[n-1] left one index.
	X*/
	Xvoid	ClipEar1( int i, int n, tPolygoni P, int labels[] )
	X{
	X	int	k;
	X	
	X	for ( k = i; k < n-1; k++ ) {
	X		PointAssign( P[k], P[k+1] );
	X		labels[k] = labels[k+1];
	X	}
	X}
	X
	Xvoid	ClipEar( int i, int n, tPolygoni P )
	X{
	X	int	k;
	X	
	X	for ( k = i; k < n-1; k++ )
	X		PointAssign( P[k], P[k+1] );
	X}
	X
	Xint	ReadPoints( tPolygoni P )
	X/*
	X        Reads in the coordinates of the vertices of a polygon from stdin,
	X        puts them into P, and returns n, the number of vertices.
	X        The input is assumed to be pairs of whitespace-separated coordinates,
	X        one pair per line.  The number of points is not part of the input.
	X*/
	X{
	X	int	n = 0;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < PMAX) && 
	X		(scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in ReadPoints:\
	X			too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P, int labels[] )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   l   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d%4d\n", 
	X			i, labels[i], P[i][0], P[i][1]);
	X}
	X/*
	X        Puts a - b into c.
	X*/
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X        int     i;
	X        
	X        for( i = 0; i < DIM; i++ )
	X                c[i] = a[i] - b[i];
	X}
	X
	X/*
	X        Returns the dot product of the two input vectors.
	X*/
	Xint     Dot( tPointi a, tPointi b )
	X{
	X        int     i;
	X        int     sum = 0;
	X
	X        for( i = 0; i < DIM; i++ )
	X                sum += a[i] * b[i];
	X        
	X        return  sum;
	X}
	X
	X/*
	X        Returns the square of the length of the vector p.
	X*/
	Xint     Length2( tPointi p )
	X{
	X        return  Dot( p, p );
	X}
SHAR_EOF
if test 9289 -ne "`wc -c < 'tri.c'`"
then
	echo shar: error transmitting "'tri.c'" '(should have been 9289 characters)'
fi
fi # end of overwriting check
echo shar: extracting "'vector.c'" '(2570 characters)'
if test -f 'vector.c'
then
	echo shar: will not over-write existing file "'vector.c'"
else
sed 's/^	X//' << \SHAR_EOF > 'vector.c'
	X/* 
	X        Returns twice the signed area of the triangle determined by a,b,c.
	X        The area is positive if a,b,c are oriented ccw, negative if cw,
	X	and zero if the points are collinear.
	X*/
	Xint     Area2( tPointi a, tPointi b, tPointi c )
	X{
	X        return
	X                a[0] * b[1] - a[1] * b[0] +
	X                a[1] * c[0] - a[0] * c[1] +
	X                b[0] * c[1] - c[0] * b[1];
	X}
	X
	X/*
	X	Exclusive or: true iff exactly one argument is true.
	X	The arguments are negated to ensure that they are 0/1
	X	values.  Then the bitwise xor operator may apply.
	X	(This idea is due to Michael Baldwin.)
	X*/
	Xbool 	xor( bool x, bool y )
	X{
	X	return	!x ^ !y;
	X}
	Xbool    Left( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) > 0;
	X}
	X
	Xbool    LeftOn( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) >= 0;
	X}
	X
	Xbool    Collinear( tPointi a, tPointi b, tPointi c )
	X{
	X        return  Area2( a, b, c ) == 0;
	X}
	X
	X/*
	X        Puts a - b into c.
	X*/
	Xvoid	SubVec( tPointi a, tPointi b, tPointi c )
	X{
	X        int     i;
	X        
	X        for( i = 0; i < DIM; i++ )
	X                c[i] = a[i] - b[i];
	X}
	X
	X/*
	X        Returns the dot product of the two input vectors.
	X*/
	Xint     Dot( tPointi a, tPointi b )
	X{
	X        int     i;
	X        int     sum = 0;
	X
	X        for( i = 0; i < DIM; i++ )
	X                sum += sum + a[i] * b[i];
	X        
	X        return  sum;
	X}
	X
	X/*
	X        Returns the square of the length of the vector p.
	X*/
	Xint     Length2( tPointi p )
	X{
	X        return  Dot( p, p );
	X}
	X
	X/*
	X	Implements a = b, assignment of points/vectors.
	X	Assignment between arrays is not possible in C.
	X*/
	Xvoid	PointAssign( tPointi a, tPointi b )
	X{
	X	int	i;
	X
	X	for ( i = 0; i < DIM; i++ )
	X		a[i] = b[i];
	X}
	Xvoid	PrintPoint( tPointi p )
	X{
	X	int	i;
	X
	X	putchar('(');
	X	for ( i = 0; i < DIM; i++ ) {
	X		printf("%d", p[i]);
	X		if ( i != DIM-1 ) putchar(',');
	X	}
	X	putchar(')');
	X}
	X
	X/*
	X	Reads in the coordinates of the vertices of a polygon from stdin,
	X	puts them into P, and returns n, the number of vertices.
	X	Formatting conventions: etc.
	X*/
	Xint	ReadPoly( tPolygoni P )
	X{
	X	int	n = 0;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
	X		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
	X		++n;
	X	}
	X	if (n < PMAX)
	X		printf("n = %3d vertices read\n",n);
	X	else	printf("Error in read_poly:  too many points; max is %d\n", PMAX);
	X	putchar('\n');
	X
	X	return	n;
	X}
	X
	Xvoid	PrintPoly( int n, tPolygoni P )
	X{
	X	int	i;
	X
	X	printf("Polygon:\n");
	X	printf("  i   x   y\n");
	X	for( i = 0; i < n; i++ )
	X		printf("%3d%4d%4d\n", i, P[i][0], P[i][1]);
	X}
SHAR_EOF
if test 2570 -ne "`wc -c < 'vector.c'`"
then
	echo shar: error transmitting "'vector.c'" '(should have been 2570 characters)'
fi
fi # end of overwriting check
#	End of shell archive
exit 0