Method IO Graphics PART04.SH

# ! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of archive 4 (of 5)."
# Contents:  2DClip/cross.c AALines/AALines.c AALines/AATables.c
#   AAPolyScan.c ConcaveScan.c Forms.c Interleave.c Sturm/sturm.c
# Wrapped by craig@weedeater on Wed Dec 12 20:45:15 1990
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f '2DClip/cross.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'2DClip/cross.c'\"
else
echo shar: Extracting \"'2DClip/cross.c'\" \(6788 characters\)
sed "s/^X//" >'2DClip/cross.c' <<'END_OF_FILE'
X
X/*
X * file cross.c:
X *	calculate the intersections
X */
X#include	<math.h>
X#include	"GraphicsGems.h"
X#include	"line.h"
X#include	"box.h"
X/*
X * cross_calc:
X *
X *	PURPOSE
X *	calculate the intersections between the polygon
X *	stored in poly and the linesegment stored in l
X *	and put the intersections into psol.
X *
X *	poly	pointer to the structure containing the polygon
X *	l	pointer to the structure containing the linesegment
X *	psol	pointer to the pointer where intersections are stored
X *	nsol	current number of intersections stored
X *	nsmax	maximum storage in psol for intersections
X *		if nsol exceeds nsmax additional storage is allocated
X *
X */
Xcross_calc(poly, l, psol, nsol, nsmax)
XCONTOUR	*poly;
XSEGMENT	*l;
XCLIST	**psol;
Xshort	*nsol, nsmax;
X{
X	SEGMENT	*p;
X	CLIST	*sol;
X	double	s;
X	long	x, y, a, b, c;
X	int	psort(), type;
X
X	sol = *psol;
X	p = poly->_s;
X	do {
X/*
X * calculate the a, b and c coefficients and determine the
X * type of intersection
X */
X
X		a = (p->_to._y - p->_from._y)*(l->_to._x - l->_from._x) -
X		    (p->_to._x - p->_from._x)*(l->_to._y - l->_from._y);
X		b = (p->_from._x - l->_from._x)*(l->_to._y - l->_from._y) -
X		    (p->_from._y - l->_from._y)*(l->_to._x - l->_from._x);
X		c = (p->_from._x - l->_from._x)*(p->_to._y - p->_from._y) -
X		    (p->_from._y - l->_from._y)*(p->_to._x - p->_from._x);
X		if(a == 0)
X			type = (b == 0) ? COINCIDE : PARALLEL;
X		else {
X			if(a > 0) {
X				if((b >= 0 && b <= a) &&
X				   (c >= 0 && c <= a))
X					type = CROSS;
X				else
X					type = NO_CROSS;
X			}
X			else {
X				if((b <= 0 && b >= a) &&
X				   (c <= 0 && c >= a))
X					type = CROSS;
X				else
X					type = NO_CROSS;
X			}
X		}
X/*
X * process the interscetion found
X */
X		switch(type) {
X			case NO_CROSS: case PARALLEL:
X				break;
X
X			case CROSS:
X				if(b == a || c == a || c == 0)
X					break;
X				if(b == 0 && 
X				   p_where(&(p->_prev->_from), &(p->_to), l) >= 0)
X					break;
X				s = (double)b/(double)a;
X				if(l->_from._x == l->_to._x)
X					x = l->_from._x;
X				else
X					x = p->_from._x + 
X					   (int)((p->_to._x - p->_from._x)*s);
X				if(l->_from._y == l->_to._y)
X					y = l->_from._y;
X				else
X					y = p->_from._y + 
X					   (int)((p->_to._y - p->_from._y)*s);
X
X				if(*nsol == nsmax) {
X					nsmax *= 2;
X					*psol = sol = (CLIST *) realloc(sol, 							nsmax*sizeof(CLIST)); 
X				}
X				sol[*nsol]._p._x = x;
X				sol[*nsol]._p._y = y;
X				sol[*nsol]._type = STD;
X				*nsol += 1;
X				break;
X
X			case COINCIDE:
X				if(p_where(&(p->_prev->_from),
X					 &(p->_next->_to), l) > 0)
X					break;
X				if(l->_from._x != l->_to._x) {
X					if((max(l->_from._x, l->_to._x) <
X					    min(p->_from._x, p->_to._x)  ) ||
X					   (min(l->_from._x, l->_to._x) >
X					    max(p->_from._x, p->_to._x))   )
X						break;
X					if(min(l->_from._x, l->_to._x) <
X					   min(p->_from._x, p->_to._x) ) {
X						if(*nsol == nsmax) {
X							nsmax *= 2;
X							*psol = sol = (CLIST *) realloc(sol, 								nsmax*sizeof(CLIST));
X						}
X						sol[*nsol]._p._x = p->_from._x;
X						sol[*nsol]._p._y = p->_from._y;
X						sol[*nsol]._type = DELAY;
X						*nsol += 1;
X					}
X					if(max(l->_from._x, l->_to._x) >
X					   max(p->_from._x, p->_to._x) ) {
X						if(*nsol == nsmax) {
X							nsmax *= 2;
X							*psol = sol = (CLIST *) realloc(sol, 								nsmax*sizeof(CLIST));
X						}
X						sol[*nsol]._p._x = p->_to._x;
X						sol[*nsol]._p._y = p->_to._y;
X						sol[*nsol]._type = DELAY;
X						*nsol += 1;
X					}
X				}
X				else {
X
X					if((max(l->_from._y, l->_to._y) <
X					    min(p->_from._y, p->_to._y)  ) ||
X					   (min(l->_from._y, l->_to._y) >
X					    max(p->_from._y, p->_to._y)) )
X						break;
X					if(min(l->_from._y, l->_to._y) <
X					   min(p->_from._y, p->_to._y) ) {
X						if(*nsol == nsmax) {
X							nsmax *= 2;
X							*psol = sol = (CLIST *) realloc(sol, 								nsmax*sizeof(CLIST));
X						}
X						sol[*nsol]._p._x = p->_from._x;
X						sol[*nsol]._p._y = p->_from._y;
X						sol[*nsol]._type = DELAY;
X						*nsol += 1;
X					}
X					if(max(l->_from._y, l->_to._y) >
X					   max(p->_from._y, p->_to._y) ) {
X						if(*nsol == nsmax) {
X							nsmax *= 2;
X							*psol = sol = (CLIST *) realloc(sol, 								nsmax*sizeof(CLIST));
X						}
X						sol[*nsol]._p._x = p->_to._x;
X						sol[*nsol]._p._y = p->_to._y;
X						sol[*nsol]._type = DELAY;
X						*nsol += 1;
X					}
X				}
X				break;
X		}
X		p = p->_next;
X	} while(p != poly->_s);
X	qsort(sol, *nsol, sizeof(CLIST), psort);
X}
X
X
X/*
X * p_where
X *
X *	PURPOSE
X *	determine position of point p1 and p2 relative to
X *	linesegment l. 
X *	Return value
X *	< 0	p1 and p2 lie at different sides of l
X *	= 0	one of both points lie on l
X *	> 0	p1 and p2 lie at same side of l
X *
X *	p1	pointer to coordinates of point
X *	p2	pointer to coordinates of point
X *	l	pointer to linesegment
X * 
X */
Xp_where(p1, p2, l)
XPOINT	*p1, *p2;
XSEGMENT	*l;
X{
X	long	dx, dy, dx1, dx2, dy1, dy2, p_1, p_2;
X
X	dx  = l->_to._x - l->_from._x;
X	dy  = l->_to._y - l->_from._y;
X	dx1 = p1->_x - l->_from._x;
X	dy1 = p1->_y - l->_from._y;
X	dx2 = p2->_x - l->_to._x;
X	dy2 = p2->_y - l->_to._y;
X	p_1 = (dx*dy1 - dy*dx1);
X	p_2 = (dx*dy2 - dy*dx2);
X	if(p_1 == 0 || p_2 == 0)
X		return(0);
X	else {
X		if((p_1 > 0 && p_2 < 0) || (p_1 < 0 && p_2 > 0))
X			return(-1);
X		else
X			return(1);
X	}
X}
X
X
X/*
X * p_inside
X *
X *	PURPOSE
X *	determine if the point stored in pt lies inside
X *	the polygon stored in p
X *	Return value:
X *	FALSE	pt lies outside p
X *	TRUE	pt lies inside  p
X *
X *	p	pointer to the polygon
X *	pt	pointer to the point
X */	
Xboolean	p_inside(p, pt)
XCONTOUR	*p;
XPOINT	*pt;
X{
X	SEGMENT	l, *sp;
X	CLIST	*sol;
X	short	nsol = 0, nsmax = 2, reduce = 0, i;
X	boolean	on_contour(), odd;
X	
X	l._from._x = p->_minx-2;
X	l._from._y = pt->_y;
X	l._to._x   = pt->_x;
X	l._to._y   = pt->_y;
X	sol = (CLIST *) calloc(2, sizeof(CLIST));
X	cross_calc(p, &l, &sol, &nsol, nsmax);
X	for(i=0; i<nsol-1; i++)
X		if(sol[i]._type == DELAY && sol[i+1]._type == DELAY)
X			reduce++;
X	free(sol);
X	odd = (nsol - reduce) & 0x01;
X	return(odd ? !on_contour(p, pt) : FALSE);
X}
X
X/*
X * function used for sorting
X */
Xpsort(p1, p2)
XCLIST	*p1, *p2;
X{
X	if(p1->_p._x != p2->_p._x)
X		return(p1->_p._x - p2->_p._x);
X	else 
X		return(p1->_p._y - p2->_p._y);
X}
X
X
X/*
X * on_contour
X *
X *	PURPOSE
X *	determine if the point pt lies on the
X *	contour p.
X *	Return value
X *	TRUE 	point lies on contour
X *	FALSE	point lies not on contour
X *
X *	p	pointer to the polygon structure
X *	pt	pointer to the point
X */
Xboolean	on_contour(p, pt)
XCONTOUR	*p;
XPOINT	*pt;
X{
X	SEGMENT	*sp;
X	long	dx1, dy1, dx2, dy2;
X
X	sp = p->_s;
X	do {
X		if((pt->_x >= min(sp->_from._x, sp->_to._x)) &&
X		   (pt->_x <= max(sp->_from._x, sp->_to._x))    ) { 
X			dx1 = pt->_x - sp->_from._x;
X			dx2 = sp->_to._x - pt->_x;
X			dy1 = pt->_y - sp->_from._y;
X			dy2 = sp->_to._y - pt->_y;
X			if(dy1*dx2 == dy2*dx1)
X				return(TRUE);
X		}
X		sp = sp->_next;
X	} while(sp != p->_s);
X	return(FALSE);
X}
X 
END_OF_FILE
if test 6788 -ne `wc -c <'2DClip/cross.c'`; then
    echo shar: \"'2DClip/cross.c'\" unpacked with wrong size!
fi
# end of '2DClip/cross.c'
fi
if test -f 'AALines/AALines.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'AALines/AALines.c'\"
else
echo shar: Extracting \"'AALines/AALines.c'\" \(5139 characters\)
sed "s/^X//" >'AALines/AALines.c' <<'END_OF_FILE'
X/*  FILENAME: AALines.c  [revised 17 AUG 90]
X
X    AUTHOR:  Kelvin Thompson
X
X    DESCRIPTION:  Code to render an anti-aliased line, from
X      "Rendering Anti-Aliased Lines" in _Graphics_Gems_.
X
X      This is almost exactly the code printed on pages 690-693
X      of _Graphics_Gems_.  An overview of the code is on pages
X      105-106.
X
X    LINK WITH:
X      AALines.h -- Shared tables, symbols, etc.
X      AAMain.c -- Calling code for this subroutine.
X      AATables.c -- Initialize lookup tables.
X*/
X
X
X#include "AALines.h"
X
X#define SWAPVARS(v1,v2) ( temp=v1, v1=v2, v2=temp )
X
X#define FIXMUL(f1,f2) \
X  ( \
X  (((f1&0x0000ffff) * (f2&0x0000ffff)) >> 16) + \
X  ((f1&0xffff0000)>>16) * (f2&0x0000ffff) + \
X  ((f2&0xffff0000)>>16) * (f1&0x0000ffff) + \
X  ((f1&0xffff0000)>>16) * (f2&0xffff0000) \
X  )
X
X/* HARDWARE ASSUMPTIONS:
X/*    * 32-bit, signed ints
X/*    * 8-bit pixels, with initialized color table
X/*    * pixels are memory mapped in a rectangular fashion */
X
X/* FIXED-POINT DATA TYPE */
X#ifndef FX_FRACBITS
X  typedef int FX;
X# define FX_FRACBITS 16  /* bits of fraction in FX format */
X# define FX_0        0   /* zero in fixed-point format */
X#endif
X
X/* ASSUMED MACROS:
X/*   SWAPVARS(v1,v2) -- swaps the contents of two variables
X/*   PIXADDR(x,y) -- returns address of pixel at (x,y)
X/*   COVERAGE(FXdist) -- lookup macro for pixel coverage 
X        given perpendicular distance; takes a fixed-point
X        integer and returns an integer in the range [0,255]
X/*   SQRTFUNC(FXval) -- lookup macro for sqrt(1/(1+FXval^2))
X        accepts and returns fixed-point numbers
X/*   FIXMUL(FX1,FX2) -- multiplies two fixed-point numbers
X        and returns the product as a fixed-point number   */
X
X/* BLENDING FUNCTION:
X/*  'cover' is coverage -- in the range [0,255]
X/*  'back' is background color -- in the range [0,255] */
X#define BLEND(cover,back) ((((255-(cover))*(back))>>8)+(cover))
X
X/* LINE DIRECTION bits and tables */
X#define DIR_STEEP  1  /* set when abs(dy) > abs(dx) */
X#define DIR_NEGY   2  /* set whey dy < 0 */
X
X
X/* pixel increment values 
X/*   -- assume PIXINC(dx,dy) is a macro such that:
X/*   PIXADDR(x0,y0) + PIXINC(dx,dy) = PIXADDR(x0+dx,y0+dy)  */
Xstatic int adj_pixinc[4] = 
X      { PIXINC(1,0), PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1) };
Xstatic int diag_pixinc[4] = 
X      { PIXINC(1,1), PIXINC(1,1), PIXINC(1,-1), PIXINC(1,-1) };
Xstatic int orth_pixinc[4] = 
X      { PIXINC(0,1), PIXINC(1,0), PIXINC(0,-1), PIXINC(1,0) };
X
X/* Global 'Pmax' is initialized elsewhere.  It is the
X   "maximum perpendicular distance" -- the sum of half the
X   line width and the effective pixel radius -- in fixed format */
XFX Pmax;
X
X
X/***************  FUNCTION ANTI_LINE  ***************/
X
Xvoid Anti_Line ( X1, Y1, X2, Y2 )
Xint X1, Y1, X2, Y2;
X{
Xint 	Bvar, 	/* decision variable for Bresenham's */
X    	Bainc,   /* adjacent-increment for 'Bvar' */
X    	Bdinc;   /* diagonal-increment for 'Bvar' */
XFX 		Pmid,  	/* perp distance at Bresenham's pixel */
X   		Pnow,  	/* perp distance at current pixel (ortho loop) */
X   		Painc, 	/* adjacent-increment for 'Pmid' */
X   		Pdinc, 	/* diagonal-increment for 'Pmid' */
X   		Poinc; 	/* orthogonal-increment for 'Pnow'--also equals 'k' */
Xchar 	*mid_addr,   /* pixel address for Bresenham's pixel */
X     	*now_addr;   /* pixel address for current pixel */
Xint 	addr_ainc,   /* adjacent pixel address offset */
X    	addr_dinc,   /* diagonal pixel address offset */
X    	addr_oinc;   /* orthogonal pixel address offset */
Xint dx,dy,dir;    	/* direction and deltas */
XFX slope;			/* slope of line */
Xint temp;
X
X/* rearrange ordering to force left-to-right */
Xif 	( X1 > X2 )
X  	{ SWAPVARS(X1,X2);  SWAPVARS(Y1,Y2); }
X
X/* init deltas */
Xdx = X2 - X1;  /* guaranteed non-negative */
Xdy = Y2 - Y1;
X
X
X/* calculate direction (slope category) */
Xdir = 0;
Xif ( dy < 0 )   { dir |= DIR_NEGY;  dy = -dy; }
Xif ( dy > dx )  { dir |= DIR_STEEP; SWAPVARS(dx,dy); }
X
X/* init address stuff */
Xmid_addr = PIXADDR(X1,Y1);
Xaddr_ainc = adj_pixinc[dir];
Xaddr_dinc = diag_pixinc[dir];
Xaddr_oinc = orth_pixinc[dir];
X
X/* perpendicular measures */
Xslope =  (dy << FX_FRACBITS) / dx;
XPoinc = SQRTFUNC( slope );
XPainc = FIXMUL( slope, Poinc );
XPdinc = Painc - Poinc;
XPmid = FX_0;
X
X/* init Bresenham's */
XBainc = dy << 1;
XBdinc = (dy-dx) << 1;
XBvar = Bainc - dx;
X
Xdo
X  	{
X  	/* do middle pixel */
X  	*mid_addr = BLEND( COVERAGE(abs(Pmid)), *mid_addr );
X
X  	/* go up orthogonally */
X  	for (
X      	Pnow = Poinc-Pmid,  now_addr = mid_addr+addr_oinc;
X      	Pnow < Pmax;
X      	Pnow += Poinc,      now_addr += addr_oinc
X      	)
X       *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
X
X  	/* go down orthogonally */
X  	for (
X      	Pnow = Poinc+Pmid,  now_addr = mid_addr-addr_oinc;
X      	Pnow < Pmax;
X      	Pnow += Poinc,      now_addr -= addr_oinc
X      	)
X       *now_addr = BLEND( COVERAGE(Pnow), *now_addr );
X
X
X  	/* update Bresenham's */
X  	if ( Bvar < 0 )
X    	{
X    	Bvar += Bainc;
X    	mid_addr = (char *) ((int)mid_addr + addr_ainc);
X    	Pmid += Painc;
X    	}
X  	else
X    	{
X    	Bvar += Bdinc;
X    	mid_addr = (char *) ((int)mid_addr + addr_dinc);
X    	Pmid += Pdinc;
X    	}
X
X  	--dx;
X  	} while ( dx >= 0 );
X}
END_OF_FILE
if test 5139 -ne `wc -c <'AALines/AALines.c'`; then
    echo shar: \"'AALines/AALines.c'\" unpacked with wrong size!
fi
# end of 'AALines/AALines.c'
fi
if test -f 'AALines/AATables.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'AALines/AATables.c'\"
else
echo shar: Extracting \"'AALines/AATables.c'\" \(5772 characters\)
sed "s/^X//" >'AALines/AATables.c' <<'END_OF_FILE'
X/*  FILENAME: AATables.c  [revised 18 AUG 90]
X
X    DESCRIPTION:  Initialization of lookup tables and frame buffer
X      for anti-aliased line rendering demo.
X
X    LINK WITH:
X      AALines.h -- Shared variables and symbols for renderer.
X      AAMain.c -- Calling routine for renderer.
X      AALines.c -- Anti-aliased line rendering code.
X*/
X
X#include <math.h>
X#include "AALines.h"
X
X/* programs in this file */
Xextern void Anti_Init();
Xstatic void Sqrt_Init();
X
X/* globals defined here */
Xchar *fbuff;
XUFX *sqrtfunc=0;
Xint sqrtcells=1024;
Xint sqrtshift;
X
X/* AA sizes */
Xfloat line_r=1.0;     /* line radius */
Xfloat pix_r=SQRT_2;   /* pixel radius */
XFX *coverage=0;
Xint covercells=128;
Xint covershift;
X
X
X
X
X/*  ************************************
X    *                                  *
X    *          Anti_Init               *
X    *                                  *
X    ************************************
X
X    DESCRIPTION:  Initialize everything for the anti-aliased line
X      renderer in 'AALines.c':  allocate the frame buffer,
X      set up lookup tables, etc.
X
X      For hints about initializing the coverage table, see 
X      "Area of Intersection: Circle and a Thick Line" on pages
X      40-42 of _Graphics_Gems_.
X*/
X
X#define FLOAT_TO_CELL(flt)  ((int) ((flt) * 255.0 + 0.5))
X#define MAXVAL_CELL         255
X
Xvoid Anti_Init ( )
X{
Xint *thiscell;
Xdouble maxdist,nowdist,incdist;
Xint tablebits,radbits;
Xint tablecells;
Xstatic int tablesize=0;
Xdouble fnear,ffar,fcover;
Xdouble half,invR,invpiRsq,invpi,Rsq;
Xdouble sum_r;
Xdouble inv_log_2;
Xextern FX Pmax;
X
X/* alloc & init frame buffer */
Xfbuff = (char *) malloc( xpix*ypix );
X  {
X  register int i;
X  for ( i=xpix*ypix-1; i>=0; --i )  fbuff[i] = 0;
X  }
X
X/* init */
Xinv_log_2 = 1.0 / log( 2.0 );
Xsum_r = line_r + pix_r;
Xtablebits = (int) ( log((double)covercells) * inv_log_2 + 0.99 );
Xradbits = (int) ( log((double)sum_r) * inv_log_2 ) + 1;
Xcovershift = FX_FRACBITS - (tablebits-radbits);
X
X/* constants */
Xhalf = 0.5;
XinvR = 1.0 / pix_r;
Xinvpi = 1.0 / PI;
XinvpiRsq = invpi * invR * invR;
XRsq = pix_r * pix_r;
X#define FRACCOVER(d) (half - d*sqrt(Rsq-d*d)*invpiRsq - invpi*asin(d*invR))
X
X/* allocate table */
XPmax = FLOAT_TO_FX(sum_r);
XPmax >>= covershift;
Xtablecells = Pmax + 2;
XPmax <<= covershift;
Xif ( coverage  &&  tablecells > tablesize )
X  { free( coverage ); coverage = 0;  tablesize = 0; }
Xif ( coverage == 0 )
X  {
X  coverage = (FX *) malloc( tablecells * sizeof(int) );
X  tablesize = tablecells;
X  }
X
X/* init for fill loops */
Xnowdist = 0.0;
Xthiscell = coverage;
Xincdist = sum_r / (double)(tablecells-2);
X
X/* fill fat portion */
Xif ( pix_r <= line_r )
X  {
X  maxdist = line_r - pix_r;
X  for ( ;
X      nowdist <= maxdist;
X      nowdist += incdist,  ++thiscell
X      )
X    {
X    *thiscell = MAXVAL_CELL;
X    }
X  }
X
X/* fill skinny portion */
Xelse
X  {
X  /* loop till edge of line, or end of skinny, whichever comes first */
X  maxdist = pix_r - line_r;
X  if ( maxdist > line_r )
X    maxdist = line_r;
X  for ( ;
X      nowdist < maxdist;
X      nowdist += incdist,  ++thiscell
X      )
X    {
X    fnear = line_r - nowdist;
X    ffar = line_r + nowdist;
X    fcover = 1.0 - FRACCOVER(fnear) - FRACCOVER(ffar);
X    *thiscell = FLOAT_TO_CELL(fcover);
X    }
X
X  /* loop till end of skinny -- only run on super-skinny */
X  maxdist = pix_r - line_r;
X  for ( ;
X      nowdist < maxdist;
X      nowdist += incdist,  ++thiscell
X      )
X    {
X    fnear = nowdist - line_r;
X    ffar = nowdist + line_r;
X    fcover = FRACCOVER(fnear) - FRACCOVER(ffar);
X    *thiscell = FLOAT_TO_CELL(fcover);
X    }
X  }
X
X/* loop till edge of line */
Xmaxdist = line_r;
Xfor ( ;
X    nowdist < maxdist;
X    nowdist += incdist,  ++thiscell
X    )
X  {
X  fnear = line_r - nowdist;
X  fcover = 1.0 - FRACCOVER(fnear);
X  *thiscell = FLOAT_TO_CELL(fcover);
X  }
X
X/* loop till max separation */
Xmaxdist = line_r + pix_r;
Xfor ( ;
X    nowdist < maxdist;
X    nowdist += incdist,  ++thiscell
X    )
X  {
X  fnear = nowdist - line_r;
X  fcover = FRACCOVER(fnear);
X  *thiscell = FLOAT_TO_CELL(fcover);
X  }
X
X/* finish off table */
X*thiscell = FLOAT_TO_CELL(0.0);
Xcoverage[tablecells-1] = FLOAT_TO_CELL(0.0);
X
XSqrt_Init();
X}
X
X
X
X
X/*  *******************************
X    *                             *
X    *       Sqrt_Init             *
X    *                             *
X    *******************************
X
X    DESCRIPTION:  Initialize the lookup table for the function
X      sqrt(1/(1+x^2))).  The table takes a shifted fixed-point
X      value as an index and returns a fixed-point value.  Input
X      values are in the range [0,1] inclusive.  The number of
X      cells in the table is a power of two plus one (the extra
X      cell provides an entry for an input value of exactly 1).
X
X    GLOBALS:
X      sqrtcells -- Number of cells to use in the table (must
X	be set before calling this routine).  This number is
X	rounded up to the nearest power of two (the global
X	variable itself is unchanged).
X      sqrtshift -- Bits to shift a fixed-point (FX) number
X	to generate a table index.
X      sqrtfunc -- Lookup table for the function.
X*/
X
Xstatic void Sqrt_Init ( )
X{
XUFX *thiscell;
Xdouble nowval,incval;
Xint tablebits;
Xint tablecells;
Xdouble one;
X
X/* init */
Xtablebits = (int) ( log((double)sqrtcells) / log(2.0) + 0.999 );
Xtablecells = (1 << tablebits) + 1;  /* one more that requested */
Xsqrtshift = FX_FRACBITS - tablebits;
Xone = 1.0;
X
X/* allocate table */
Xif ( sqrtfunc )
X  free( sqrtfunc );
Xsqrtfunc = (UFX *) malloc( tablecells * sizeof(int) );
X
X/* init for fill loop */
Xincval = one / (double)(tablecells-1);  /* a negative power of two */
X
Xfor (
X    nowval = 0.0,      thiscell = sqrtfunc;
X    nowval < 1.0;
X    nowval += incval,  ++thiscell
X    )
X  {
X  *thiscell = FLOAT_TO_FX( sqrt(one/(one+nowval*nowval)) );
X  }
X
Xsqrtfunc[tablecells-1] = FLOAT_TO_FX( sqrt(0.5) );
X}
END_OF_FILE
if test 5772 -ne `wc -c <'AALines/AATables.c'`; then
    echo shar: \"'AALines/AATables.c'\" unpacked with wrong size!
fi
# end of 'AALines/AATables.c'
fi
if test -f 'AAPolyScan.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'AAPolyScan.c'\"
else
echo shar: Extracting \"'AAPolyScan.c'\" \(7519 characters\)
sed "s/^X//" >'AAPolyScan.c' <<'END_OF_FILE'
X/* 
XFast Anti-Aliasing Polygon Scan Converstion
Xby Jack Morrison
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*
X* Anti-aliased polygon scan conversion by Jack Morrison
X*
X* This code renders a polygon, computing subpixel coverage at
X* 8 times Y and 16 times X display resolution for anti-aliasing.
X* One optimization left out for clarity is the use of incremental
X* interpolations. X coordinate interpolation in particular can be
X* with integers. See Dan Field's article in ACM Transactions on
X* Graphics, January 1985 for a fast incremental interpolator.
X*/
X#include "GraphicsGems.h"
X
X#define	SUBYRES	8		/* subpixel Y resolution per scanline */
X#define	SUBXRES	16		/* subpixel X resolution per pixel */
X#define	MAX_AREA	(SUBYRES*SUBXRES)
X#define	MODRES(y)	((y) & 7)		/*subpixel Y modulo */
X#define MAX_X	0x7FFF	/* subpixel X beyond right edge */
X
Xtypedef struct SurfaceStruct {  /* object shading surface info */
X	int	red, green, blue;		   /* color components */
X	} Surface;
X/*
X* In  real life, SurfaceStruct will contain many more parameters as
X* required by the shading and rendering programs, such as diffuse
X* and specular factors, texture information, transparency, etc.
X*/
X
Xtypedef struct VertexStruct	{	/* polygon vertex */
X	Vector3	model, world,		/* geometric information */
X		    normal, image;
X	int y;					/* subpixel display coordinate */
X	} Vertex;
X
XVertex *Vleft, *VnextLeft;		/* current left edge */
XVertex *Vright, *VnextRight;	/* current right edge */
X
Xstruct	SubPixel  {			/* subpixel extents for scanline */
X	int xLeft, xRight;
X	} sp[SUBYRES];
X
Xint	xLmin, xLmax;		/* subpixel x extremes for scanline */
Xint	xRmax, xRmin;		/* (for optimization shortcut) */
X
X/* Compute sub-pixel x coordinate for vertex */
Xextern int screenX(/* Vertex *v */);
X
X/* Interpolate vertex information */
Xextern void vLerp(/* double alpha, Vertex *Va, *Vb, *Vout */);
X
X/* Render polygon for one pixel, given coverage area */
X/*  and bitmask */
Xextern void renderPixel(/* int x, y, Vertex *V,
X						int area, unsigned mask[], 
X						Surface *object */);
X
X/*
X * Render shaded polygon
X */
XdrawPolygon(polygon, numVertex, object)
X	Vertex	polygon[];		/*clockwise clipped vertex list */
X	int	numVertex;			/*number of vertices in polygon */
X
X	Surface *object;			/* shading parms for this object */
X{
X	Vertex *endPoly;			/* end of polygon vertex list */
X	Vertex VscanLeft, VscanRight;	/* interpolated vertices */ 								/* at scanline */
X	double aLeft, aRight;			/* interpolation ratios */
X	struct SubPixel *sp_ptr;		/* current subpixel info */
X	int xLeft, xNextLeft;			/* subpixel coordinates for */
X	int  xRight, xNextRight;		/* active polygon edges */
X	int i,y;						
X
X/* find vertex with minimum y (display coordinate) */
XVleft = polygon;
Xfor  (i=1; i<numVertex; i++)
X	if  (polygon[i].y < Vleft->y)
X		Vleft = &polygon[i];
XendPoly = &polygon[numVertex-1];
X
X/* initialize scanning edges */
XVright = VnextRight = VnextLeft = Vleft;
X
X/* prepare bottom of initial scanline - no coverage by polygon */
Xfor (i=0; i<SUBYRES; i++)
X	sp[i].xLeft = sp[i].xRight = -1;
XxLmin = xRmin = MAX_X;
XxLmax = xRmax = -1;
X
X/* scan convert for each subpixel from bottom to top */
Xfor (y=Vleft->y; ; y++)	{
X
X	while (y == VnextLeft->y)	{	/* reached next left vertex */
X		VnextLeft = (Vleft=VnextLeft) + 1; 	/* advance */
X		if (VnextLeft > endPoly)			/* (wraparound) */
X			VnextLeft = polygon;
X		if (VnextLeft == Vright)	/* all y's same?  */
X			return;				/* (null polygon) */ 
X		xLeft = screenX(Vleft);
X		xNextLeft = screenX(VnextLeft);
X	}
X
X	while (y == VnextRight->y)  { /*reached next right vertex */
X		VnextRight = (Vright=VnextRight) -1;
X		if (VnextRight < polygon)			/* (wraparound) */
X			VnextRight = endPoly;
X		xRight = screenX(Vright);
X		xNextRight = screenX(VnextRight);
X	}
X
X	if (y>VnextLeft->y || y>VnextRight->y)	{
X				/* done, mark uncovered part of last scanline */
X		for (; MODRES(y); y++)
X			sp[MODRES(y)].xLeft = sp[MODRES(y)].xRight = -1;
X		renderScanline(Vleft, Vright, y/SUBYRES, object);
X		return;
X	}
X
X/*
X * Interpolate sub-pixel x endpoints at this y,
X * and update extremes for pixel coherence optimization
X */
X	
X	sp_ptr = &sp[MODRES(y)];
X	aLeft = (double)(y - Vleft->y) / (VnextLeft->y - Vleft->y);
X	sp_ptr->xLeft = LERP(aLeft, xLeft, xNextLeft);
X	if (sp_ptr->xLeft < xLmin)
X		xLmin = sp_ptr->xLeft;
X	if (sp_ptr->xLeft > xLmax)
X		xLmax = sp_ptr->xLeft;
X
X	aRight = (double)(y - Vright->y) / (VnextRight->y 
X					- Vright->y);
X	sp_ptr->xRight = LERP(aRight, xRight, xNextRight);
X	if (sp_ptr->xRight < xRmin)
X		xRmin = sp_ptr->xRight;
X	if (sp_ptr->xRight > xRmax)
X		xRmax = sp_ptr->xRight;
X
X	if (MODRES(y) == SUBYRES-1)	{	/* end of scanline */
X			/* interpolate edges to this scanline */
X		vLerp(aLeft, Vleft, VnextLeft, &VscanLeft);
X		vLerp(aRight, Vright, VnextRight, &VscanRight);
X		renderScanline(&VscanLeft, &VscanRight, y/SUBYRES, object);
X		xLmin = xRmin = MAX_X; 		/* reset extremes */
X		xLmax = xRmax = -1;
X	}
X  }
X}
X
X/*
X * Render one scanline of polygon
X */
X
XrenderScanline(Vl, Vr, y, object)
X	Vertex *Vl, *Vr; 	/* polygon vertices interpolated */
X					/* at scanline */   
X	int y;			/* scanline coordinate */
X	Surface *object;	/* shading parms for this object */
X{
X	Vertex Vpixel;	/*object info interpolated at one pixel */
X	unsigned mask[SUBYRES];	/*pixel coverage bitmask */
X	int x;			/* leftmost subpixel of current pixel */
X
X	for (x=SUBXRES*FLOOR(xLmin/SUBXRES); x<=xRmax; x+=SUBXRES) {
X		vLerp((double)(x-xLmin)/(xRmax-xLmin), Vl, Vr, &Vpixel);
X		computePixelMask(x, mask);
X		renderPixel(x/SUBXRES, y, &Vpixel,
X					/*computePixel*/Coverage(x), mask, object);
X	}
X}
X
X/*
X * Compute number of subpixels covered by polygon at current pixel
X*/
X/*computePixel*/Coverage(x)
X	int x;			/* left subpixel of pixel */
X{
X	int  area;			/* total covered area */
X	int partialArea;	  /* covered area for current subpixel y */
X	int xr = x+SUBXRES-1;	/*right subpixel of pixel */
X	int y;
X
X	/* shortcut for common case of fully covered pixel */
X	if (x>xLmax && x<xRmin)
X		return MAX_AREA;
X	
X	for (area=y=0; y<SUBYRES; y++) {
X		partialArea = MIN(sp[y].xRight, xr)
X			 - MAX(sp[y].xLeft, x) + 1;
X		if (partialArea > 0)
X			area += partialArea;
X	}
X	return area;
X}
X
X/* Compute bitmask indicating which subpixels are covered by 
X * polygon at current pixel. (Not all hidden-surface methods
X * need this mask. )
X*/
XcomputePixelMask(x, mask)
X	int x;			/* left subpixel of pixel */
X	unsigned mask[];	/* output bitmask */
X{
X	static unsigned leftMaskTable[] =
X		{ 0xFFFF, 0x7FFF, 0x3FFF, 0x1FFF, 0x0FFF, 0x07FF, 0x03FF,
X		  0x01FF, 0x00FF, 0x007F, 0x003F, 0x001F, 0x000F, 0x0007,
X		  0x0003, 0x0001  };
X	static unsigned rightMaskTable[] = 
X		{ 0x8000, 0xC000, 0xE000, 0xF000, 0xF800, 0xFC00, 
X		  0xFE00, 0xFF00, 0xFF80, 0xFFC0, 0xFFE0, 0xFFF0,
X		  0xFFF8, 0xFFFC, 0xFFFE, 0xFFFF   };
X	unsigned leftMask, rightMask; 		/* partial masks */
X	int xr = x+SUBXRES-1;			/* right subpixel of pixel */
X	int y;
X
X/* shortcut for common case of fully covered pixel */
X	if (x>xLmax && x<xRmin) 	{
X		for (y=0; y<SUBYRES; y++)
X			mask[y] = 0xFFFF;
X	} else 	{
X		for (y=0; y<SUBYRES; y++)	{
X			if (sp[y].xLeft < x) 	/* completely left of pixel*/
X				leftMask = 0xFFFF;
X			else if (sp[y].xLeft > xr)  /* completely right */	
X				leftMask = 0;
X			else
X				leftMask = leftMaskTable[sp[y].xLeft -x];
X
X			if (sp[y].xRight > xr)  	/* completely  */
X							/* right of pixel*/
X				rightMask = 0xFFFF;
X			else if (sp[y].xRight < x)	/*completely left */
X				rightMask = 0;
X			else
X				rightMask = rightMaskTable[sp[y].xRight -x];
X			mask[y] = leftMask & rightMask;
X		}
X	}
X}
END_OF_FILE
if test 7519 -ne `wc -c <'AAPolyScan.c'`; then
    echo shar: \"'AAPolyScan.c'\" unpacked with wrong size!
fi
# end of 'AAPolyScan.c'
fi
if test -f 'ConcaveScan.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'ConcaveScan.c'\"
else
echo shar: Extracting \"'ConcaveScan.c'\" \(4931 characters\)
sed "s/^X//" >'ConcaveScan.c' <<'END_OF_FILE'
X/*
X * Concave Polygon Scan Conversion
X * by Paul Heckbert
X * from "Graphics Gems", Academic Press, 1990
X */
X
X/*
X * concave: scan convert nvert-sided concave non-simple polygon with vertices at
X * (point[i].x, point[i].y) for i in [0..nvert-1] within the window win by
X * calling spanproc for each visible span of pixels.
X * Polygon can be clockwise or counterclockwise.
X * Algorithm does uniform point sampling at pixel centers.
X * Inside-outside test done by Jordan's rule: a point is considered inside if
X * an emanating ray intersects the polygon an odd number of times.
X * drawproc should fill in pixels from xl to xr inclusive on scanline y,
X * e.g:
X *	drawproc(y, xl, xr)
X *	int y, xl, xr;
X *	{
X *	    int x;
X *	    for (x=xl; x<=xr; x++)
X *		pixel_write(x, y, pixelvalue);
X *	}
X *
X *  Paul Heckbert	30 June 81, 18 Dec 89
X */
X
X#include <stdio.h>
X#include <math.h>
X#include "GraphicsGems.h"
X
X#define ALLOC(ptr, type, n)  ASSERT(ptr = (type *)malloc((n)*sizeof(type)))
X
Xtypedef struct {		/* window: a discrete 2-D rectangle */
X    int x0, y0;			/* xmin and ymin */
X    int x1, y1;			/* xmax and ymax (inclusive) */
X} Window;
X
Xtypedef struct {		/* a polygon edge */
X    double x;	/* x coordinate of edge's intersection with current scanline */
X    double dx;	/* change in x with respect to y */
X    int i;	/* edge number: edge i goes from pt[i] to pt[i+1] */
X} Edge;
X
Xstatic int n;			/* number of vertices */
Xstatic Point2 *pt;		/* vertices */
X
Xstatic int nact;		/* number of active edges */
Xstatic Edge *active;		/* active edge list:edges crossing scanline y */
X
Xint compare_ind(), compare_active();
X
Xconcave(nvert, point, win, spanproc)
Xint nvert;			/* number of vertices */
XPoint2 *point;			/* vertices of polygon */
XWindow *win;			/* screen clipping window */
Xvoid (*spanproc)();		/* called for each span of pixels */
X{
X    int k, y0, y1, y, i, j, xl, xr;
X    int *ind;		/* list of vertex indices, sorted by pt[ind[j]].y */
X
X    n = nvert;
X    pt = point;
X    if (n<=0) return;
X    ALLOC(ind, int, n);
X    ALLOC(active, Edge, n);
X
X    /* create y-sorted array of indices ind[k] into vertex list */
X    for (k=0; k<n; k++)
X	ind[k] = k;
X    qsort(ind, n, sizeof ind[0], compare_ind);	/* sort ind by pt[ind[k]].y */
X
X    nact = 0;				/* start with empty active list */
X    k = 0;				/* ind[k] is next vertex to process */
X    y0 = MAX(win->y0, ceil(pt[ind[0]].y-.5));		/* ymin of polygon */
X    y1 = MIN(win->y1, floor(pt[ind[n-1]].y-.5));	/* ymax of polygon */
X
X    for (y=y0; y<=y1; y++) {		/* step through scanlines */
X	/* scanline y is at y+.5 in continuous coordinates */
X
X	/* check vertices between previous scanline and current one, if any */
X	for (; k<n && pt[ind[k]].y<=y+.5; k++) {
X	    /* to simplify, if pt.y=y+.5, pretend it's above */
X	    /* invariant: y-.5 < pt[i].y <= y+.5 */
X	    i = ind[k];	
X	    /*
X	     * insert or delete edges before and after vertex i (i-1 to i,
X	     * and i to i+1) from active list if they cross scanline y
X	     */
X	    j = i>0 ? i-1 : n-1;	/* vertex previous to i */
X	    if (pt[j].y <= y-.5)	/* old edge, remove from active list */
X		delete(j);
X	    else if (pt[j].y > y+.5)	/* new edge, add to active list */
X		insert(j, y);
X	    j = i<n-1 ? i+1 : 0;	/* vertex next after i */
X	    if (pt[j].y <= y-.5)	/* old edge, remove from active list */
X		delete(i);
X	    else if (pt[j].y > y+.5)	/* new edge, add to active list */
X		insert(i, y);
X	}
X
X	/* sort active edge list by active[j].x */
X	qsort(active, nact, sizeof active[0], compare_active);
X
X	/* draw horizontal segments for scanline y */
X	for (j=0; j<nact; j+=2) {	/* draw horizontal segments */
X	    /* span 'tween j & j+1 is inside, span tween j+1 & j+2 is outside */
X	    xl = ceil(active[j].x-.5);		/* left end of span */
X	    if (xl<win->x0) xl = win->x0;
X	    xr = floor(active[j+1].x-.5);	/* right end of span */
X	    if (xr>win->x1) xr = win->x1;
X	    if (xl<=xr)
X		(*spanproc)(y, xl, xr);		/* draw pixels in span */
X	    active[j].x += active[j].dx;	/* increment edge coords */
X	    active[j+1].x += active[j+1].dx;
X	}
X    }
X}
X
Xstatic delete(i)		/* remove edge i from active list */
Xint i;
X{
X    int j;
X
X    for (j=0; j<nact && active[j].i!=i; j++);
X    if (j>=nact) return;	/* edge not in active list; happens at win->y0*/
X    nact--;
X    bcopy(&active[j+1], &active[j], (nact-j)*sizeof active[0]);
X}
X
Xstatic insert(i, y)		/* append edge i to end of active list */
Xint i, y;
X{
X    int j;
X    double dx;
X    Point2 *p, *q;
X
X    j = i<n-1 ? i+1 : 0;
X    if (pt[i].y < pt[j].y) {p = &pt[i]; q = &pt[j];}
X    else		   {p = &pt[j]; q = &pt[i];}
X    /* initialize x position at intersection of edge with scanline y */
X    active[nact].dx = dx = (q->x-p->x)/(q->y-p->y);
X    active[nact].x = dx*(y+.5-p->y)+p->x;
X    active[nact].i = i;
X    nact++;
X}
X
X/* comparison routines for qsort */
Xcompare_ind(u, v) int *u, *v; {return pt[*u].y <= pt[*v].y ? -1 : 1;}
Xcompare_active(u, v) Edge *u, *v; {return u->x <= v->x ? -1 : 1;}
END_OF_FILE
if test 4931 -ne `wc -c <'ConcaveScan.c'`; then
    echo shar: \"'ConcaveScan.c'\" unpacked with wrong size!
fi
# end of 'ConcaveScan.c'
fi
if test -f 'Forms.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Forms.c'\"
else
echo shar: Extracting \"'Forms.c'\" \(5341 characters\)
sed "s/^X//" >'Forms.c' <<'END_OF_FILE'
X/*
XForms, Vectors,and Transforms
Xby Bob Wallis
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X/*----------------------------------------------------------------------
XThe main program below is set up to solve the Bezier subdivision problem
Xin "Forms, Vectors, and Transforms".  The subroutines are useful in
Xsolving general problems which require manipulating matrices via exact
Xinteger arithmetic.  The intended application is validating or avoiding
Xtedious algebraic calculations.  As such, no thought was given to
Xefficiency.  
X----------------------------------------------------------------------*/
X#define ABS(x) ((x)>(0)? (x):(-x))
X#define N 4			/* size of matrices to deal with */
Xint     M[N][N] =		/* Bezier weights */
X{
X      1,   0,   0,   0,
X     -3,   3,   0,   0,
X      3,  -6,   3,   0,
X     -1,   3,  -3,   1,
X};
Xint     T[N][N] =	/* re-parameterization xform for top half */
X{
X      1,  -1,   1,  -1,
X      0,   2,  -4,   6,
X      0,   0,   4, -12,
X      0,   0,   0,   8
X};
Xmain ()
X{
X    int     i,
X            j,
X            scale,
X            gcd,
X            C[N][N],
X            S[N][N],
X            Madj[N][N],
X            Tadj[N][N],
X            Mdet,
X            Tdet;
X
X
X    Tdet = adjoint (T, Tadj);   /* inverse without division by */
X    Mdet = adjoint (M, Madj);   /* determinant of T and M */
X    matmult (Madj, Tadj, C);
X    matmult (C, M, S);		/* Madj*Tadj*M -> S */
X    scale = gcd = Mdet * Tdet;  /* scale factors of both determinants */
X    for (i = 0; i < N; i++)	/* find the greatest common */
X    {				/* demoninator of S and determinants */
X		for (j = 0; j < N; j++)
X	    	gcd = Gcd (gcd, S[i][j]);
X    }
X    scale /= gcd;		/* divide everything by gcd to get */
X    for (i = 0; i < N; i++)	/* matrix and scale factor in lowest */
X    {				/* integer terms possible */
X		for (j = 0; j < N; j++)
X	    	S[i][j] /= gcd;
X    }
X    printf ("scale factor = 1/%d  ", scale);
X    print_mat ("M=", M, N);     /* display the results */
X    print_mat ("T=", T, N);
X    print_mat ("S=", S, N);     /* subdivision matrix */
X    exit (0);
X}
XGcd (a, b)			/*returns greatest common demoninator */
Xint     a,			/* of (a,b) */
X       	b;
X{
X    int		i,
X            	r;
X    a = ABS (a);		/* force positive */
X    b = ABS (b);
X    if (a < b)			/* exchange so that a >= b */
X    {
X		i = b;
X		b = a;
X		a = i;
X    }
X    if 	(b == 0)
X		return (a);	/* finished */
X    r = a % b;			/* remainder */
X    if (r == 0)
X		return (b);	/* finished */
X    else			/* recursive call */
X		return (Gcd (b, r));
X}
X
Xadjoint (A, Aadj)			/* returns determinant of A */
Xint     A[N][N],			/* input matrix */
X        Aadj[N][N];			/* output = adjoint of A */
X{					/* must have N >= 3 */
X    int	   i,
X            j,
X            I[N],			/* arrays of row and column indices */
X            J[N],
X            Isub[N],			/* sub-arrays of the above */
X            Jsub[N],
X            cofactor,
X            det;
X    if (N < 3)
X    {
X		printf ("must have N >= 3\n");
X		exit (1);
X    }
X    for (i = 0; i < N; i++)
X    {					/* lookup tables to select a */
X		I[i] = i;		/* particular subset of */
X		J[i] = i;		/* rows and columns */
X    }
X    det = 0;
X    for (i = 0; i < N; i++)
X    {					/* delete ith row */
X		subarray (I, Isub, N, i);
X		for (j = 0; j < N; j++)
X		{				/* delete jth column */
X	    		subarray (J, Jsub, N, j);
X	    		cofactor = determinant (A, Isub, Jsub, N - 1,(i + j) & 1);
X	    		if (j == 0)		/* use 0th column for det */
X				det += cofactor * A[i][0];
X			Aadj[j][i] = cofactor;
X		}
X   	}
X    return (det);
X}
Xdeterminant (A, I, J, n, parity)/* actually gets a sub-determinant */
Xint     A[N][N],			/* input = entire matrix */
X       	I[N],				/* row sub-array we want */
X        J[N],				/* col sub-array we want */
X        parity,				/* 1-> flip polarity */
X        n;				/* # elements in subarrays */
X
X{
X    int     i,
X            j,
X            det,
X            j_,
X            Jsub[N];
X    if (n <= 2)			/* call ourselves till we get down to */
X    {				/* a 2x2 matrix */
X		det =
X	    	(A[I[0]][J[0]] * A[I[1]][J[1]]) -
X	    	(A[I[1]][J[0]] * A[I[0]][J[1]]);
X		if (parity)
X	    det = -det;
X		return (det);
X    }					/* if (n <= 2) */
X    det = 0;				/* n > 2; call recursively */
X    i = I[0];				/* strike out 0th row */
X    for (j_ = 0; j_ < n; j_++)
X    {					/* strike out jth column */
X		subarray (J, Jsub, n, j_);
X		j = J[j_];		/* I + 1 => struck out 0th row */
X		det += A[i][j] * determinant (A, I + 1, Jsub, n - 1, j_ & 1);
X    }
X    if (parity)
X		det = -det;
X    return (det);
X}
Xsubarray (src, dest, n, k)		/* strike out kth row/column */
Xint  	*src,				/* 	source array of n indices */
X     	*dest,				/* dest array formed by deleting k  */
X       	 n,
X       	 k;
X{
X    int     i;
X    for (i = 0; i < n; i++, src++)
X		if (i != k)			/* skip over k */
X	    	*dest++ = *src;
X}
X
Xmatmult (A, B, C)				/* C = A*B */
Xint     A[N][N],
X        B[N][N],
X        C[N][N];
X{
X    int		i,
X		j,
X		k,
X		sum;
X    for (i = 0; i < N; i++)
X    {
X		for (k = 0; k < N; k++)
X		{
X	    	sum = 0;
X	    	for (j = 0; j < N; j++)
X				sum += A[i][j] * B[j][k];
X	    	C[i][k] = sum;
X		}
X    }
X}
Xprint_mat (string, mat, n)
Xchar   *string;
Xint     mat[N][N],
X         n;
X{
X    int     i,
X             j;
X    printf ("%s\n", string);
X    for (i = 0; i < n; i++)
X    {
X		for (j = 0; j < n; j++)
X	    printf (" %8ld", mat[i][j]);
X	printf ("\n");
X    }
X}
END_OF_FILE
if test 5341 -ne `wc -c <'Forms.c'`; then
    echo shar: \"'Forms.c'\" unpacked with wrong size!
fi
# end of 'Forms.c'
fi
if test -f 'Interleave.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Interleave.c'\"
else
echo shar: Extracting \"'Interleave.c'\" \(6497 characters\)
sed "s/^X//" >'Interleave.c' <<'END_OF_FILE'
X/*
XBit Interleaving for Quad- or Octrees
Xby Clifford A. Shaffer
Xfrom "Graphics Gems", Academic Press, 1990
X*/
X
X#include "GraphicsGems.h"
X#define B_MAX_DEPTH 14	/* maximum depth allowed */
X
X/* byteval is the lookup table for coordinate interleaving.  Given a
X   4 bit portion of the (x, y) coordinates, return the bit interleaving.
X   Notice that this table looks like the order in which the pixels of
X   a 16 X 16 pixel image would be visited. */
Xint byteval[16][16] = 
X    {  0,  1,  4,  5, 16, 17, 20, 21, 64, 65, 68, 69, 80, 81, 84, 85,
X       2,  3,  6,  7, 18, 19, 22, 23, 66, 67, 70, 71, 82, 83, 86, 87,
X       8,  9, 12, 13, 24, 25, 28, 29, 72, 73, 76, 77, 88, 89, 92, 93,
X      10, 11, 14, 15, 26, 27, 30, 31, 74, 75, 78, 79, 90, 91, 94, 95,
X      32, 33, 36, 37, 48, 49, 52, 53, 96, 97,100,101,112,113,116,117,
X      34, 35, 38, 39, 50, 51, 54, 55, 98, 99,102,103,114,115,118,119,
X      40, 41, 44, 45, 56, 57, 60, 61,104,105,108,109,120,121,124,125,
X      42, 43, 46, 47, 58, 59, 62, 63,106,107,110,111,122,123,126,127,
X     128,129,132,133,144,145,148,149,192,193,196,197,208,209,212,213,
X     130,131,134,135,146,147,150,151,194,195,198,199,210,211,214,215,
X     136,137,140,141,152,153,156,157,200,201,204,205,216,217,220,221,
X     138,139,142,143,154,155,158,159,202,203,206,207,218,219,222,223,
X     160,161,164,165,176,177,180,181,224,225,228,229,240,241,244,245,
X     162,163,166,167,178,179,182,183,226,227,230,231,242,243,246,247,
X     168,169,172,173,184,185,188,189,232,233,236,237,248,249,252,253,
X     170,171,174,175,186,187,190,191,234,235,238,239,250,251,254,255};
X
X/* bytemask is the mask for byte interleaving - masks out the 
X   non-significant bit positions.  This is determined by the
X   depth of the node. For example, a node of depth 0 is at the root.
X  Thus, there are no branchs and no bits are significant. 
X  The bottom 4 bits (the depth) are always retained. 
X  Values are in octal notation. */
Xint bytemask[B_MAX_DEPTH + 1] = {017,
X     030000000017,036000000017,037400000017,037700000017,
X     037760000017,037774000017,037777000017,037777600017,
X     037777740017,037777770017,037777776017,037777777417,
X     037777777717,037777777777};
X
X
Xlong *interleave(addr, x, y, depth, max_depth)
X/* Return the interleaved code for a quadtree node at depth depth 
Xwhose upper left hand corner has coordinates (x, y) in a tree with maximum
Xdepth max_depth.  This function receives and returns a 
Xpointer to addr, which is either a long interger or (more typically)
Xan array of long integers whose first integer contains the 
Xinterleaved code. */
X long *addr;
X int max_depth, depth;
X int x, y;
X{
X
X/* Scale x, y values to be consistent with maximum coord size */
X/* and depth of tree */
X x <<= (B_MAX_DEPTH - max_depth);
X y <<= (B_MAX_DEPTH - max_depth);
X
X/* calculate the bit interleaving of the x, y values that have now
X   been appropriately shifted, and place this interleave in the address
X   portion of addr.  Note that the binary representations of x and y are
X   being processed from right to left */
X
X *addr = depth;
X *addr |= byteval[y & 03][x & 03] << 4;
X *addr |= byteval[(y >> 2) & 017][(x >> 2) & 017] << 8;
X *addr |= byteval[(y >> 6) & 017][(x >> 6) & 017] << 16;
X *addr |= byteval[(y >> 10) & 017][(x >> 10) & 017] << 24;
X *addr &= bytemask[depth];
X
X/* if there were unused portions of the x and y addresses then  */
X/* the address was too large for the depth values given.  */
X/*  Return address built */
X return (addr);
X}
X
X
X
X/* The next two arrays are used in calculating the (x, y) coordinates
X   of the upper left-hand corner of a node from its bit interleaved
X   address.  Given an 8 bit number, the arrays return the effect of
X   removing every other bit (the y bits preceed the x bits). */
X
Xint xval[256] = { 0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
X		         4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
X		         0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
X		         4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
X		         8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
X		        12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
X		         8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
X		        12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
X		         0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
X		         4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
X		         0, 1, 0, 1, 2, 3, 2, 3, 0, 1, 0, 1, 2, 3, 2, 3,
X		         4, 5, 4, 5, 6, 7, 6, 7, 4, 5, 4, 5, 6, 7, 6, 7,
X		         8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
X		        12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15,
X		         8, 9, 8, 9,10,11,10,11, 8, 9, 8, 9,10,11,10,11,
X		        12,13,12,13,14,15,14,15,12,13,12,13,14,15,14,15};
X
X
Xint yval[256] = { 0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
X		      0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
X		      4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
X		      4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
X		      0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
X		      0, 0, 1, 1, 0, 0, 1, 1, 2, 2, 3, 3, 2, 2, 3, 3,
X		      4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
X		      4, 4, 5, 5, 4, 4, 5, 5, 6, 6, 7, 7, 6, 6, 7, 7,
X		      8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
X		      8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
X		     12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
X		     12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
X		      8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
X		      8, 8, 9, 9, 8, 8, 9, 9,10,10,11,11,10,10,11,11,
X		     12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15,
X		     12,12,13,13,12,12,13,13,14,14,15,15,14,14,15,15};
X
X
X
X
X
Xint getx(addr, max_depth)
X/* Return the x coordinate of the upper left hand corner of addr for a
X   tree with maximum depth max_depth. */
X long *addr;
X int max_depth;
X{
X register x;
X
X x = xval[(*addr >> 4) & 017];
X x |= xval[(*addr >> 8) & 0377] << 2;
X x |= xval[(*addr >> 16) & 0377] << 6;
X x |= xval[(*addr >> 24) & 0377] << 10;
X x >>= B_MAX_DEPTH - max_depth;
X return (x);
X}
X
X
X
Xint QKy(addr, max_depth)
X/* Return the y coordinate of the upper left hand corner of addr for a
X   tree with maximum depth max_depth. */
X
X long *addr;
X int max_depth;
X{
X register y;
X
X y = yval[(*addr >> 4) & 017];
X y |= yval[(*addr >> 8) & 0377] << 2;
X y |= yval[(*addr >> 16) & 0377] << 6;
X y |= yval[(*addr >> 24) & 0377] << 10;
X y >>= B_MAX_DEPTH - max_depth;
X return (y);
X}
X
Xint getdepth(addr)
X/* Return the depth of the node.  Simply return the bottom 4 bits. */
X
X long *addr;
X{
X
X return(*addr & 017);
X}
END_OF_FILE
if test 6497 -ne `wc -c <'Interleave.c'`; then
    echo shar: \"'Interleave.c'\" unpacked with wrong size!
fi
# end of 'Interleave.c'
fi
if test -f 'Sturm/sturm.c' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'Sturm/sturm.c'\"
else
echo shar: Extracting \"'Sturm/sturm.c'\" \(5198 characters\)
sed "s/^X//" >'Sturm/sturm.c' <<'END_OF_FILE'
X
X/*
X * sturm.c
X *
X *	the functions to build and evaluate the Sturm sequence
X */
X#include <math.h>
X#include <stdio.h>
X#include "solve.h"
X
X/*
X * modp
X *
X *	calculates the modulus of u(x) / v(x) leaving it in r, it
X *  returns 0 if r(x) is a constant.
X *  note: this function assumes the leading coefficient of v 
X *	is 1 or -1
X */
Xstatic int
Xmodp(u, v, r)
X	poly	*u, *v, *r;
X{
X	int		k, j;
X	double	*nr, *end, *uc;
X
X	nr = r->coef;
X	end = &u->coef[u->ord];
X
X	uc = u->coef;
X	while (uc <= end)
X			*nr++ = *uc++;
X
X	if (v->coef[v->ord] < 0.0) {
X
X
X			for (k = u->ord - v->ord - 1; k >= 0; k -= 2)
X				r->coef[k] = -r->coef[k];
X
X			for (k = u->ord - v->ord; k >= 0; k--)
X				for (j = v->ord + k - 1; j >= k; j--)
X					r->coef[j] = -r->coef[j] - r->coef[v->ord + k]
X					* v->coef[j - k];
X	} else {
X			for (k = u->ord - v->ord; k >= 0; k--)
X				for (j = v->ord + k - 1; j >= k; j--)
X				r->coef[j] -= r->coef[v->ord + k] * v->coef[j - k];
X	}
X
X	k = v->ord - 1;
X	while (k >= 0 && fabs(r->coef[k]) < SMALL_ENOUGH) {
X		r->coef[k] = 0.0;
X		k--;
X	}
X
X	r->ord = (k < 0) ? 0 : k;
X
X	return(r->ord);
X}
X
X/*
X * buildsturm
X *
X *	build up a sturm sequence for a polynomial in smat, returning
X * the number of polynomials in the sequence
X */
Xint
Xbuildsturm(ord, sseq)
X	int		ord;
X	poly	*sseq;
X{
X	int		i;
X	double	f, *fp, *fc;
X	poly	*sp;
X
X	sseq[0].ord = ord;
X	sseq[1].ord = ord - 1;
X
X
X	/*
X	 * calculate the derivative and normalise the leading
X	 * coefficient.
X	 */
X	f = fabs(sseq[0].coef[ord] * ord);
X	fp = sseq[1].coef;
X	fc = sseq[0].coef + 1;
X	for (i = 1; i <= ord; i++)
X			*fp++ = *fc++ * i / f;
X
X	/*
X	 * construct the rest of the Sturm sequence
X	 */
X	for (sp = sseq + 2; modp(sp - 2, sp - 1, sp); sp++) {
X
X		/*
X		 * reverse the sign and normalise
X		 */
X		f = -fabs(sp->coef[sp->ord]);
X		for (fp = &sp->coef[sp->ord]; fp >= sp->coef; fp--)
X				*fp /= f;
X	}
X
X	sp->coef[0] = -sp->coef[0];	/* reverse the sign */
X
X	return(sp - sseq);
X}
X
X/*
X * numroots
X *
X *	return the number of distinct real roots of the polynomial
X * described in sseq.
X */
Xint
Xnumroots(np, sseq, atneg, atpos)
X		int		np;
X		poly	*sseq;
X		int		*atneg, *atpos;
X{
X		int		atposinf, atneginf;
X		poly	*s;
X		double	f, lf;
X
X		atposinf = atneginf = 0;
X
X
X	/*
X	 * changes at positve infinity
X	 */
X	lf = sseq[0].coef[sseq[0].ord];
X
X	for (s = sseq + 1; s <= sseq + np; s++) {
X			f = s->coef[s->ord];
X			if (lf == 0.0 || lf * f < 0)
X				atposinf++;
X		lf = f;
X	}
X
X	/*
X	 * changes at negative infinity
X	 */
X	if (sseq[0].ord & 1)
X			lf = -sseq[0].coef[sseq[0].ord];
X	else
X			lf = sseq[0].coef[sseq[0].ord];
X
X	for (s = sseq + 1; s <= sseq + np; s++) {
X			if (s->ord & 1)
X				f = -s->coef[s->ord];
X			else
X				f = s->coef[s->ord];
X			if (lf == 0.0 || lf * f < 0)
X				atneginf++;
X			lf = f;
X	}
X
X	*atneg = atneginf;
X	*atpos = atposinf;
X
X	return(atneginf - atposinf);
X}
X
X/*
X * numchanges
X *
X *	return the number of sign changes in the Sturm sequence in
X * sseq at the value a.
X */
Xint
Xnumchanges(np, sseq, a)
X	int		np;
X	poly	*sseq;
X	double	a;
X
X{
X	int		changes;
X	double	f, lf;
X	poly	*s;
X
X	changes = 0;
X
X	lf = evalpoly(sseq[0].ord, sseq[0].coef, a);
X
X	for (s = sseq + 1; s <= sseq + np; s++) {
X			f = evalpoly(s->ord, s->coef, a);
X			if (lf == 0.0 || lf * f < 0)
X				changes++;
X			lf = f;
X	}
X
X	return(changes);
X}
X
X/*
X * sbisect
X *
X *	uses a bisection based on the sturm sequence for the polynomial
X * described in sseq to isolate intervals in which roots occur,
X * the roots are returned in the roots array in order of magnitude.
X */
Xsbisect(np, sseq, min, max, atmin, atmax, roots)
X	int		np;
X	poly	*sseq;
X	double	min, max;
X	int		atmin, atmax;
X	double	*roots;
X{
X	double	mid;
X	int		n1, n2, its, atmid, nroot;
X
X	if ((nroot = atmin - atmax) == 1) {
X
X		/*
X		 * first try a less expensive technique.
X		 */
X		if (modrf(sseq->ord, sseq->coef, min, max, &roots[0]))
X			return;
X
X
X		/*
X		 * if we get here we have to evaluate the root the hard
X		 * way by using the Sturm sequence.
X		 */
X		for (its = 0; its < MAXIT; its++) {
X				mid = (min + max) / 2;
X
X				atmid = numchanges(np, sseq, mid);
X
X				if (fabs(mid) > RELERROR) {
X					if (fabs((max - min) / mid) < RELERROR) {
X						roots[0] = mid;
X						return;
X					}
X				} else if (fabs(max - min) < RELERROR) {
X					roots[0] = mid;
X					return;
X				}
X
X				if ((atmin - atmid) == 0)
X					min = mid;
X				else
X					max = mid;
X			}
X
X		if (its == MAXIT) {
X				fprintf(stderr, "sbisect: overflow min %f max %f\
X					diff %e nroot %d n1 %d n2 %d\n",
X					min, max, max - min, nroot, n1, n2);
X			roots[0] = mid;
X		}
X
X		return;
X	}
X
X	/*
X	 * more than one root in the interval, we have to bisect...
X	 */
X	for (its = 0; its < MAXIT; its++) {
X
X			mid = (min + max) / 2;
X
X			atmid = numchanges(np, sseq, mid);
X
X			n1 = atmin - atmid;
X			n2 = atmid - atmax;
X
X
X			if (n1 != 0 && n2 != 0) {
X				sbisect(np, sseq, min, mid, atmin, atmid, roots);
X				sbisect(np, sseq, mid, max, atmid, atmax, &roots[n1]);
X				break;
X			}
X
X			if (n1 == 0)
X				min = mid;
X			else
X				max = mid;
X	}
X
X	if (its == MAXIT) {
X			fprintf(stderr, "sbisect: roots too close together\n");
X			fprintf(stderr, "sbisect: overflow min %f max %f diff %e\
X				nroot %d n1 %d n2 %d\n",
X				min, max, max - min, nroot, n1, n2);
X			for (n1 = atmax; n1 < atmin; n1++)
X			roots[n1 - atmax] = mid;
X	}
X}
END_OF_FILE
if test 5198 -ne `wc -c <'Sturm/sturm.c'`; then
    echo shar: \"'Sturm/sturm.c'\" unpacked with wrong size!
fi
# end of 'Sturm/sturm.c'
fi
echo shar: End of archive 4 \(of 5\).
cp /dev/null ark4isdone
MISSING=""
for I in 1 2 3 4 5 ; do
    if test ! -f ark${I}isdone ; then
	MISSING="${MISSING} ${I}"
    fi
done
if test "${MISSING}" = "" ; then
    echo You have unpacked all 5 archives.
    rm -f ark[1-9]isdone
else
    echo You still need to unpack the following archives:
    echo "        " ${MISSING}
fi
##  End of shell archive.
exit 0

  .