From 285f5bec5bb03d4e825e5d866e94008088dd6155 Mon Sep 17 00:00:00 2001
From: xleroy <xleroy@fca1b0fc-160b-0410-b1d3-a4f43f01ea2e>
Date: Sat, 9 Aug 2008 08:06:33 +0000
Subject: [PATCH] Ajout nouveaux tests

git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@708 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e
---
 test/compression/.depend     |   15 +
 test/compression/Makefile    |   59 ++
 test/compression/arcode.c    |  926 ++++++++++++++++++++++++++++++
 test/compression/arcode.h    |   69 +++
 test/compression/armain.c    |  266 +++++++++
 test/compression/bitfile.c   | 1043 ++++++++++++++++++++++++++++++++++
 test/compression/bitfile.h   |  121 ++++
 test/compression/lzdecode.c  |  282 +++++++++
 test/compression/lzencode.c  |  300 ++++++++++
 test/compression/lzhash.c    |  358 ++++++++++++
 test/compression/lzlocal.h   |  123 ++++
 test/compression/lzss.h      |   96 ++++
 test/compression/lzssmain.c  |  278 +++++++++
 test/compression/lzvars.c    |   76 +++
 test/compression/lzw.h       |   70 +++
 test/compression/lzwdecode.c |  306 ++++++++++
 test/compression/lzwencode.c |  465 +++++++++++++++
 test/compression/lzwmain.c   |  260 +++++++++
 test/compression/optlist.c   |  228 ++++++++
 test/compression/optlist.h   |   74 +++
 test/raytracer/.depend       |   20 +
 test/raytracer/Makefile      |   43 ++
 test/raytracer/arrays.c      |   30 +
 test/raytracer/arrays.h      |   31 +
 test/raytracer/config.h      |   27 +
 test/raytracer/eval.c        |  672 ++++++++++++++++++++++
 test/raytracer/eval.h        |   18 +
 test/raytracer/gml.h         |   73 +++
 test/raytracer/gmllexer.c    |  297 ++++++++++
 test/raytracer/gmllexer.h    |   21 +
 test/raytracer/gmlparser.c   |  102 ++++
 test/raytracer/gmlparser.h   |    3 +
 test/raytracer/intersect.c   |  416 ++++++++++++++
 test/raytracer/intersect.h   |   11 +
 test/raytracer/kal.gml       |   78 +++
 test/raytracer/light.c       |  126 ++++
 test/raytracer/light.h       |   33 ++
 test/raytracer/main.c        |   13 +
 test/raytracer/matrix.c      |  102 ++++
 test/raytracer/matrix.h      |   18 +
 test/raytracer/memory.c      |   60 ++
 test/raytracer/object.c      |  214 +++++++
 test/raytracer/object.h      |   36 ++
 test/raytracer/point.h       |   18 +
 test/raytracer/render.c      |  128 +++++
 test/raytracer/render.h      |    9 +
 test/raytracer/simplify.c    |  177 ++++++
 test/raytracer/simplify.h    |    3 +
 test/raytracer/surface.c     |  148 +++++
 test/raytracer/surface.h     |    7 +
 test/raytracer/vector.c      |   74 +++
 test/raytracer/vector.h      |   23 +
 52 files changed, 8446 insertions(+)
 create mode 100644 test/compression/.depend
 create mode 100644 test/compression/Makefile
 create mode 100755 test/compression/arcode.c
 create mode 100755 test/compression/arcode.h
 create mode 100755 test/compression/armain.c
 create mode 100644 test/compression/bitfile.c
 create mode 100644 test/compression/bitfile.h
 create mode 100644 test/compression/lzdecode.c
 create mode 100644 test/compression/lzencode.c
 create mode 100644 test/compression/lzhash.c
 create mode 100644 test/compression/lzlocal.h
 create mode 100644 test/compression/lzss.h
 create mode 100644 test/compression/lzssmain.c
 create mode 100644 test/compression/lzvars.c
 create mode 100644 test/compression/lzw.h
 create mode 100644 test/compression/lzwdecode.c
 create mode 100644 test/compression/lzwencode.c
 create mode 100644 test/compression/lzwmain.c
 create mode 100644 test/compression/optlist.c
 create mode 100644 test/compression/optlist.h
 create mode 100644 test/raytracer/.depend
 create mode 100644 test/raytracer/Makefile
 create mode 100644 test/raytracer/arrays.c
 create mode 100644 test/raytracer/arrays.h
 create mode 100644 test/raytracer/config.h
 create mode 100644 test/raytracer/eval.c
 create mode 100644 test/raytracer/eval.h
 create mode 100644 test/raytracer/gml.h
 create mode 100644 test/raytracer/gmllexer.c
 create mode 100644 test/raytracer/gmllexer.h
 create mode 100644 test/raytracer/gmlparser.c
 create mode 100644 test/raytracer/gmlparser.h
 create mode 100644 test/raytracer/intersect.c
 create mode 100644 test/raytracer/intersect.h
 create mode 100644 test/raytracer/kal.gml
 create mode 100644 test/raytracer/light.c
 create mode 100644 test/raytracer/light.h
 create mode 100644 test/raytracer/main.c
 create mode 100644 test/raytracer/matrix.c
 create mode 100644 test/raytracer/matrix.h
 create mode 100644 test/raytracer/memory.c
 create mode 100644 test/raytracer/object.c
 create mode 100644 test/raytracer/object.h
 create mode 100644 test/raytracer/point.h
 create mode 100644 test/raytracer/render.c
 create mode 100644 test/raytracer/render.h
 create mode 100644 test/raytracer/simplify.c
 create mode 100644 test/raytracer/simplify.h
 create mode 100644 test/raytracer/surface.c
 create mode 100644 test/raytracer/surface.h
 create mode 100644 test/raytracer/vector.c
 create mode 100644 test/raytracer/vector.h

diff --git a/test/compression/.depend b/test/compression/.depend
new file mode 100644
index 000000000..e0dcdaec1
--- /dev/null
+++ b/test/compression/.depend
@@ -0,0 +1,15 @@
+arcode.o: arcode.c arcode.h bitfile.h
+arcode.light.o: arcode.light.c
+armain.o: armain.c optlist.h arcode.h
+bitfile.o: bitfile.c bitfile.h
+bitfile.light.o: bitfile.light.c
+lzdecode.o: lzdecode.c lzlocal.h bitfile.h
+lzencode.o: lzencode.c lzlocal.h bitfile.h
+lzhash.o: lzhash.c lzlocal.h
+lzssmain.o: lzssmain.c lzss.h optlist.h
+lzvars.o: lzvars.c lzlocal.h
+lzwdecode.o: lzwdecode.c lzw.h bitfile.h
+lzwencode.o: lzwencode.c lzw.h bitfile.h
+lzwmain.o: lzwmain.c optlist.h lzw.h
+optlist.o: optlist.c optlist.h
+optlist.light.o: optlist.light.c
diff --git a/test/compression/Makefile b/test/compression/Makefile
new file mode 100644
index 000000000..ba83c871c
--- /dev/null
+++ b/test/compression/Makefile
@@ -0,0 +1,59 @@
+CC=../../ccomp 
+CFLAGS=-U__GNUC__ -stdlib ../../runtime -dclight -dasm
+LIBS=
+TIME=xtime -o /dev/null -mintime 1.0
+
+EXE=arcode lzw lzss
+
+COMMON_OBJS=optlist.o bitfile.o
+
+all: $(EXE)
+
+ARCODE_OBJS=$(COMMON_OBJS) arcode.o armain.o
+
+arcode: $(ARCODE_OBJS)
+	$(CC) $(CFLAGS) -o $@ $(ARCODE_OBJS) $(LIBS)
+
+LZW_OBJS=$(COMMON_OBJS) lzwencode.o lzwdecode.o lzwmain.o
+
+lzw: $(LZW_OBJS)
+	$(CC) $(CFLAGS) -o $@ $(LZW_OBJS) $(LIBS)
+
+LZSS_OBJS=$(COMMON_OBJS) lzvars.o lzhash.o lzencode.o lzdecode.o lzssmain.o
+
+lzss: $(LZSS_OBJS)
+	$(CC) $(CFLAGS) -o $@ $(LZSS_OBJS) $(LIBS)
+
+TESTFILE=/mach_kernel
+TESTCOMPR=/tmp/testcompr.out
+TESTEXPND=/tmp/testexpnd.out
+
+test:
+	rm -f $(TESTCOMPR) $(TESTEXPND)
+	@for i in $(EXE); do \
+          echo "$$i: compression..."; \
+          ./$$i -c -i $(TESTFILE) -o $(TESTCOMPR); \
+          echo "$$i: decompression..."; \
+          ./$$i -d -i $(TESTCOMPR) -o $(TESTEXPND); \
+          if cmp $(TESTFILE) $(TESTEXPND); \
+          then echo "$$i: passed"; \
+          else echo "$$i: FAILED"; \
+          fi; \
+        done
+	rm -f $(TESTCOMPR) $(TESTEXPND)
+
+bench:
+	rm -f $(TESTCOMPR)
+	@for i in $(EXE); do \
+          echo -n "$$i "; \
+          $(TIME) sh -c "./$$i -c -i $(TESTFILE) -o $(TESTCOMPR) && ./$$i -d -i $(TESTCOMPR) -o /dev/null"; \
+         done
+	rm -f $(TESTCOMPR)
+
+include .depend
+
+clean:
+	rm -f *.o *.light.c *.s $(EXE)
+
+depend:
+	gcc -MM *.c > .depend
diff --git a/test/compression/arcode.c b/test/compression/arcode.c
new file mode 100755
index 000000000..f915cc257
--- /dev/null
+++ b/test/compression/arcode.c
@@ -0,0 +1,926 @@
+/***************************************************************************
+*                 Arithmetic Encoding and Decoding Library
+*
+*   File    : arcode.c
+*   Purpose : Use arithmetic coding to compress/decompress file streams
+*   Author  : Michael Dipperstein
+*   Date    : April 2, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: arcode.c,v 1.5 2007/09/08 15:47:02 michael Exp $
+*   $Log: arcode.c,v $
+*   Revision 1.5  2007/09/08 15:47:02  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.4  2006/03/02 06:43:37  michael
+*   Expanded tabs
+*
+*   Revision 1.3  2006/01/12 07:39:24  michael
+*   Use BitFileGetBitsIntBit and FilePutBitsInt for reading and writing the
+*   header.  This makes the code a little cleaner, but the new header is not
+*   compatible with the old header.
+*
+*   Revision 1.2  2004/08/13 13:10:27  michael
+*   Add support for adaptive encoding
+*
+*   Use binary search when trying to find decoded symbol
+*
+*   Revision 1.1.1.1  2004/04/04 14:54:13  michael
+*   Initial version
+*
+****************************************************************************
+*
+* Arcode: An ANSI C Arithmetic Encoding/Decoding Routines
+* Copyright (C) 2004, 2006-2007 by Michael Dipperstein (mdipper@cs.ucsb.edu)
+*
+* This file is part of the arcode library.
+*
+* The arcode library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The arcode library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+#include "arcode.h"
+#include "bitfile.h"
+
+/* compile-time options */
+#undef BUILD_DEBUG_OUTPUT                   /* debugging output */
+
+#if !(USHRT_MAX < ULONG_MAX)
+#error "Implementation requires USHRT_MAX < ULONG_MAX"
+#endif
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+typedef unsigned short probability_t;       /* probability count type */
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#define EOF_CHAR    (UCHAR_MAX + 1)
+
+/* number of bits used to compute running code values */
+#define PRECISION           (8 * sizeof(probability_t))
+
+/* 2 bits less than precision. keeps lower and upper bounds from crossing. */
+#define MAX_PROBABILITY     (1 << (PRECISION - 2))
+
+/***************************************************************************
+*                                  MACROS
+***************************************************************************/
+/* set bit x to 1 in probability_t.  Bit 0 is MSB */
+#define MASK_BIT(x) (probability_t)(1 << (PRECISION - (1 + (x))))
+
+/* indices for a symbol's lower and upper cumulative probability ranges */
+#define LOWER(c)        (c)
+#define UPPER(c)    ((c) + 1)
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+probability_t lower;          /* lower bound of current code range */
+probability_t upper;          /* upper bound of current code range */
+
+probability_t code;           /* current MSBs of encode input stream */
+
+unsigned char underflowBits;      /* current underflow bit count */
+
+/* probability ranges for each symbol: [ranges[LOWER(c)], ranges[UPPER(c)]) */
+probability_t ranges[UPPER(EOF_CHAR) + 1];
+probability_t cumulativeProb;   /* cumulative probability  of all ranges */
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+/* read write file headers */
+void WriteHeader(bit_file_t *bfpOut);
+int ReadHeader(bit_file_t *bfpIn);
+
+/* applies symbol's ranges to current upper and lower range bounds */
+void ApplySymbolRange(int symbol, char staticModel);
+
+/* routines for encoding*/
+void WriteEncodedBits(bit_file_t *bfpOut);
+void WriteRemaining(bit_file_t *bfpOut);
+int BuildProbabilityRangeList(FILE *fpIn);
+void InitializeAdaptiveProbabilityRangeList(void);
+
+/* routines for decoding */
+void InitializeDecoder(bit_file_t *bfpOut, char staticModel);
+probability_t GetUnscaledCode(void);
+int GetSymbolFromProbability(probability_t probability);
+void ReadEncodedBits(bit_file_t *bfpIn);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/***************************************************************************
+*   Function   : ArEncodeFile
+*   Description: This routine generates a list of arithmetic code ranges for
+*                a file and then uses them to write out an encoded version
+*                of that file.
+*   Parameters : inFile - Name of file to encode
+*                outFile - Name of file to write encoded output to
+*                staticModel - TRUE if encoding with a static model
+*   Effects    : File is arithmetically encoded
+*   Returned   : TRUE for success, otherwise FALSE.
+***************************************************************************/
+int ArEncodeFile(char *inFile, char *outFile, char staticModel)
+{
+    int c;
+    FILE *fpIn;                         /* uncoded input */
+    bit_file_t *bfpOut;                 /* encoded output */
+
+    /* open input and output files */
+    if ((fpIn = fopen(inFile, "rb")) == NULL)
+    {
+        perror(inFile);
+        return FALSE;
+    }
+
+    if (outFile == NULL)
+    {
+        bfpOut = MakeBitFile(stdout, BF_WRITE);
+    }
+    else
+    {
+        if ((bfpOut = BitFileOpen(outFile, BF_WRITE)) == NULL)
+        {
+            fclose(fpIn);
+            perror(outFile);
+            return FALSE;
+        }
+    }
+
+    if (staticModel)
+    {
+        /* count symbols in file and come up with a list of probability ranges */
+        if (!BuildProbabilityRangeList(fpIn))
+        {
+            fclose(fpIn);
+            BitFileClose(bfpOut);
+            fprintf(stderr, "Error determining frequency ranges.\n");
+            return FALSE;
+        }
+
+        rewind(fpIn);
+
+        /* write information required to decode file to encoded file */
+        WriteHeader(bfpOut);
+    }
+    else
+    {
+        /* initialize probability ranges asumming uniform distribution */
+        InitializeAdaptiveProbabilityRangeList();
+    }
+
+    /* initialize coder start with full probability range [0%, 100%) */
+    lower = 0;
+    upper = ~0;                     /* all ones */
+    underflowBits = 0;
+
+    /* encode symbols one at a time */
+    while ((c = fgetc(fpIn)) != EOF)
+    {
+        ApplySymbolRange(c, staticModel);
+        WriteEncodedBits(bfpOut);
+    }
+
+    fclose(fpIn);
+
+    ApplySymbolRange(EOF_CHAR, staticModel);    /* encode an EOF */
+    WriteEncodedBits(bfpOut);
+
+    WriteRemaining(bfpOut);         /* write out least significant bits */
+    BitFileClose(bfpOut);
+
+    return TRUE;
+}
+
+/***************************************************************************
+*   Function   : SymbolCountToProbabilityRanges
+*   Description: This routine converts the ranges array containing only
+*                symbol counts to an array containing the upper and lower
+*                probability ranges for each symbol.
+*   Parameters : None
+*   Effects    : ranges array containing symbol counts in the upper field
+*                for each symbol is converted to a list of upper and lower
+*                probability bounds for each symbol.
+*   Returned   : None
+***************************************************************************/
+void SymbolCountToProbabilityRanges(void)
+{
+    int c;
+
+    ranges[0] = 0;                              /* absolute lower bound is 0 */
+    ranges[UPPER(EOF_CHAR)] = 1;        /* add 1 EOF character */
+    cumulativeProb++;
+
+    /* assign upper and lower probability ranges */
+    for (c = 1; c <= UPPER(EOF_CHAR); c++)
+    {
+        ranges[c] += ranges[c - 1];
+    }
+
+#ifdef BUILD_DEBUG_OUTPUT
+    /* dump list of ranges */
+    for (c = 0; c < UPPER(EOF_CHAR); c++)
+    {
+        printf("%02X\t%d\t%d\n", c, ranges[LOWER(c)], ranges[UPPER(c)]);
+    }
+#endif
+
+    return;
+}
+
+/***************************************************************************
+*   Function   : BuildProbabilityRangeList
+*   Description: This routine reads the input file and builds the global
+*                list of upper and lower probability ranges for each
+*                symbol.
+*   Parameters : fpIn - file to build range list for
+*   Effects    : ranges array is made to contain probability ranges for
+*                each symbol.
+*   Returned   : TRUE for success, otherwise FALSE.
+***************************************************************************/
+int BuildProbabilityRangeList(FILE *fpIn)
+{
+    int c;
+
+    /***********************************************************************
+    * unsigned long is used to hold the largest counts we can have without
+    * any rescaling.  Rescaling may take place before probability ranges
+    * are computed.
+    ***********************************************************************/
+    unsigned long countArray[EOF_CHAR];
+    unsigned long totalCount = 0;
+    unsigned long rescaleValue;
+
+    if (fpIn == NULL)
+    {
+        return FALSE;
+    }
+
+        /* start with no symbols counted */
+    for (c = 0; c < EOF_CHAR; c++)
+    {
+        countArray[c] = 0;
+    }
+
+    while ((c = fgetc(fpIn)) != EOF)
+    {
+        if (totalCount == ULONG_MAX)
+        {
+            fprintf(stderr, "Error: file too large\n");
+            return FALSE;
+        }
+
+        countArray[c]++;
+        totalCount++;
+    }
+
+    /* rescale counts to be < MAX_PROBABILITY */
+    if (totalCount >= MAX_PROBABILITY)
+    {
+        rescaleValue = (totalCount / MAX_PROBABILITY) + 1;
+
+        for (c = 0; c < EOF_CHAR; c++)
+        {
+            if (countArray[c] > rescaleValue)
+            {
+                countArray[c] /= rescaleValue;
+            }
+            else if (countArray[c] != 0)
+            {
+                countArray[c] = 1;
+            }
+        }
+    }
+
+    /* copy scaled symbol counts to range list */
+    ranges[0] = 0;
+    cumulativeProb = 0;
+    for (c = 0; c < EOF_CHAR; c++)
+    {
+        ranges[UPPER(c)] = countArray[c];
+        cumulativeProb += countArray[c];
+    }
+
+    /* convert counts to a range of probabilities */
+    SymbolCountToProbabilityRanges();
+    return TRUE;
+}
+
+/***************************************************************************
+*   Function   : WriteHeader
+*   Description: This function writes each symbol contained in the encoded
+*                file as well as its rescaled number of occurrences.  A
+*                decoding algorithm may use these numbers to reconstruct
+*                the probability range list used to encode the file.
+*   Parameters : bfpOut - pointer to open binary file to write to.
+*   Effects    : Symbol values and symbol counts are written to a file.
+*   Returned   : None
+***************************************************************************/
+void WriteHeader(bit_file_t *bfpOut)
+{
+    int c;
+    probability_t previous = 0;         /* symbol count so far */
+
+#if BUILD_DEBUG_OUTPUT
+    printf("HEADER:\n");
+#endif
+
+    for(c = 0; c <= (EOF_CHAR - 1); c++)
+    {
+        if (ranges[UPPER(c)] > previous)
+        {
+            /* some of these symbols will be encoded */
+            BitFilePutChar((char)c, bfpOut);
+            previous = (ranges[UPPER(c)] - previous);   /* symbol count */
+
+#if BUILD_DEBUG_OUTPUT
+            printf("%02X\t%d\n", c, previous);
+#endif
+
+            /* write out PRECISION - 2 bit count */
+            BitFilePutBitsInt(bfpOut, &previous, (PRECISION - 2),
+                sizeof(probability_t));
+
+            /* current upper range is previous for the next character */
+            previous = ranges[UPPER(c)];
+        }
+    }
+
+    /* now write end of table (zero count) */
+    BitFilePutChar(0x00, bfpOut);
+    previous = 0;
+    BitFilePutBits(bfpOut, (void *)&previous, PRECISION - 2);
+}
+
+/***************************************************************************
+*   Function   : InitializeAdaptiveProbabilityRangeList
+*   Description: This routine builds the initial global list of upper and
+*                lower probability ranges for each symbol.  This routine
+*                is used by both adaptive encoding and decoding.
+*                Currently it provides a uniform symbol distribution.
+*                Other distributions might be better suited for known data
+*                types (such as English text).
+*   Parameters : NONE
+*   Effects    : ranges array is made to contain initial probability ranges
+*                for each symbol.
+*   Returned   : NONE
+***************************************************************************/
+void InitializeAdaptiveProbabilityRangeList(void)
+{
+    int c;
+
+    cumulativeProb = 0;
+    ranges[0] = 0;          /* absolute lower range */
+
+    /* assign upper and lower probability ranges assuming */
+    for (c = 1; c <= UPPER(EOF_CHAR); c++)
+    {
+        ranges[c] = ranges[c - 1] + 1;
+        cumulativeProb++;
+    }
+
+#ifdef BUILD_DEBUG_OUTPUT
+    /* dump list of ranges */
+    for (c = 0; c < UPPER(EOF_CHAR); c++)
+    {
+        printf("%02X\t%d\t%d\n", c, ranges[LOWER(c)], ranges[UPPER(c)]);
+    }
+#endif
+
+    return;
+}
+
+/***************************************************************************
+*   Function   : ApplySymbolRange
+*   Description: This function is used for both encoding and decoding.  It
+*                applies the range restrictions of a new symbol to the
+*                current upper and lower range bounds of an encoded stream.
+*                If an adaptive model is being used, the probability range
+*                list will be updated after the effect of the symbol is
+*                applied.
+*   Parameters : symbol - The symbol to be added to the current code range
+*                staticModel - TRUE if encoding/decoding with a static
+*                              model.
+*   Effects    : The current upper and lower range bounds are adjusted to
+*                include the range effects of adding another symbol to the
+*                encoded stream.  If an adaptive model is being used, the
+*                probability range list will be updated.
+*   Returned   : None
+***************************************************************************/
+void ApplySymbolRange(int symbol, char staticModel)
+{
+    unsigned long range;        /* must be able to hold max upper + 1 */
+    unsigned long rescaled;     /* range rescaled for range of new symbol */
+                                /* must hold range * max upper */
+
+    /* for updating dynamic models */
+    int i;
+    probability_t original;     /* range value prior to rescale */
+    probability_t delta;        /* range for individual symbol */
+
+    /***********************************************************************
+    * Calculate new upper and lower ranges.  Since the new upper range is
+    * dependant of the old lower range, compute the upper range first.
+    ***********************************************************************/
+    range = (unsigned long)(upper - lower) + 1;         /* current range */
+
+    /* scale upper range of the symbol being coded */
+    rescaled = (unsigned long)ranges[UPPER(symbol)] * range;
+    rescaled /= (unsigned long)cumulativeProb;
+
+    /* new upper = old lower + rescaled new upper - 1*/
+    upper = lower + (probability_t)rescaled - 1;
+
+    /* scale lower range of the symbol being coded */
+    rescaled = (unsigned long)ranges[LOWER(symbol)] * range;
+    rescaled /= (unsigned long)cumulativeProb;
+
+    /* new lower = old lower + rescaled new upper */
+    lower = lower + (probability_t)rescaled;
+
+    if (!staticModel)
+    {
+        /* add new symbol to model */
+        cumulativeProb++;
+        for (i = UPPER(symbol); i <= UPPER(EOF_CHAR); i++)
+        {
+            ranges[i] += 1;
+        }
+
+        /* half current values if cumulativeProb is too large */
+        if (cumulativeProb >= MAX_PROBABILITY)
+        {
+            cumulativeProb = 0;
+            original = 0;
+
+            for (i = 1; i <= UPPER(EOF_CHAR); i++)
+            {
+                delta = ranges[i] - original;
+                if (delta <= 2)
+                {
+                    /* prevent probability from being 0 */
+                    original = ranges[i];
+                    ranges[i] = ranges[i - 1] + 1;
+                }
+                else
+                {
+                    original = ranges[i];
+                    ranges[i] = ranges[i - 1] + (delta / 2);
+                }
+
+                cumulativeProb += (ranges[i] - ranges[i - 1]);
+            }
+        }
+    }
+
+#ifdef BUILD_DEBUG_OUTPUT
+    if (lower > upper)
+    {
+        /* compile this in when testing new models. */
+        fprintf(stderr, "Panic: lower (%X)> upper (%X)\n", lower, upper);
+    }
+#endif
+}
+
+/***************************************************************************
+*   Function   : WriteEncodedBits
+*   Description: This function attempts to shift out as many code bits as
+*                possible, writing the shifted bits to the encoded output
+*                file.  Only bits that will be unchanged when additional
+*                symbols are encoded may be written out.
+*
+*                If the n most significant bits of the lower and upper range
+*                bounds match, they will not be changed when additional
+*                symbols are encoded, so they may be shifted out.
+*
+*                Adjustments are also made to prevent possible underflows
+*                that occur when the upper and lower ranges are so close
+*                that encoding another symbol won't change their values.
+*   Parameters : bfpOut - pointer to open binary file to write to.
+*   Effects    : The upper and lower code bounds are adjusted so that they
+*                only contain only bits that may be affected by the
+*                addition of a new symbol to the encoded stream.
+*   Returned   : None
+***************************************************************************/
+void WriteEncodedBits(bit_file_t *bfpOut)
+{
+    for (;;)
+    {
+        if ((upper & MASK_BIT(0)) == (lower & MASK_BIT(0)))
+        {
+            /* MSBs match, write them to output file */
+            BitFilePutBit((upper & MASK_BIT(0)) != 0, bfpOut);
+
+            /* we can write out underflow bits too */
+            while (underflowBits > 0)
+            {
+                BitFilePutBit((upper & MASK_BIT(0)) == 0, bfpOut);
+                underflowBits--;
+            }
+        }
+        else if ((lower & MASK_BIT(1)) && !(upper & MASK_BIT(1)))
+        {
+            /***************************************************************
+            * Possible underflow condition: neither MSBs nor second MSBs
+            * match.  It must be the case that lower and upper have MSBs of
+            * 01 and 10.  Remove 2nd MSB from lower and upper.
+            ***************************************************************/
+            underflowBits += 1;
+            lower &= ~(MASK_BIT(0) | MASK_BIT(1));
+            upper |= MASK_BIT(1);
+
+            /***************************************************************
+            * The shifts below make the rest of the bit removal work.  If
+            * you don't believe me try it yourself.
+            ***************************************************************/
+        }
+        else
+        {
+            /* nothing left to do */
+            return ;
+        }
+
+        /*******************************************************************
+        * Shift out old MSB and shift in new LSB.  Remember that lower has
+        * all 0s beyond it's end and upper has all 1s beyond it's end.
+        *******************************************************************/
+        lower <<= 1;
+        upper <<= 1;
+        upper |= 1;
+    }
+}
+
+/***************************************************************************
+*   Function   : WriteRemaining
+*   Description: This function writes out all remaining significant bits
+*                in the upper and lower ranges and the underflow bits once
+*                the last symbol has been encoded.
+*   Parameters : bfpOut - pointer to open binary file to write to.
+*   Effects    : Remaining significant range bits are written to the output
+*                file.
+*   Returned   : None
+***************************************************************************/
+void WriteRemaining(bit_file_t *bfpOut)
+{
+    BitFilePutBit((lower & MASK_BIT(1)) != 0, bfpOut);
+
+    /* write out any unwritten underflow bits */
+    for (underflowBits++; underflowBits > 0; underflowBits--)
+    {
+        BitFilePutBit((lower & MASK_BIT(1)) == 0, bfpOut);
+    }
+}
+
+/***************************************************************************
+*   Function   : ArDecodeFile
+*   Description: This routine opens an arithmetically encoded file, reads
+*                it's header, and builds a list of probability ranges which
+*                it then uses to decode the rest of the file.
+*   Parameters : inFile - Name of file to decode
+*                outFile - Name of file to write decoded output to
+*                staticModel - TRUE if decoding with a static model
+*   Effects    : Encoded file is decoded
+*   Returned   : TRUE for success, otherwise FALSE.
+***************************************************************************/
+int ArDecodeFile(char *inFile, char *outFile, char staticModel)
+{
+    int c;
+    probability_t unscaled;
+    bit_file_t *bfpIn;
+    FILE *fpOut;
+
+    /* open input and output files */
+    if ((bfpIn = BitFileOpen(inFile, BF_READ)) == NULL)
+    {
+        perror(inFile);
+        return FALSE;
+    }
+
+    if (outFile == NULL)
+    {
+        fpOut = stdout;
+    }
+    else
+    {
+        if ((fpOut = fopen(outFile, "wb")) == NULL)
+        {
+            BitFileClose(bfpIn);
+            perror(outFile);
+            return FALSE;
+        }
+    }
+
+    if (staticModel)
+    {
+        /* build probability ranges from header in file */
+        if (ReadHeader(bfpIn) == FALSE)
+        {
+            BitFileClose(bfpIn);
+            fclose(fpOut);
+            return FALSE;
+        }
+    }
+
+    /* read start of code and initialize bounds, and adaptive ranges */
+    InitializeDecoder(bfpIn, staticModel);
+
+    /* decode one symbol at a time */
+    for (;;)
+    {
+        /* get the unscaled probability of the current symbol */
+        unscaled = GetUnscaledCode();
+
+        /* figure out which symbol has the above probability */
+        if((c = GetSymbolFromProbability(unscaled)) == -1)
+        {
+            /* error: unknown symbol */
+            break;
+        }
+
+        if (c == EOF_CHAR)
+        {
+            /* no more symbols */
+            break;
+        }
+
+        fputc((char)c, fpOut);
+
+        /* factor out symbol */
+        ApplySymbolRange(c, staticModel);
+        ReadEncodedBits(bfpIn);
+    }
+
+    fclose(fpOut);
+    BitFileClose(bfpIn);
+
+    return TRUE;
+}
+
+/****************************************************************************
+*   Function   : ReadHeader
+*   Description: This function reads the header information stored by
+*                WriteHeader.  The header can then be used to build a
+*                probability range list matching the list that was used to
+*                encode the file.
+*   Parameters : bfpIn - file to read from
+*   Effects    : Probability range list is built.
+*   Returned   : TRUE for success, otherwise FALSE
+****************************************************************************/
+int ReadHeader(bit_file_t *bfpIn)
+{
+    int c;
+    probability_t count;
+
+#if BUILD_DEBUG_OUTPUT
+    printf("HEADER:\n");
+#endif
+
+    cumulativeProb = 0;
+
+    for (c = 0; c <= UPPER(EOF_CHAR); c++)
+    {
+        ranges[UPPER(c)] = 0;
+    }
+
+    /* read [character, probability] sets */
+    for (;;)
+    {
+        c = BitFileGetChar(bfpIn);
+        count = 0;
+
+        /* read (PRECISION - 2) bit count */
+       if (BitFileGetBitsInt(bfpIn, &count, (PRECISION - 2),
+            sizeof(probability_t)) == EOF)
+        {
+            /* premature EOF */
+            fprintf(stderr, "Error: unexpected EOF\n");
+            return FALSE;
+        }
+
+#if BUILD_DEBUG_OUTPUT
+        printf("%02X\t%d\n", c, count);
+#endif
+
+        if (count == 0)
+        {
+            /* 0 count means end of header */
+            break;
+        }
+
+        ranges[UPPER(c)] = count;
+        cumulativeProb += count;
+    }
+
+    /* convert counts to range list */
+    SymbolCountToProbabilityRanges();
+    return TRUE;
+}
+
+/****************************************************************************
+*   Function   : InitializeDecoder
+*   Description: This function starts the upper and lower ranges at their
+*                max/min values and reads in the most significant encoded
+*                bits.
+*   Parameters : bfpIn - file to read from
+*                staticModel - TRUE if decoding using a staticModel
+*   Effects    : upper, lower, and code are initialized.  The probability
+*                range list will also be initialized if an adaptive model
+*                will be used.
+*   Returned   : TRUE for success, otherwise FALSE
+****************************************************************************/
+void InitializeDecoder(bit_file_t *bfpIn, char staticModel)
+{
+    int i;
+
+    if (!staticModel)
+    {
+        /* initialize ranges for adaptive model */
+        InitializeAdaptiveProbabilityRangeList();
+    }
+
+    code = 0;
+
+    /* read PERCISION MSBs of code one bit at a time */
+    for (i = 0; i < PRECISION; i++)
+    {
+        code <<= 1;
+
+        /* treat EOF like 0 */
+        if(BitFileGetBit(bfpIn) == 1)
+        {
+            code |= 1;
+        }
+    }
+
+    /* start with full probability range [0%, 100%) */
+    lower = 0;
+    upper = ~0;         /* all ones */
+}
+
+/****************************************************************************
+*   Function   : GetUnscaledCode
+*   Description: This function undoes the scaling that ApplySymbolRange
+*                performed before bits were shifted out.  The value returned
+*                is the probability of the encoded symbol.
+*   Parameters : None
+*   Effects    : None
+*   Returned   : The probability of the current symbol
+****************************************************************************/
+probability_t GetUnscaledCode(void)
+{
+    unsigned long range;        /* must be able to hold max upper + 1 */
+    unsigned long unscaled;
+
+    range = (unsigned long)(upper - lower) + 1;
+
+    /* reverse the scaling operations from ApplySymbolRange */
+    unscaled = (unsigned long)(code - lower) + 1;
+    unscaled = unscaled * (unsigned long)cumulativeProb - 1;
+    unscaled /= range;
+
+    return ((probability_t)unscaled);
+}
+
+/****************************************************************************
+*   Function   : GetSymbolFromProbability
+*   Description: Given a probability, this function will return the symbol
+*                whose range includes that probability.  Symbol is found
+*                binary search on probability ranges.
+*   Parameters : probability - probability of symbol.
+*   Effects    : None
+*   Returned   : -1 for failure, otherwise encoded symbol
+****************************************************************************/
+int GetSymbolFromProbability(probability_t probability)
+{
+    int first, last, middle;    /* indicies for binary search */
+
+    first = 0;
+    last = UPPER(EOF_CHAR);
+    middle = last / 2;
+
+    /* binary search */
+    while (last >= first)
+    {
+        if (probability < ranges[LOWER(middle)])
+        {
+            /* lower bound is higher than probability */
+            last = middle - 1;
+            middle = first + ((last - first) / 2);
+            continue;
+        }
+
+        if (probability >= ranges[UPPER(middle)])
+        {
+            /* upper bound is lower than probability */
+            first = middle + 1;
+            middle = first + ((last - first) / 2);
+            continue;
+        }
+
+        /* we must have found the right value */
+        return middle;
+    }
+
+    /* error: none of the ranges include the probability */
+    fprintf(stderr, "Unknown Symbol: %d (max: %d)\n", probability,
+        ranges[UPPER(EOF_CHAR)]);
+    return -1;
+}
+
+/***************************************************************************
+*   Function   : ReadEncodedBits
+*   Description: This function attempts to shift out as many code bits as
+*                possible, as bits are shifted out the coded input is
+*                populated with bits from the encoded file.  Only bits
+*                that will be unchanged when additional symbols are decoded
+*                may be shifted out.
+*
+*                If the n most significant bits of the lower and upper range
+*                bounds match, they will not be changed when additional
+*                symbols are decoded, so they may be shifted out.
+*
+*                Adjustments are also made to prevent possible underflows
+*                that occur when the upper and lower ranges are so close
+*                that decoding another symbol won't change their values.
+*   Parameters : bfpOut - pointer to open binary file to read from.
+*   Effects    : The upper and lower code bounds are adjusted so that they
+*                only contain only bits that will be affected by the
+*                addition of a new symbol.  Replacements are read from the
+*                encoded stream.
+*   Returned   : None
+***************************************************************************/
+void ReadEncodedBits(bit_file_t *bfpIn)
+{
+    int nextBit;        /* next bit from encoded input */
+
+    for (;;)
+    {
+        if (( upper & MASK_BIT(0)) == (lower & MASK_BIT(0)))
+        {
+                        /* MSBs match, allow them to be shifted out*/
+        }
+        else if ((lower & MASK_BIT(1)) && !(upper & MASK_BIT(1)))
+        {
+            /***************************************************************
+            * Possible underflow condition: neither MSBs nor second MSBs
+            * match.  It must be the case that lower and upper have MSBs of
+            * 01 and 10.  Remove 2nd MSB from lower and upper.
+            ***************************************************************/
+                        lower   &= ~(MASK_BIT(0) | MASK_BIT(1));
+            upper  |= MASK_BIT(1);
+            code ^= MASK_BIT(1);
+
+            /* the shifts below make the rest of the bit removal work */
+        }
+        else
+        {
+            /* nothing to shift out */
+            return;
+        }
+
+        /*******************************************************************
+        * Shift out old MSB and shift in new LSB.  Remember that lower has
+        * all 0s beyond it's end and upper has all 1s beyond it's end.
+        *******************************************************************/
+                lower <<= 1;
+        upper <<= 1;
+        upper |= 1;
+        code <<= 1;
+
+        if ((nextBit = BitFileGetBit(bfpIn)) == EOF)
+        {
+            /* either all bits are shifted out or error occurred */
+        }
+        else
+        {
+            code |= nextBit;                /* add next encoded bit to code */
+        }
+    }
+
+    return;
+}
diff --git a/test/compression/arcode.h b/test/compression/arcode.h
new file mode 100755
index 000000000..aac320802
--- /dev/null
+++ b/test/compression/arcode.h
@@ -0,0 +1,69 @@
+/***************************************************************************
+*            Header for Arithmetic Encoding and Decoding Library
+*
+*   File    : arcode.h
+*   Purpose : Provides prototypes for functions that use arithmetic coding
+*             to encode/decode files.
+*   Author  : Michael Dipperstein
+*   Date    : April 2, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: arcode.h,v 1.3 2007/09/08 15:47:02 michael Exp $
+*   $Log: arcode.h,v $
+*   Revision 1.3  2007/09/08 15:47:02  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.2  2004/08/13 13:09:46  michael
+*   Add support for adaptive encoding
+*
+*   Revision 1.1.1.1  2004/04/04 14:54:13  michael
+*   Initial version
+*
+****************************************************************************
+*
+* Arcode: An ANSI C Arithmetic Encoding/Decoding Routines
+* Copyright (C) 2004, 2006-2007 by Michael Dipperstein (mdipper@cs.ucsb.edu)
+*
+* This file is part of the arcode library.
+*
+* The arcode library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The arcode library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+#ifndef _ARCODE_H_
+#define _ARCODE_H_
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#ifndef FALSE
+#define FALSE       0
+#endif
+
+#ifndef TRUE
+#define TRUE        1
+#endif
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+ /* encode inFile */
+int ArEncodeFile(char *inFile, char *outFile, char staticModel);
+
+/* decode inFile*/
+int ArDecodeFile(char *inFile, char *outFile, char staticModel);
+
+#endif  /* ndef _ARCODE_H_ */
diff --git a/test/compression/armain.c b/test/compression/armain.c
new file mode 100755
index 000000000..8f37c4ff9
--- /dev/null
+++ b/test/compression/armain.c
@@ -0,0 +1,266 @@
+/***************************************************************************
+*              Sample Program Using Arithmetic Encoding Library
+*
+*   File    : sample.c
+*   Purpose : Demonstrate usage of arithmetic encoding library
+*   Author  : Michael Dipperstein
+*   Date    : March 10, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: sample.c,v 1.3 2007/09/08 15:48:14 michael Exp $
+*   $Log: sample.c,v $
+*   Revision 1.3  2007/09/08 15:48:14  michael
+*   Replace getopt with optlist.
+*   Changes required for LGPL v3.
+*
+*   Revision 1.2  2004/08/13 13:08:43  michael
+*   Add support for adaptive encoding
+*
+*   Use executable name in help messages
+*
+*   Revision 1.1.1.1  2004/04/04 14:54:13  michael
+*   Initial version
+*
+*
+****************************************************************************
+*
+* SAMPLE: Sample usage of the arcode Arithmetic Encoding Library
+* Copyright (C) 2004, 2007 by Michael Dipperstein (mdipper@cs.ucsb.edu)
+*
+* This file is part of the arcode library.
+*
+* The arcode library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The arcode library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "optlist.h"
+#include "arcode.h"
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+char *RemovePath(char *fullPath);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : main
+*   Description: This is the main function for this program, it validates
+*                the command line input and, if valid, it will call
+*                functions to encode or decode a file using the arithmetic
+*                coding algorithm.
+*   Parameters : argc - number of parameters
+*                argv - parameter list
+*   Effects    : Encodes/Decodes input file
+*   Returned   : EXIT_SUCCESS for success, otherwise EXIT_FAILURE.
+****************************************************************************/
+int main(int argc, char *argv[])
+{
+    option_t *optList, *thisOpt;
+    char *inFile, *outFile; /* name of input & output files */
+    char encode;            /* encode/decode */
+    char staticModel;       /* static/adaptive model*/
+
+    /* initialize data */
+    inFile = NULL;
+    outFile = NULL;
+    encode = TRUE;
+    staticModel = TRUE;
+
+    /* parse command line */
+    optList = GetOptList(argc, argv, "acdi:o:h?");
+    thisOpt = optList;
+
+    while (thisOpt != NULL)
+    {
+        switch(thisOpt->option)
+        {
+            case 'a':       /* adaptive model vs. static */
+                staticModel = FALSE;
+                break;
+
+            case 'c':       /* compression mode */
+                encode = TRUE;
+                break;
+
+            case 'd':       /* decompression mode */
+                encode = FALSE;
+                break;
+
+            case 'i':       /* input file name */
+                if (inFile != NULL)
+                {
+                    fprintf(stderr, "Multiple input files not allowed.\n");
+                    free(inFile);
+
+                    if (outFile != NULL)
+                    {
+                        free(outFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                else if ((inFile =
+                    (char *)malloc(strlen(thisOpt->argument) + 1)) == NULL)
+                {
+                    perror("Memory allocation");
+
+                    if (outFile != NULL)
+                    {
+                        free(outFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                strcpy(inFile, thisOpt->argument);
+                break;
+
+            case 'o':       /* output file name */
+                if (outFile != NULL)
+                {
+                    fprintf(stderr, "Multiple output files not allowed.\n");
+                    free(outFile);
+
+                    if (inFile != NULL)
+                    {
+                        free(inFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                else if ((outFile =
+                    (char *)malloc(strlen(thisOpt->argument) + 1)) == NULL)
+                {
+                    perror("Memory allocation");
+
+                    if (inFile != NULL)
+                    {
+                        free(inFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                strcpy(outFile, thisOpt->argument);
+                break;
+
+            case 'h':
+            case '?':
+                printf("Usage: %s <options>\n\n", RemovePath(argv[0]));
+                printf("options:\n");
+                printf("  -c : Encode input file to output file.\n");
+                printf("  -d : Decode input file to output file.\n");
+                printf("  -i <filename> : Name of input file.\n");
+                printf("  -o <filename> : Name of output file.\n");
+                printf("  -a : Use adaptive model instead of static.\n");
+                printf("  -h | ?  : Print out command line options.\n\n");
+                printf("Default: %s -c\n", RemovePath(argv[0]));
+
+                FreeOptList(optList);
+                return(EXIT_SUCCESS);
+        }
+
+        optList = thisOpt->next;
+        free(thisOpt);
+        thisOpt = optList;
+    }
+
+    /* validate command line */
+    if (inFile == NULL)
+    {
+        fprintf(stderr, "Input file must be provided\n");
+        fprintf(stderr, "Enter \"%s -?\" for help.\n", RemovePath(argv[0]));
+
+        if (outFile != NULL)
+        {
+            free(outFile);
+        }
+
+        exit (EXIT_FAILURE);
+    }
+    else if (outFile == NULL)
+    {
+        fprintf(stderr, "Output file must be provided\n");
+        fprintf(stderr, "Enter \"%s -?\" for help.\n", RemovePath(argv[0]));
+
+        if (inFile != NULL)
+        {
+            free(inFile);
+        }
+
+        exit (EXIT_FAILURE);
+    }
+
+    /* we have valid parameters encode or decode */
+    if (encode)
+    {
+        ArEncodeFile(inFile, outFile, staticModel);
+    }
+    else
+    {
+        ArDecodeFile(inFile, outFile, staticModel);
+    }
+
+    free(inFile);
+    free(outFile);
+    return EXIT_SUCCESS;
+}
+
+/****************************************************************************
+*   Function   : RemovePath
+*   Description: This is function accepts a pointer to the name of a file
+*                along with path information and returns a pointer to the
+*                character that is not part of the path.
+*   Parameters : fullPath - pointer to an array of characters containing
+*                           a file name and possible path modifiers.
+*   Effects    : None
+*   Returned   : Returns a pointer to the first character after any path
+*                information.
+****************************************************************************/
+char *RemovePath(char *fullPath)
+{
+    int i;
+    char *start, *tmp;                          /* start of file name */
+    const char delim[3] = {'\\', '/', ':'};     /* path deliminators */
+
+    start = fullPath;
+
+    /* find the first character after all file path delimiters */
+    for (i = 0; i < 3; i++)
+    {
+        tmp = strrchr(start, delim[i]);
+
+        if (tmp != NULL)
+        {
+            start = tmp + 1;
+        }
+    }
+
+    return start;
+}
diff --git a/test/compression/bitfile.c b/test/compression/bitfile.c
new file mode 100644
index 000000000..7480ce999
--- /dev/null
+++ b/test/compression/bitfile.c
@@ -0,0 +1,1043 @@
+/***************************************************************************
+*                        Bit Stream File Implementation
+*
+*   File    : bitfile.c
+*   Purpose : This file implements a simple library of I/O functions for
+*             files that contain data in sizes that aren't integral bytes.
+*             An attempt was made to make the functions in this library
+*             analogous to functions provided to manipulate byte streams.
+*             The functions contained in this library were created with
+*             compression algorithms in mind, but may be suited to other
+*             applications.
+*   Author  : Michael Dipperstein
+*   Date    : January 9, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: bitfile.c,v 1.10 2007/08/26 21:53:48 michael Exp $
+*   $Log: bitfile.c,v $
+*   Revision 1.10  2007/08/26 21:53:48  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.9  2007/07/10 05:34:07  michael
+*   Remove ',' after last element in the enum endian_t.
+*
+*   Revision 1.8  2007/02/06 06:22:07  michael
+*   Used trim program to remove trailing spaces.
+*
+*   Revision 1.7  2006/06/03 19:32:38  michael
+*   Corrected error reporetd anonymous.  The allocation of constants used to
+*   open underlying read/write/append files did not account for a terminating
+*   null.
+*
+*   Used spell checker to correct spelling.
+*
+*   Revision 1.6  2005/12/08 06:56:55  michael
+*   Minor text corrections.
+*
+*   Revision 1.5  2005/12/06 15:06:37  michael
+*   Added BitFileGetBitsInt and BitFilePutBitsInt for integer types.
+*
+*   Revision 1.4  2005/06/23 04:34:18  michael
+*   Prevent BitfileGetBits/PutBits from accessing an extra byte when given
+*   an integral number of bytes.
+*
+*   Revision 1.3  2004/11/09 14:16:58  michael
+*   Added functions to convert open bit_file_t to FILE and to
+*   align open bit_file_t to the next byte.
+*
+*   Revision 1.2  2004/06/15 13:15:58  michael
+*   Use incomplete type to hide definition of bitfile structure
+*
+*   Revision 1.1.1.1  2004/02/09 05:31:42  michael
+*   Initial release
+*
+*
+****************************************************************************
+*
+* Bitfile: Bit stream File I/O Routines
+* Copyright (C) 2004-2007 by Michael Dipperstein (mdipper@cs.ucsb.edu)
+*
+* This file is part of the bit file library.
+*
+* The bit file library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The bit file library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdlib.h>
+#include <errno.h>
+#include "bitfile.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+typedef enum
+{
+    BF_UNKNOWN_ENDIAN,
+    BF_LITTLE_ENDIAN,
+    BF_BIG_ENDIAN
+} endian_t;
+
+struct bit_file_t
+{
+    FILE *fp;                   /* file pointer used by stdio functions */
+    endian_t endian;            /* endianess of architecture */
+    unsigned char bitBuffer;    /* bits waiting to be read/written */
+    unsigned char bitCount;     /* number of bits in bitBuffer */
+    BF_MODES mode;              /* open for read, write, or append */
+};
+
+/* union used to test for endianess */
+typedef union
+{
+    unsigned long word;
+    unsigned char bytes[sizeof(unsigned long)];
+} endian_test_t;
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+endian_t DetermineEndianess(void);
+
+int BitFilePutBitsLE(bit_file_t *stream, void *bits, const unsigned int count);
+int BitFilePutBitsBE(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size);
+
+int BitFileGetBitsLE(bit_file_t *stream, void *bits, const unsigned int count);
+int BitFileGetBitsBE(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/***************************************************************************
+*   Function   : BitFileOpen
+*   Description: This function opens a bit file for reading, writing,
+*                or appending.  If successful, a bit_file_t data
+*                structure will be allocated and a pointer to the
+*                structure will be returned.
+*   Parameters : fileName - NULL terminated string containing the name of
+*                           the file to be opened.
+*                mode - The mode of the file to be opened
+*   Effects    : The specified file will be opened and file structure will
+*                be allocated.
+*   Returned   : Pointer to the bit_file_t structure for the bit file
+*                opened, or NULL on failure.  errno will be set for all
+*                failure cases.
+***************************************************************************/
+bit_file_t *BitFileOpen(const char *fileName, const BF_MODES mode)
+{
+    char modes[3][3] = {"rb", "wb", "ab"};  /* binary modes for fopen */
+    bit_file_t *bf;
+
+    bf = (bit_file_t *)malloc(sizeof(bit_file_t));
+
+    if (bf == NULL)
+    {
+        /* malloc failed */
+        errno = ENOMEM;
+    }
+    else
+    {
+        bf->fp = fopen(fileName, modes[mode]);
+
+        if (bf->fp == NULL)
+        {
+            /* fopen failed */
+            free(bf);
+            bf = NULL;
+        }
+        else
+        {
+            /* fopen succeeded fill in remaining bf data */
+            bf->bitBuffer = 0;
+            bf->bitCount = 0;
+            bf->mode = mode;
+
+            /***************************************************************
+            *  TO DO: Consider using the last byte in a file to indicate
+            * the number of bits in the previous byte that actually have
+            * data.  If I do that, I'll need special handling of files
+            * opened with a mode of BF_APPEND.
+            ***************************************************************/
+        }
+    }
+
+    bf->endian = DetermineEndianess();
+    return (bf);
+}
+
+/***************************************************************************
+*   Function   : MakeBitFile
+*   Description: This function naively wraps a standard file in a
+*                bit_file_t structure.  ANSI-C doesn't support file status
+*                functions commonly found in other C variants, so the
+*                caller must be passed as a parameter.
+*   Parameters : stream - pointer to the standard file being wrapped.
+*                mode - The mode of the file being wrapped.
+*   Effects    : A bit_file_t structure will be created for the stream
+*                passed as a parameter.
+*   Returned   : Pointer to the bit_file_t structure for the bit file
+*                or NULL on failure.  errno will be set for all failure
+*                cases.
+***************************************************************************/
+bit_file_t *MakeBitFile(FILE *stream, const BF_MODES mode)
+{
+    bit_file_t *bf;
+
+    if (stream == NULL)
+    {
+        /* can't wrapper empty steam */
+        errno = EBADF;
+        bf = NULL;
+    }
+    else
+    {
+        bf = (bit_file_t *)malloc(sizeof(bit_file_t));
+
+        if (bf == NULL)
+        {
+            /* malloc failed */
+            errno = ENOMEM;
+        }
+        else
+        {
+            /* set structure data */
+            bf->fp = stream;
+            bf->bitBuffer = 0;
+            bf->bitCount = 0;
+            bf->mode = mode;
+        }
+    }
+
+    bf->endian = DetermineEndianess();
+
+    return (bf);
+}
+
+/***************************************************************************
+*   Function   : DetermineEndianess
+*   Description: This function determines the endianess of the current
+*                hardware architecture.  An unsigned long is set to 1.  If
+*                the 1st byte of the unsigned long gets the 1, this is a
+*                little endian machine.  If the last byte gets the 1, this
+*                is a big endian machine.
+*   Parameters : None
+*   Effects    : None
+*   Returned   : endian_t for current machine architecture
+***************************************************************************/
+endian_t DetermineEndianess(void)
+{
+    endian_t endian;
+    endian_test_t endianTest;
+
+    endianTest.word = 1;
+
+    if (endianTest.bytes[0] == 1)
+    {
+        /* LSB is 1st byte (little endian)*/
+        endian = BF_LITTLE_ENDIAN;
+    }
+    else if (endianTest.bytes[sizeof(unsigned long) - 1] == 1)
+    {
+        /* LSB is last byte (big endian)*/
+        endian = BF_BIG_ENDIAN;
+    }
+    else
+    {
+        endian = BF_UNKNOWN_ENDIAN;
+    }
+
+    return endian;
+}
+
+/***************************************************************************
+*   Function   : BitFileClose
+*   Description: This function closes a bit file and frees all associated
+*                data.
+*   Parameters : stream - pointer to bit file stream being closed
+*   Effects    : The specified file will be closed and the file structure
+*                will be freed.
+*   Returned   : 0 for success or EOF for failure.
+***************************************************************************/
+int BitFileClose(bit_file_t *stream)
+{
+    int returnValue = 0;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    if ((stream->mode == BF_WRITE) || (stream->mode == BF_APPEND))
+    {
+        /* write out any unwritten bits */
+        if (stream->bitCount != 0)
+        {
+            (stream->bitBuffer) <<= 8 - (stream->bitCount);
+            fputc(stream->bitBuffer, stream->fp);   /* handle error? */
+        }
+    }
+
+    /***********************************************************************
+    *  TO DO: Consider writing an additional byte indicating the number of
+    *  valid bits (bitCount) in the previous byte.
+    ***********************************************************************/
+
+    /* close file */
+    returnValue = fclose(stream->fp);
+
+    /* free memory allocated for bit file */
+    free(stream);
+
+    return(returnValue);
+}
+
+/***************************************************************************
+*   Function   : BitFileToFILE
+*   Description: This function flushes and frees the bitfile structure,
+*                returning a pointer to a stdio file.
+*   Parameters : stream - pointer to bit file stream being closed
+*   Effects    : The specified bitfile will be made usable as a stdio
+*                FILE.
+*   Returned   : Pointer to FILE.  NULL for failure.
+***************************************************************************/
+FILE *BitFileToFILE(bit_file_t *stream)
+{
+    FILE *fp = NULL;
+
+    if (stream == NULL)
+    {
+        return(NULL);
+    }
+
+    if ((stream->mode == BF_WRITE) || (stream->mode == BF_APPEND))
+    {
+        /* write out any unwritten bits */
+        if (stream->bitCount != 0)
+        {
+            (stream->bitBuffer) <<= 8 - (stream->bitCount);
+            fputc(stream->bitBuffer, stream->fp);   /* handle error? */
+        }
+    }
+
+    /***********************************************************************
+    *  TO DO: Consider writing an additional byte indicating the number of
+    *  valid bits (bitCount) in the previous byte.
+    ***********************************************************************/
+
+    /* close file */
+    fp = stream->fp;
+
+    /* free memory allocated for bit file */
+    free(stream);
+
+    return(fp);
+}
+
+/***************************************************************************
+*   Function   : BitFileByteAlign
+*   Description: This function aligns the bitfile to the nearest byte.  For
+*                output files, this means writing out the bit buffer with
+*                extra bits set to 0.  For input files, this means flushing
+*                the bit buffer.
+*   Parameters : stream - pointer to bit file stream to align
+*   Effects    : Flushes out the bit buffer.
+*   Returned   : EOF if stream is NULL.  Otherwise, the contents of the
+*                bit buffer.
+***************************************************************************/
+int BitFileByteAlign(bit_file_t *stream)
+{
+    int returnValue;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    returnValue = stream->bitBuffer;
+
+    if ((stream->mode == BF_WRITE) || (stream->mode == BF_APPEND))
+    {
+        /* write out any unwritten bits */
+        if (stream->bitCount != 0)
+        {
+            (stream->bitBuffer) <<= 8 - (stream->bitCount);
+            fputc(stream->bitBuffer, stream->fp);   /* handle error? */
+        }
+    }
+
+    stream->bitBuffer = 0;
+    stream->bitCount = 0;
+
+    return (returnValue);
+}
+
+/***************************************************************************
+*   Function   : BitFileGetChar
+*   Description: This function returns the next byte from the file passed as
+*                a parameter.
+*   Parameters : stream - pointer to bit file stream to read from
+*   Effects    : Reads next byte from file and updates buffer accordingly.
+*   Returned   : EOF if a whole byte cannot be obtained.  Otherwise,
+*                the character read.
+***************************************************************************/
+int BitFileGetChar(bit_file_t *stream)
+{
+    int returnValue;
+    unsigned char tmp;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    returnValue = fgetc(stream->fp);
+
+    if (stream->bitCount == 0)
+    {
+        /* we can just get byte from file */
+        return returnValue;
+    }
+
+    /* we have some buffered bits to return too */
+    if (returnValue != EOF)
+    {
+        /* figure out what to return */
+        tmp = ((unsigned char)returnValue) >> (stream->bitCount);
+        tmp |= ((stream->bitBuffer) << (8 - (stream->bitCount)));
+
+        /* put remaining in buffer. count shouldn't change. */
+        stream->bitBuffer = returnValue;
+
+        returnValue = tmp;
+    }
+
+    return returnValue;
+}
+
+/***************************************************************************
+*   Function   : BitFilePutChar
+*   Description: This function writes the byte passed as a parameter to the
+*                file passed a parameter.
+*   Parameters : c - the character to be written
+*                stream - pointer to bit file stream to write to
+*   Effects    : Writes a byte to the file and updates buffer accordingly.
+*   Returned   : On success, the character written, otherwise EOF.
+***************************************************************************/
+int BitFilePutChar(const int c, bit_file_t *stream)
+{
+    unsigned char tmp;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    if (stream->bitCount == 0)
+    {
+        /* we can just put byte from file */
+        return fputc(c, stream->fp);
+    }
+
+    /* figure out what to write */
+    tmp = ((unsigned char)c) >> (stream->bitCount);
+    tmp = tmp | ((stream->bitBuffer) << (8 - stream->bitCount));
+
+    if (fputc(tmp, stream->fp) != EOF)
+    {
+        /* put remaining in buffer. count shouldn't change. */
+        stream->bitBuffer = c;
+    }
+    else
+    {
+        return EOF;
+    }
+
+    return tmp;
+}
+
+/***************************************************************************
+*   Function   : BitFileGetBit
+*   Description: This function returns the next bit from the file passed as
+*                a parameter.  The bit value returned is the msb in the
+*                bit buffer.
+*   Parameters : stream - pointer to bit file stream to read from
+*   Effects    : Reads next bit from bit buffer.  If the buffer is empty,
+*                a new byte will be read from the file.
+*   Returned   : 0 if bit == 0, 1 if bit == 1, and EOF if operation fails.
+***************************************************************************/
+int BitFileGetBit(bit_file_t *stream)
+{
+    int returnValue;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    if (stream->bitCount == 0)
+    {
+        /* buffer is empty, read another character */
+        if ((returnValue = fgetc(stream->fp)) == EOF)
+        {
+            return EOF;
+        }
+        else
+        {
+            stream->bitCount = 8;
+            stream->bitBuffer = returnValue;
+        }
+    }
+
+    /* bit to return is msb in buffer */
+    stream->bitCount--;
+    returnValue = (stream->bitBuffer) >> (stream->bitCount);
+
+    return (returnValue & 0x01);
+}
+
+/***************************************************************************
+*   Function   : BitFilePutBit
+*   Description: This function writes the bit passed as a parameter to the
+*                file passed a parameter.
+*   Parameters : c - the bit value to be written
+*                stream - pointer to bit file stream to write to
+*   Effects    : Writes a bit to the bit buffer.  If the buffer has a byte,
+*                the buffer is written to the file and cleared.
+*   Returned   : On success, the bit value written, otherwise EOF.
+***************************************************************************/
+int BitFilePutBit(const int c, bit_file_t *stream)
+{
+    int returnValue = c;
+
+    if (stream == NULL)
+    {
+        return(EOF);
+    }
+
+    stream->bitCount++;
+    stream->bitBuffer <<= 1;
+
+    if (c != 0)
+    {
+        stream->bitBuffer |= 1;
+    }
+
+    /* write bit buffer if we have 8 bits */
+    if (stream->bitCount == 8)
+    {
+        if (fputc(stream->bitBuffer, stream->fp) == EOF)
+        {
+            returnValue = EOF;
+        }
+
+        /* reset buffer */
+        stream->bitCount = 0;
+        stream->bitBuffer = 0;
+    }
+
+    return returnValue;
+}
+
+/***************************************************************************
+*   Function   : BitFileGetBits
+*   Description: This function reads the specified number of bits from the
+*                file passed as a parameter and writes them to the
+*                requested memory location (msb to lsb).
+*   Parameters : stream - pointer to bit file stream to read from
+*                bits - address to store bits read
+*                count - number of bits to read
+*   Effects    : Reads bits from the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.
+*   Returned   : EOF for failure, otherwise the number of bits read.  If
+*                an EOF is reached before all the bits are read, bits
+*                will contain every bit through the last complete byte.
+***************************************************************************/
+int BitFileGetBits(bit_file_t *stream, void *bits, const unsigned int count)
+{
+    unsigned char *bytes, shifts;
+    int offset, remaining, returnValue;
+
+    bytes = (unsigned char *)bits;
+
+    if ((stream == NULL) || (bits == NULL))
+    {
+        return(EOF);
+    }
+
+    offset = 0;
+    remaining = count;
+
+    /* read whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFileGetChar(stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        bytes[offset] = (unsigned char)returnValue;
+        remaining -= 8;
+        offset++;
+    }
+
+    if (remaining != 0)
+    {
+        /* read remaining bits */
+        shifts = 8 - remaining;
+
+        while (remaining > 0)
+        {
+            returnValue = BitFileGetBit(stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            bytes[offset] <<= 1;
+            bytes[offset] |= (returnValue & 0x01);
+            remaining--;
+        }
+
+        /* shift last bits into position */
+        bytes[offset] <<= shifts;
+    }
+
+    return count;
+}
+
+/***************************************************************************
+*   Function   : BitFilePutBits
+*   Description: This function writes the specified number of bits from the
+*                memory location passed as a parameter to the file passed
+*                as a parameter.   Bits are written msb to lsb.
+*   Parameters : stream - pointer to bit file stream to write to
+*                bits - pointer to bits to write
+*                count - number of bits to write
+*   Effects    : Writes bits to the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.
+*   Returned   : EOF for failure, otherwise the number of bits written.  If
+*                an error occurs after a partial write, the partially
+*                written bits will not be unwritten.
+***************************************************************************/
+int BitFilePutBits(bit_file_t *stream, void *bits, const unsigned int count)
+{
+    unsigned char *bytes, tmp;
+    int offset, remaining, returnValue;
+
+    bytes = (unsigned char *)bits;
+
+    if ((stream == NULL) || (bits == NULL))
+    {
+        return(EOF);
+    }
+
+    offset = 0;
+    remaining = count;
+
+    /* write whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFilePutChar(bytes[offset], stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        remaining -= 8;
+        offset++;
+    }
+
+    if (remaining != 0)
+    {
+        /* write remaining bits */
+        tmp = bytes[offset];
+
+        while (remaining > 0)
+        {
+            returnValue = BitFilePutBit((tmp & 0x80), stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            tmp <<= 1;
+            remaining--;
+        }
+    }
+
+    return count;
+}
+
+/***************************************************************************
+*   Function   : BitFileGetBitsInt
+*   Description: This function provides a machine independent layer that
+*                allows a single function call to stuff an arbitrary number
+*                of bits into an integer type variable.
+*   Parameters : stream - pointer to bit file stream to read from
+*                bits - address to store bits read
+*                count - number of bits to read
+*                size - sizeof type containing "bits"
+*   Effects    : Calls a function that reads bits from the bit buffer and
+*                file stream.  The bit buffer will be modified as necessary.
+*                the bits will be written to "bits" from least significant
+*                byte to most significant byte.
+*   Returned   : EOF for failure, otherwise the number of bits read by the
+*                called function.
+***************************************************************************/
+int BitFileGetBitsInt(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size)
+{
+    int returnValue;
+
+    if ((stream == NULL) || (bits == NULL))
+    {
+        return(EOF);
+    }
+
+    if (stream->endian == BF_LITTLE_ENDIAN)
+    {
+        returnValue = BitFileGetBitsLE(stream, bits, count);
+    }
+    else if (stream->endian == BF_BIG_ENDIAN)
+    {
+        returnValue = BitFileGetBitsBE(stream, bits, count, size);
+    }
+    else
+    {
+        returnValue = EOF;
+    }
+
+    return returnValue;
+}
+
+/***************************************************************************
+*   Function   : BitFileGetBitsLE   (Little Endian)
+*   Description: This function reads the specified number of bits from the
+*                file passed as a parameter and writes them to the
+*                requested memory location (LSB to MSB).
+*   Parameters : stream - pointer to bit file stream to read from
+*                bits - address to store bits read
+*                count - number of bits to read
+*   Effects    : Reads bits from the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.  bits is treated as
+*                a little endian integer of length >= (count/8) + 1.
+*   Returned   : EOF for failure, otherwise the number of bits read.  If
+*                an EOF is reached before all the bits are read, bits
+*                will contain every bit through the last successful read.
+***************************************************************************/
+int BitFileGetBitsLE(bit_file_t *stream, void *bits, const unsigned int count)
+{
+    unsigned char *bytes, shifts;
+    int offset, remaining, returnValue;
+
+    bytes = (unsigned char *)bits;
+
+    offset = 0;
+    remaining = count;
+
+    /* read whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFileGetChar(stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        bytes[offset] = (unsigned char)returnValue;
+        remaining -= 8;
+        offset++;
+    }
+
+    if (remaining != 0)
+    {
+        /* read remaining bits */
+        shifts = 8 - remaining;
+
+        while (remaining > 0)
+        {
+            returnValue = BitFileGetBit(stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            bytes[offset] <<= 1;
+            bytes[offset] |= (returnValue & 0x01);
+            remaining--;
+        }
+
+    }
+
+    return count;
+}
+
+/***************************************************************************
+*   Function   : BitFileGetBitsBE   (Big Endian)
+*   Description: This function reads the specified number of bits from the
+*                file passed as a parameter and writes them to the
+*                requested memory location (LSB to MSB).
+*   Parameters : stream - pointer to bit file stream to read from
+*                bits - address to store bits read
+*                count - number of bits to read
+*                size - sizeof type containing "bits"
+*   Effects    : Reads bits from the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.  bits is treated as
+*                a big endian integer of length size.
+*   Returned   : EOF for failure, otherwise the number of bits read.  If
+*                an EOF is reached before all the bits are read, bits
+*                will contain every bit through the last successful read.
+***************************************************************************/
+int BitFileGetBitsBE(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size)
+{
+    unsigned char *bytes, shifts;
+    int offset, remaining, returnValue;
+
+    if (count > (size * 8))
+    {
+        /* too many bits to read */
+        return EOF;
+    }
+
+    bytes = (unsigned char *)bits;
+
+    offset = size - 1;
+    remaining = count;
+
+    /* read whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFileGetChar(stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        bytes[offset] = (unsigned char)returnValue;
+        remaining -= 8;
+        offset--;
+    }
+
+    if (remaining != 0)
+    {
+        /* read remaining bits */
+        shifts = 8 - remaining;
+
+        while (remaining > 0)
+        {
+            returnValue = BitFileGetBit(stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            bytes[offset] <<= 1;
+            bytes[offset] |= (returnValue & 0x01);
+            remaining--;
+        }
+
+    }
+
+    return count;
+}
+
+/***************************************************************************
+*   Function   : BitFilePutBitsInt
+*   Description: This function provides a machine independent layer that
+*                allows a single function call to write an arbitrary number
+*                of bits from an integer type variable into a file.
+*   Parameters : stream - pointer to bit file stream to write to
+*                bits - pointer to bits to write
+*                count - number of bits to write
+*                size - sizeof type containing "bits"
+*   Effects    : Calls a function that writes bits to the bit buffer and
+*                file stream.  The bit buffer will be modified as necessary.
+*                the bits will be written to the file stream from least
+*                significant byte to most significant byte.
+*   Returned   : EOF for failure, otherwise the number of bits written.  If
+*                an error occurs after a partial write, the partially
+*                written bits will not be unwritten.
+***************************************************************************/
+int BitFilePutBitsInt(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size)
+{
+    int returnValue;
+
+    if ((stream == NULL) || (bits == NULL))
+    {
+        return(EOF);
+    }
+
+    if (stream->endian == BF_LITTLE_ENDIAN)
+    {
+        returnValue = BitFilePutBitsLE(stream, bits, count);
+    }
+    else if (stream->endian == BF_BIG_ENDIAN)
+    {
+        returnValue = BitFilePutBitsBE(stream, bits, count, size);
+    }
+    else
+    {
+        returnValue = EOF;
+    }
+
+    return returnValue;
+}
+
+/***************************************************************************
+*   Function   : BitFilePutBitsLE   (Little Endian)
+*   Description: This function writes the specified number of bits from the
+*                memory location passed as a parameter to the file passed
+*                as a parameter.   Bits are written LSB to MSB.
+*   Parameters : stream - pointer to bit file stream to write to
+*                bits - pointer to bits to write
+*                count - number of bits to write
+*   Effects    : Writes bits to the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.  bits is treated as
+*                a little endian integer of length >= (count/8) + 1.
+*   Returned   : EOF for failure, otherwise the number of bits written.  If
+*                an error occurs after a partial write, the partially
+*                written bits will not be unwritten.
+***************************************************************************/
+int BitFilePutBitsLE(bit_file_t *stream, void *bits, const unsigned int count)
+{
+    unsigned char *bytes, tmp;
+    int offset, remaining, returnValue;
+
+    bytes = (unsigned char *)bits;
+    offset = 0;
+    remaining = count;
+
+    /* write whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFilePutChar(bytes[offset], stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        remaining -= 8;
+        offset++;
+    }
+
+    if (remaining != 0)
+    {
+        /* write remaining bits */
+        tmp = bytes[offset];
+        tmp <<= (8 - remaining);
+
+        while (remaining > 0)
+        {
+            returnValue = BitFilePutBit((tmp & 0x80), stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            tmp <<= 1;
+            remaining--;
+        }
+    }
+
+    return count;
+}
+
+/***************************************************************************
+*   Function   : BitFilePutBitsBE   (Big Endian)
+*   Description: This function writes the specified number of bits from the
+*                memory location passed as a parameter to the file passed
+*                as a parameter.   Bits are written LSB to MSB.
+*   Parameters : stream - pointer to bit file stream to write to
+*                bits - pointer to bits to write
+*                count - number of bits to write
+*   Effects    : Writes bits to the bit buffer and file stream.  The bit
+*                buffer will be modified as necessary.  bits is treated as
+*                a big endian integer of length size.
+*   Returned   : EOF for failure, otherwise the number of bits written.  If
+*                an error occurs after a partial write, the partially
+*                written bits will not be unwritten.
+***************************************************************************/
+int BitFilePutBitsBE(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size)
+{
+    unsigned char *bytes, tmp;
+    int offset, remaining, returnValue;
+
+    if (count > (size * 8))
+    {
+        /* too many bits to write */
+        return EOF;
+    }
+
+    bytes = (unsigned char *)bits;
+    offset = size - 1;
+    remaining = count;
+
+    /* write whole bytes */
+    while (remaining >= 8)
+    {
+        returnValue = BitFilePutChar(bytes[offset], stream);
+
+        if (returnValue == EOF)
+        {
+            return EOF;
+        }
+
+        remaining -= 8;
+        offset--;
+    }
+
+    if (remaining != 0)
+    {
+        /* write remaining bits */
+        tmp = bytes[offset];
+        tmp <<= (8 - remaining);
+
+        while (remaining > 0)
+        {
+            returnValue = BitFilePutBit((tmp & 0x80), stream);
+
+            if (returnValue == EOF)
+            {
+                return EOF;
+            }
+
+            tmp <<= 1;
+            remaining--;
+        }
+    }
+
+    return count;
+}
diff --git a/test/compression/bitfile.h b/test/compression/bitfile.h
new file mode 100644
index 000000000..20d102c7e
--- /dev/null
+++ b/test/compression/bitfile.h
@@ -0,0 +1,121 @@
+/***************************************************************************
+*                            Bit Stream File Header
+*
+*   File    : bitfile.h
+*   Purpose : Provides definitions and prototypes for a simple library of
+*             I/O functions for files that contain data in sizes that aren't
+*             integral bytes.
+*             An attempt was made to make the functions in this library
+*             analogous to functions provided to manipulate byte streams.
+*             The functions contained in this library were created with
+*             compression algorithms in mind, but may be suited to other
+*             applications.
+*   Author  : Michael Dipperstein
+*   Date    : January 9, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: bitfile.h,v 1.6 2007/08/26 21:53:48 michael Exp $
+*   $Log: bitfile.h,v $
+*   Revision 1.6  2007/08/26 21:53:48  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.5  2006/06/03 19:33:11  michael
+*   Used spell checker to correct spelling.
+*
+*   Revision 1.4  2005/12/06 15:06:37  michael
+*   Added BitFileGetBitsInt and BitFilePutBitsInt for integer types.
+*
+*   Revision 1.3  2004/11/09 14:16:58  michael
+*   Added functions to convert open bit_file_t to FILE and to
+*   align open bit_file_t to the next byte.
+*
+*   Revision 1.2  2004/06/15 13:16:10  michael
+*   Use incomplete type to hide definition of bitfile structure
+*
+*   Revision 1.1.1.1  2004/02/09 05:31:42  michael
+*   Initial release
+*
+*
+****************************************************************************
+*
+* Bitfile: Bit stream File I/O Routines
+* Copyright (C) 2004-2007 by Michael Dipperstein (mdipper@cs.ucsb.edu)
+*
+* This file is part of the bit file library.
+*
+* The bit file library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The bit file library is distributed in the hope that it will be useful,
+* but WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+#ifndef _BITFILE_H_
+#define _BITFILE_H_
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+typedef enum
+{
+    BF_READ = 0,
+    BF_WRITE = 1,
+    BF_APPEND= 2,
+    BF_NO_MODE
+} BF_MODES;
+
+/* incomplete type to hide implementation */
+struct bit_file_t;
+typedef struct bit_file_t bit_file_t;
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/* open/close file */
+bit_file_t *BitFileOpen(const char *fileName, const BF_MODES mode);
+bit_file_t *MakeBitFile(FILE *stream, const BF_MODES mode);
+int BitFileClose(bit_file_t *stream);
+FILE *BitFileToFILE(bit_file_t *stream);
+
+/* toss spare bits and byte align file */
+int BitFileByteAlign(bit_file_t *stream);
+
+/* get/put character */
+int BitFileGetChar(bit_file_t *stream);
+int BitFilePutChar(const int c, bit_file_t *stream);
+
+/* get/put single bit */
+int BitFileGetBit(bit_file_t *stream);
+int BitFilePutBit(const int c, bit_file_t *stream);
+
+/* get/put number of bits (most significant bit to least significat bit) */
+int BitFileGetBits(bit_file_t *stream, void *bits, const unsigned int count);
+int BitFilePutBits(bit_file_t *stream, void *bits, const unsigned int count);
+
+/***************************************************************************
+* get/put number of bits from integer types (short, int, long, ...)
+* machine endiness is accounted for.
+* size is the size of the data structure pointer to by bits.
+***************************************************************************/
+int BitFileGetBitsInt(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size);
+int BitFilePutBitsInt(bit_file_t *stream, void *bits, const unsigned int count,
+    const size_t size);
+
+#endif /* _BITFILE_H_ */
diff --git a/test/compression/lzdecode.c b/test/compression/lzdecode.c
new file mode 100644
index 000000000..cb46301e3
--- /dev/null
+++ b/test/compression/lzdecode.c
@@ -0,0 +1,282 @@
+/***************************************************************************
+*                 Lempel, Ziv, Storer, and Szymanski Decoding
+*
+*   File    : lzdecode.c
+*   Purpose : Use lzss coding (Storer and Szymanski's modified LZ77) to
+*             decompress lzss  encoded files.
+*   Author  : Michael Dipperstein
+*   Date    : November 07, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   Date        Change
+*   12/10/03    Changed handling of sliding window to better match standard
+*               algorithm description.
+*   12/11/03    Remebered to copy encoded characters to the sliding window
+*               even when there are no more characters in the input stream.
+*
+*
+*   Revision 1.2  2004/02/22 17:14:26  michael
+*   - Separated encode/decode, match finding, and main.
+*   - Use bitfiles for reading/writing files
+*   - Use traditional LZSS encoding where the coded/uncoded bits
+*     precede the symbol they are associated with, rather than
+*     aggregating the bits.
+*
+*   Revision 1.1.1.1  2004/01/21 06:25:49  michael
+*   Initial version
+*
+*   11/07/04    Separated encode and decode functions for improved
+*               modularity.
+*
+*   $Id: lzdecode.c,v 1.7 2007/09/20 04:34:25 michael Exp $
+*   $Log: lzdecode.c,v $
+*   Revision 1.7  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.6  2007/03/25 05:11:32  michael
+*   Corrected file closure error reported by "Carl@Yahoo" .  Now  both input
+*   and output files are closed.
+*
+*   Revision 1.5  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.4  2006/12/26 02:05:00  michael
+*   Corrected bug identified by Andrej Sinicyn which resulted in stdin being
+*   used as the default output for decoded data.
+*
+*   Revision 1.3  2005/12/28 06:03:30  michael
+*   Use slower but clearer Get/PutBitsInt for reading/writing bits.
+*   Replace mod with conditional Wrap macro.
+*
+*   Revision 1.2  2004/11/13 22:51:00  michael
+*   Provide distinct names for by file and by name functions and add some
+*   comments to make their usage clearer.
+*
+*   Revision 1.1  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*
+*
+****************************************************************************
+*
+* LZDecode: An ANSI C LZSS Decoding Routines
+* Copyright (C) 2003-2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "lzlocal.h"
+#include "bitfile.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+/* cyclic buffer sliding window of already read characters */
+extern unsigned char slidingWindow[];
+extern unsigned char uncodedLookahead[];
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : DecodeLZSSByFile
+*   Description: This function will read an LZSS encoded input file and
+*                write an output file.  This algorithm encodes strings as 16
+*                bits (a 12 bit offset + a 4 bit length).
+*   Parameters : fpIn - pointer to the open binary file to decode
+*                fpOut - pointer to the open binary file to write decoded
+*                       output
+*   Effects    : fpIn is decoded and written to fpOut.  Neither file is
+*                closed after exit.
+*   Returned   : EXIT_SUCCESS or EXIT_FAILURE
+****************************************************************************/
+int DecodeLZSSByFile(FILE *fpIn, FILE *fpOut)
+{
+    bit_file_t *bfpIn;
+
+    int c;
+    unsigned int i, nextChar;
+    encoded_string_t code;              /* offset/length code for string */
+
+    if (fpIn == NULL)
+    {
+        /* use stdin if no input file */
+        bfpIn = MakeBitFile(stdin, BF_READ);
+    }
+    else
+    {
+        /* convert output file to bitfile */
+        bfpIn = MakeBitFile(fpIn, BF_READ);
+    }
+
+    /* use stdout if no output file */
+    if (fpOut == NULL)
+    {
+        fpOut = stdout;
+    }
+
+    /************************************************************************
+    * Fill the sliding window buffer with some known vales.  EncodeLZSS must
+    * use the same values.  If common characters are used, there's an
+    * increased chance of matching to the earlier strings.
+    ************************************************************************/
+    memset(slidingWindow, ' ', WINDOW_SIZE * sizeof(unsigned char));
+
+    nextChar = 0;
+
+    while (TRUE)
+    {
+        if ((c = BitFileGetBit(bfpIn)) == EOF)
+        {
+            /* we hit the EOF */
+            break;
+        }
+
+        if (c == UNCODED)
+        {
+            /* uncoded character */
+            if ((c = BitFileGetChar(bfpIn)) == EOF)
+            {
+                break;
+            }
+
+            /* write out byte and put it in sliding window */
+            putc(c, fpOut);
+            slidingWindow[nextChar] = c;
+            nextChar = Wrap((nextChar + 1), WINDOW_SIZE);
+        }
+        else
+        {
+            /* offset and length */
+            code.offset = 0;
+            code.length = 0;
+
+            if ((BitFileGetBitsInt(bfpIn, &code.offset, OFFSET_BITS,
+                sizeof(unsigned int))) == EOF)
+            {
+                break;
+            }
+
+            if ((BitFileGetBitsInt(bfpIn, &code.length, LENGTH_BITS,
+                sizeof(unsigned int))) == EOF)
+            {
+                break;
+            }
+
+            code.length += MAX_UNCODED + 1;
+
+            /****************************************************************
+            * Write out decoded string to file and lookahead.  It would be
+            * nice to write to the sliding window instead of the lookahead,
+            * but we could end up overwriting the matching string with the
+            * new string if abs(offset - next char) < match length.
+            ****************************************************************/
+            for (i = 0; i < code.length; i++)
+            {
+                c = slidingWindow[Wrap((code.offset + i), WINDOW_SIZE)];
+                putc(c, fpOut);
+                uncodedLookahead[i] = c;
+            }
+
+            /* write out decoded string to sliding window */
+            for (i = 0; i < code.length; i++)
+            {
+                slidingWindow[(nextChar + i) % WINDOW_SIZE] =
+                    uncodedLookahead[i];
+            }
+
+            nextChar = Wrap((nextChar + code.length), WINDOW_SIZE);
+        }
+    }
+
+    /* we've decoded everything, free bitfile structure */
+    BitFileToFILE(bfpIn);
+
+    return (EXIT_SUCCESS);
+}
+
+/****************************************************************************
+*   Function   : DecodeLZSSByName
+*   Description: This function is provided to maintain compatibility with
+*                older versions of my LZSS implementation which used file
+*                names instead of file pointers.
+*   Parameters : inFile - name of file to decode
+*                outFile - name of file to write decoded output
+*   Effects    : inFile is decoded and written to outFile
+*   Returned   : EXIT_SUCCESS or EXIT_FAILURE
+****************************************************************************/
+int DecodeLZSSByName(char *inFile, char *outFile)
+{
+    FILE *fpIn, *fpOut;
+    int returnValue;
+
+    if (inFile == NULL)
+    {
+        fpIn = stdin;
+    }
+    if ((fpIn = fopen(inFile, "rb")) == NULL)
+    {
+        perror(inFile);
+        return (EXIT_FAILURE);
+    }
+
+    if (outFile == NULL)
+    {
+        fpOut = stdout;
+    }
+    else
+    {
+        if ((fpOut = fopen(outFile, "wb")) == NULL)
+        {
+            fclose(fpIn);
+            perror(outFile);
+            return (EXIT_FAILURE);
+        }
+    }
+
+    returnValue = DecodeLZSSByFile(fpIn, fpOut);
+
+    /* close files */
+    fclose(fpIn);
+    fclose(fpOut);
+
+    return (returnValue);
+}
diff --git a/test/compression/lzencode.c b/test/compression/lzencode.c
new file mode 100644
index 000000000..35f8b8345
--- /dev/null
+++ b/test/compression/lzencode.c
@@ -0,0 +1,300 @@
+/***************************************************************************
+*                 Lempel, Ziv, Storer, and Szymanski Encoding
+*
+*   File    : lzdecode.c
+*   Purpose : Use lzss coding (Storer and Szymanski's modified LZ77) to
+*             compress lzss data files.
+*   Author  : Michael Dipperstein
+*   Date    : November 07, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   Date        Change
+*   12/10/03    Changed handling of sliding window to better match standard
+*               algorithm description.
+*   12/11/03    Remebered to copy encoded characters to the sliding window
+*               even when there are no more characters in the input stream.
+*
+*
+*   Revision 1.2  2004/02/22 17:14:26  michael
+*   - Separated encode/decode, match finding, and main.
+*   - Use bitfiles for reading/writing files
+*   - Use traditional LZSS encoding where the coded/uncoded bits
+*     precede the symbol they are associated with, rather than
+*     aggregating the bits.
+*
+*   Revision 1.1.1.1  2004/01/21 06:25:49  michael
+*   Initial version
+*
+*   11/07/04    Separated encode and decode functions for improved
+*               modularity.
+*
+*   $Id: lzencode.c,v 1.6 2007/09/20 04:34:25 michael Exp $
+*   $Log: lzencode.c,v $
+*   Revision 1.6  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.5  2007/03/25 05:11:32  michael
+*   Corrected file closure error reported by "Carl@Yahoo" .  Now  both input
+*   and output files are closed.
+*
+*   Revision 1.4  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.3  2005/12/28 06:03:30  michael
+*   Use slower but clearer Get/PutBitsInt for reading/writing bits.
+*   Replace mod with conditional Wrap macro.
+*
+*   Revision 1.2  2004/11/13 22:51:00  michael
+*   Provide distinct names for by file and by name functions and add some
+*   comments to make their usage clearer.
+*
+*   Revision 1.1  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*
+*
+****************************************************************************
+*
+* LZEncode: An ANSI C LZSS Encoding Routines
+* Copyright (C) 2003-2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "lzlocal.h"
+#include "bitfile.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+/* cyclic buffer sliding window of already read characters */
+extern unsigned char slidingWindow[];
+extern unsigned char uncodedLookahead[];
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : EncodeLZSSByFile
+*   Description: This function will read an input file and write an output
+*                file encoded according to the traditional LZSS algorithm.
+*                This algorithm encodes strings as 16 bits (a 12 bit offset
+*                + a 4 bit length).
+*   Parameters : fpIn - pointer to the open binary file to encode
+*                fpOut - pointer to the open binary file to write encoded
+*                       output
+*   Effects    : fpIn is encoded and written to fpOut.  Neither file is
+*                closed after exit.
+*   Returned   : EXIT_SUCCESS or EXIT_FAILURE
+****************************************************************************/
+int EncodeLZSSByFile(FILE *fpIn, FILE *fpOut)
+{
+    bit_file_t *bfpOut;
+
+    encoded_string_t matchData;
+    unsigned int i, c;
+    unsigned int len;                       /* length of string */
+
+    /* head of sliding window and lookahead */
+    unsigned int windowHead, uncodedHead;
+
+    /* use stdin if no input file */
+    if (fpIn == NULL)
+    {
+        fpIn = stdin;
+    }
+
+    if (fpOut == NULL)
+    {
+        /* use stdout if no output file */
+        bfpOut = MakeBitFile(stdout, BF_WRITE);
+    }
+    else
+    {
+        /* convert output file to bitfile */
+        bfpOut = MakeBitFile(fpOut, BF_WRITE);
+    }
+
+    windowHead = 0;
+    uncodedHead = 0;
+
+    /************************************************************************
+    * Fill the sliding window buffer with some known vales.  DecodeLZSS must
+    * use the same values.  If common characters are used, there's an
+    * increased chance of matching to the earlier strings.
+    ************************************************************************/
+    memset(slidingWindow, ' ', WINDOW_SIZE * sizeof(unsigned char));
+
+    /************************************************************************
+    * Copy MAX_CODED bytes from the input file into the uncoded lookahead
+    * buffer.
+    ************************************************************************/
+    for (len = 0; len < MAX_CODED; len++) {
+        c = getc(fpIn);
+        if (c == EOF) break;
+        uncodedLookahead[len] = c;
+    }
+
+    if (len == 0)
+    {
+        return (EXIT_SUCCESS);   /* inFile was empty */
+    }
+
+    /* Look for matching string in sliding window */
+    InitializeSearchStructures();
+    FindMatch(&matchData, windowHead, uncodedHead);
+
+    /* now encoded the rest of the file until an EOF is read */
+    while (len > 0)
+    {
+        if (matchData.length > len)
+        {
+            /* garbage beyond last data happened to extend match length */
+            matchData.length = len;
+        }
+
+        if (matchData.length <= MAX_UNCODED)
+        {
+            /* not long enough match.  write uncoded flag and character */
+            BitFilePutBit(UNCODED, bfpOut);
+            BitFilePutChar(uncodedLookahead[uncodedHead], bfpOut);
+
+            matchData.length = 1;   /* set to 1 for 1 byte uncoded */
+        }
+        else
+        {
+            unsigned int adjustedLen;
+
+            /* adjust the length of the match so minimun encoded len is 0*/
+            adjustedLen = matchData.length - (MAX_UNCODED + 1);
+
+            /* match length > MAX_UNCODED.  Encode as offset and length. */
+            BitFilePutBit(ENCODED, bfpOut);
+            BitFilePutBitsInt(bfpOut, &matchData.offset, OFFSET_BITS,
+                sizeof(unsigned int));
+            BitFilePutBitsInt(bfpOut, &adjustedLen, LENGTH_BITS,
+                sizeof(unsigned int));
+        }
+
+        /********************************************************************
+        * Replace the matchData.length worth of bytes we've matched in the
+        * sliding window with new bytes from the input file.
+        ********************************************************************/
+        i = 0;
+        while ((i < matchData.length) && ((c = getc(fpIn)) != EOF))
+        {
+            /* add old byte into sliding window and new into lookahead */
+            ReplaceChar(windowHead, uncodedLookahead[uncodedHead]);
+            uncodedLookahead[uncodedHead] = c;
+            windowHead = Wrap((windowHead + 1), WINDOW_SIZE);
+            uncodedHead = Wrap((uncodedHead + 1), MAX_CODED);
+            i++;
+        }
+
+        /* handle case where we hit EOF before filling lookahead */
+        while (i < matchData.length)
+        {
+            ReplaceChar(windowHead, uncodedLookahead[uncodedHead]);
+            /* nothing to add to lookahead here */
+            windowHead = Wrap((windowHead + 1), WINDOW_SIZE);
+            uncodedHead = Wrap((uncodedHead + 1), MAX_CODED);
+            len--;
+            i++;
+        }
+
+        /* find match for the remaining characters */
+        FindMatch(&matchData, windowHead, uncodedHead);
+    }
+
+    /* we've decoded everything, free bitfile structure */
+    BitFileToFILE(bfpOut);
+
+   return (EXIT_SUCCESS);
+}
+
+/****************************************************************************
+*   Function   : EncodeLZSSByName
+*   Description: This function is provided to maintain compatibility with
+*                older versions of my LZSS implementation which used file
+*                names instead of file pointers.
+*   Parameters : inFile - name of file to encode
+*                outFile - name of file to write encoded output
+*   Effects    : inFile is encoded and written to outFile
+*   Returned   : EXIT_SUCCESS or EXIT_FAILURE
+****************************************************************************/
+int EncodeLZSSByName(char *inFile, char *outFile)
+{
+    FILE *fpIn, *fpOut;
+    int returnValue;
+
+    /* open binary input and output files */
+    if (inFile == NULL)
+    {
+        fpIn = stdin;
+    }
+    if ((fpIn = fopen(inFile, "rb")) == NULL)
+    {
+        perror(inFile);
+        return (EXIT_FAILURE);
+    }
+
+    if (outFile == NULL)
+    {
+        fpOut = stdout;
+    }
+    else
+    {
+        if ((fpOut = fopen(outFile, "wb")) == NULL)
+        {
+            fclose(fpIn);
+            perror(outFile);
+            return (EXIT_FAILURE);
+        }
+    }
+
+    returnValue = EncodeLZSSByFile(fpIn, fpOut);
+
+    /* close files */
+    fclose(fpIn);
+    fclose(fpOut);
+
+    return (returnValue);
+}
diff --git a/test/compression/lzhash.c b/test/compression/lzhash.c
new file mode 100644
index 000000000..83cbe94df
--- /dev/null
+++ b/test/compression/lzhash.c
@@ -0,0 +1,358 @@
+#include <stdio.h>
+/***************************************************************************
+*          Lempel, Ziv, Storer, and Szymanski Encoding and Decoding
+*
+*   File    : hash.c
+*   Purpose : Implement hash table optimized matching of uncoded strings
+*             for LZSS algorithm.
+*   Author  : Michael Dipperstein
+*   Date    : February 21, 2004
+*
+****************************************************************************
+*   UPDATES
+*   $Id: hash.c,v 1.7 2007/09/20 04:34:25 michael Exp $
+*   $Log: hash.c,v $
+*   Revision 1.7  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.6  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.5  2005/12/29 14:37:56  michael
+*   Remove debug statements.
+*
+*   Revision 1.4  2005/12/29 14:32:56  michael
+*   Correct algorithm for replacing characters in dictionary.
+*
+*   Revision 1.3  2005/12/29 07:06:24  michael
+*   Let the hash table size vary with the size of the sliding window.
+*
+*   Revision 1.2  2005/12/28 06:03:30  michael
+*   Use slower but clearer Get/PutBitsInt for reading/writing bits.
+*   Replace mod with conditional Wrap macro.
+*
+*   Revision 1.1  2004/02/22 17:23:02  michael
+*   Initial revision of hash table search.  Mostly code from lzhash.c.
+*
+*
+****************************************************************************
+*
+* Hash: Hash table optimized matching routines used by LZSS
+*       Encoding/Decoding Routine
+* Copyright (C) 2004-2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include "lzlocal.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#define NULL_INDEX      (WINDOW_SIZE + 1)
+
+#define HASH_SIZE       (WINDOW_SIZE >> 2)  /* size of hash table */
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+/* cyclic buffer sliding window of already read characters */
+extern unsigned char slidingWindow[];
+extern unsigned char uncodedLookahead[];
+
+unsigned int hashTable[HASH_SIZE];          /* list head for each hash key */
+unsigned int next[WINDOW_SIZE];             /* indices of next in hash list */
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : HashKey
+*   Description: This function generates a hash key for a (MAX_UNCODED + 1)
+*                long string located either in the sliding window or the
+*                uncoded lookahead.  The key generation algorithm is
+*                supposed to be based on the algorithm used by gzip.  As
+*                reported in K. Sadakane, H. Imai. "Improving the Speed of
+*                LZ77 Compression by Hashing and Suffix Sorting". IEICE
+*                Trans. Fundamentals, Vol. E83-A, No. 12 (December 2000)
+*   Parameters : offset - offset into either the uncoded lookahead or the
+*                         sliding window.
+*                lookahead - TRUE iff offset is an offset into the uncoded
+*                            lookahead buffer.
+*   Effects    : NONE
+*   Returned   : The sliding window index where the match starts and the
+*                length of the match.  If there is no match a length of
+*                zero will be returned.
+****************************************************************************/
+unsigned int HashKey(unsigned int offset, unsigned char lookahead)
+{
+    int i, hashKey;
+
+    hashKey = 0;
+
+    if (lookahead)
+    {
+        /* string is in the lookahead buffer */
+        for (i = 0; i < (MAX_UNCODED + 1); i++)
+        {
+            hashKey = (hashKey << 5) ^ (uncodedLookahead[offset]);
+            hashKey %= HASH_SIZE;
+            offset = Wrap((offset + 1), MAX_CODED);
+        }
+    }
+    else
+    {
+        /* string is in the sliding window */
+        for (i = 0; i < (MAX_UNCODED + 1); i++)
+        {
+            hashKey = (hashKey << 5) ^ (slidingWindow[offset]);
+            hashKey %= HASH_SIZE;
+            offset = Wrap((offset + 1), WINDOW_SIZE);
+        }
+    }
+
+    return hashKey;
+}
+
+/****************************************************************************
+*   Function   : InitializeSearchStructures
+*   Description: This function initializes structures used to speed up the
+*                process of mathcing uncoded strings to strings in the
+*                sliding window.  For hashed searches, this means that a
+*                hash table pointing to linked lists is initialized.
+*   Parameters : None
+*   Effects    : The hash table and next array are initialized.
+*   Returned   : None
+*
+*   NOTE: This function assumes that the sliding window is initially filled
+*         with all identical characters.
+****************************************************************************/
+void InitializeSearchStructures()
+{
+    unsigned int i;
+
+    /************************************************************************
+    * Since the encode routine only fills the sliding window with one
+    * character, there is only one hash key for the entier sliding window.
+    * That means all positions are in the same linked list.
+    ************************************************************************/
+    for (i = 0; i < (WINDOW_SIZE - 1); i++)
+    {
+        next[i] = i + 1;
+    }
+
+    /* there's no next for the last character */
+    next[WINDOW_SIZE - 1] = NULL_INDEX;
+
+    /* the only list right now is the "   " list */
+    for (i = 0; i < HASH_SIZE; i++)
+    {
+        hashTable[i] = NULL_INDEX;
+    }
+
+    hashTable[HashKey(0, FALSE)] = 0;
+}
+
+/****************************************************************************
+*   Function   : FindMatch
+*   Description: This function will search through the slidingWindow
+*                dictionary for the longest sequence matching the MAX_CODED
+*                long string stored in uncodedLookahead.
+*   Parameters : uncodedHead - head of uncoded lookahead buffer
+*   Effects    : NONE
+*   Returned   : The sliding window index where the match starts and the
+*                length of the match.  If there is no match a length of
+*                zero will be returned.
+****************************************************************************/
+/* XL: no struct return */
+void FindMatch(encoded_string_t * matchData, 
+               unsigned int windowHead, unsigned int uncodedHead)
+{
+    unsigned int i, j;
+
+    matchData->length = 0;
+    matchData->offset = 0;
+
+    i = hashTable[HashKey(uncodedHead, TRUE)];  /* start of proper list */
+    j = 0;
+
+    while (i != NULL_INDEX)
+    {
+        if (slidingWindow[i] == uncodedLookahead[uncodedHead])
+        {
+            /* we matched one how many more match? */
+            j = 1;
+
+            while(slidingWindow[Wrap((i + j), WINDOW_SIZE)] ==
+                uncodedLookahead[Wrap((uncodedHead + j), MAX_CODED)])
+            {
+                if (j >= MAX_CODED)
+                {
+                    break;
+                }
+                j++;
+            }
+
+            if (j > matchData->length)
+            {
+                matchData->length = j;
+                matchData->offset = i;
+            }
+        }
+
+        if (j >= MAX_CODED)
+        {
+            matchData->length = MAX_CODED;
+            break;
+        }
+
+        i = next[i];    /* try next in list */
+    }
+}
+
+/****************************************************************************
+*   Function   : AddString
+*   Description: This function adds the (MAX_UNCODED + 1) long string
+*                starting at slidingWindow[charIndex] to the hash table's
+*                linked list associated with its hash key.
+*   Parameters : charIndex - sliding window index of the string to be
+*                            added to the linked list.
+*   Effects    : The string starting at slidingWindow[charIndex] is appended
+*                to the end of the appropriate linked list.
+*   Returned   : NONE
+****************************************************************************/
+void AddString(unsigned int charIndex)
+{
+    unsigned int i, hashKey;
+
+    /* inserted character will be at the end of the list */
+    next[charIndex] = NULL_INDEX;
+
+    hashKey = HashKey(charIndex, FALSE);
+
+    if (hashTable[hashKey] == NULL_INDEX)
+    {
+        /* this is the only character in it's list */
+        hashTable[hashKey] = charIndex;
+        return;
+    }
+
+    /* find the end of the list */
+    i = hashTable[hashKey];
+
+    while(next[i] != NULL_INDEX)
+    {
+        i = next[i];
+    }
+
+    /* add new character to the list end */
+    next[i] = charIndex;
+}
+
+/****************************************************************************
+*   Function   : RemoveString
+*   Description: This function removes the (MAX_UNCODED + 1) long string
+*                starting at slidingWindow[charIndex] from the hash table's
+*                linked list associated with its hash key.
+*   Parameters : charIndex - sliding window index of the string to be
+*                            removed from the linked list.
+*   Effects    : The string starting at slidingWindow[charIndex] is removed
+*                from its linked list.
+*   Returned   : NONE
+****************************************************************************/
+void RemoveString(unsigned int charIndex)
+{
+    unsigned int i, hashKey;
+    unsigned int nextIndex;
+
+    nextIndex = next[charIndex];        /* remember where this points to */
+    next[charIndex] = NULL_INDEX;
+
+    hashKey = HashKey(charIndex, FALSE);
+
+    if (hashTable[hashKey] == charIndex)
+    {
+        /* we're deleting a list head */
+        hashTable[hashKey] = nextIndex;
+        return;
+    }
+
+    /* find character pointing to ours */
+    i = hashTable[hashKey];
+
+    while(next[i] != charIndex)
+    {
+        i = next[i];
+    }
+
+    /* point the previous next */
+    next[i] = nextIndex;
+}
+
+/****************************************************************************
+*   Function   : ReplaceChar
+*   Description: This function replaces the character stored in
+*                slidingWindow[charIndex] with the one specified by
+*                replacement.  The hash table entries effected by the
+*                replacement are also corrected.
+*   Parameters : charIndex - sliding window index of the character to be
+*                            removed from the linked list.
+*   Effects    : slidingWindow[charIndex] is replaced by replacement.  Old
+*                hash entries for strings containing slidingWindow[charIndex]
+*                are removed and new ones are added.
+*   Returned   : NONE
+****************************************************************************/
+void ReplaceChar(unsigned int charIndex, unsigned char replacement)
+{
+    unsigned int firstIndex, i;
+
+    if (charIndex < MAX_UNCODED)
+    {
+        firstIndex = (WINDOW_SIZE + charIndex) - MAX_UNCODED;
+    }
+    else
+    {
+        firstIndex = charIndex - MAX_UNCODED;
+    }
+
+    /* remove all hash entries containing character at char index */
+    for (i = 0; i < (MAX_UNCODED + 1); i++)
+    {
+        RemoveString(Wrap((firstIndex + i), WINDOW_SIZE));
+    }
+
+    slidingWindow[charIndex] = replacement;
+
+    /* add all hash entries containing character at char index */
+    for (i = 0; i < (MAX_UNCODED + 1); i++)
+    {
+        AddString(Wrap((firstIndex + i), WINDOW_SIZE));
+    }
+}
diff --git a/test/compression/lzlocal.h b/test/compression/lzlocal.h
new file mode 100644
index 000000000..876de9172
--- /dev/null
+++ b/test/compression/lzlocal.h
@@ -0,0 +1,123 @@
+/***************************************************************************
+*          Lempel, Ziv, Storer, and Szymanski Encoding and Decoding
+*
+*   File    : lzlocal.h
+*   Purpose : Internal headers for LZSS encode and decode routines.
+*             Contains the prototypes to be used by the different match
+*             finding algorithms.
+*   Author  : Michael Dipperstein
+*   Date    : February 18, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzlocal.h,v 1.5 2007/09/20 04:34:25 michael Exp $
+*   $Log: lzlocal.h,v $
+*   Revision 1.5  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.4  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.3  2005/12/28 05:58:52  michael
+*   Add Wrap macro to replace mod when value is less than twice the limit.
+*
+*   Revision 1.2  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*   Revision 1.1  2004/02/22 17:32:40  michael
+*   Initial revision of header files for sliding window search implementations.
+*
+*
+****************************************************************************
+*
+* LZSS: An ANSI C LZSS Encoding/Decoding Routine
+* Copyright (C) 2004-2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+#ifndef _LZSS_LOCAL_H
+#define _LZSS_LOCAL_H
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <limits.h>
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#define FALSE   0
+#define TRUE    1
+
+#define OFFSET_BITS     12
+#define LENGTH_BITS     4
+
+#if (((1 << (OFFSET_BITS + LENGTH_BITS)) - 1) > UINT_MAX)
+#error "Size of encoded data must not exceed the size of an unsigned int"
+#endif
+
+/* We want a sliding window*/
+#define WINDOW_SIZE     (1 << OFFSET_BITS)
+
+/* maximum match length not encoded and maximum length encoded (4 bits) */
+#define MAX_UNCODED     2
+#define MAX_CODED       ((1 << LENGTH_BITS) + MAX_UNCODED)
+
+#define ENCODED     0       /* encoded string */
+#define UNCODED     1       /* unencoded character */
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+* This data structure stores an encoded string in (offset, length) format.
+* The actual encoded string is stored using OFFSET_BITS for the offset and
+* LENGTH_BITS for the length.
+***************************************************************************/
+typedef struct encoded_string_t
+{
+    unsigned int offset;    /* offset to start of longest match */
+    unsigned int length;    /* length of longest match */
+} encoded_string_t;
+
+/***************************************************************************
+*                                 MACROS
+***************************************************************************/
+/* wraps array index within array bounds (assumes value < 2 * limit) */
+#if 0
+#define Wrap(value, limit)      (((value) < (limit)) ? (value) : ((value) - (limit)))
+#else
+#define Wrap(value, limit)   ((unsigned int)(value) % (unsigned int)(limit))
+#endif
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+void InitializeSearchStructures(void);
+
+/* XL: no struct return */
+void FindMatch(encoded_string_t * /*out*/,
+               unsigned int windowHead, unsigned int uncodedHead);
+
+void ReplaceChar(unsigned int charIndex, unsigned char replacement);
+
+#endif      /* ndef _LZSS_LOCAL_H */
diff --git a/test/compression/lzss.h b/test/compression/lzss.h
new file mode 100644
index 000000000..990d893e0
--- /dev/null
+++ b/test/compression/lzss.h
@@ -0,0 +1,96 @@
+/***************************************************************************
+*          Lempel, Ziv, Storer, and Szymanski Encoding and Decoding
+*
+*   File    : lzss.h
+*   Purpose : Header for LZSS encode and decode routines.  Contains the
+*             prototypes to be used by programs linking to the LZSS
+*             library.
+*   Author  : Michael Dipperstein
+*   Date    : February 21, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzss.h,v 1.5 2007/09/20 04:34:25 michael Exp $
+*   $Log: lzss.h,v $
+*   Revision 1.5  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.4  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.3  2004/11/13 22:51:00  michael
+*   Provide distinct names for by file and by name functions and add some
+*   comments to make their usage clearer.
+*
+*   Revision 1.2  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*   Revision 1.1  2004/02/22 17:37:50  michael
+*   Initial revision of headers for encode and decode routines.
+*
+*
+****************************************************************************
+*
+* LZSS: An ANSI C LZSS Encoding/Decoding Routine
+* Copyright (C) 2004, 2006, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+#ifndef _LZSS_H
+#define _LZSS_H
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+
+/***************************************************************************
+*                                 MACROS
+***************************************************************************/
+/* macros for compatibility with older library */
+#define EncodeLZSS(in, out)     EncodeLZSSByName((in), (out))
+#define DecodeLZSS(in, out)     DecodeLZSSByName((in), (out))
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+/***************************************************************************
+* LZSS encoding and decoding prototypes for functions with file name
+* parameters.  Provide these functions with name of the file to be
+* encoded/decoded (inFile) and the name of the target file (outFile).
+* These functions return EXIT_SUCCESS or EXIT_FAILURE.
+***************************************************************************/
+int EncodeLZSSByName(char *inFile, char *outFile);
+int DecodeLZSSByName(char *inFile, char *outFile);
+
+/***************************************************************************
+* LZSS encoding and decoding prototypes for functions with file pointer
+* parameters.  Provide these functions with a pointer to the open binary
+* file to be encoded/decoded (fpIn) and pointer to the open binary target
+* file (fpOut).  It is the job of the function caller to open the files
+* prior to callings these functions and to close the file after these
+* functions have been called.
+* These functions return EXIT_SUCCESS or EXIT_FAILURE.
+***************************************************************************/
+int EncodeLZSSByFile(FILE *fpIn, FILE *fpOut);
+int DecodeLZSSByFile(FILE *fpIn, FILE *fpOut);
+
+#endif      /* ndef _LZSS_H */
diff --git a/test/compression/lzssmain.c b/test/compression/lzssmain.c
new file mode 100644
index 000000000..30fc5ab24
--- /dev/null
+++ b/test/compression/lzssmain.c
@@ -0,0 +1,278 @@
+/***************************************************************************
+*                     Sample Program Using LZSS Library
+*
+*   File    : sample.c
+*   Purpose : Demonstrate usage of LZSS library
+*   Author  : Michael Dipperstein
+*   Date    : February 21, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   Revision 1.1  2004/02/22 17:36:30  michael
+*   Initial revision.  Mostly code form old lzss.c.
+*
+*   11/07/04    Name changed to sample.c
+*
+*   $Id: sample.c,v 1.5 2007/09/20 04:34:45 michael Exp $
+*   $Log: sample.c,v $
+*   Revision 1.5  2007/09/20 04:34:45  michael
+*   Replace getopt with optlist.
+*   Changes required for LGPL v3.
+*
+*   Revision 1.4  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.3  2004/11/13 22:51:01  michael
+*   Provide distinct names for by file and by name functions and add some
+*   comments to make their usage clearer.
+*
+*   Revision 1.2  2004/11/11 14:37:26  michael
+*   Open input and output files as binary.
+*
+*   Revision 1.1  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*
+****************************************************************************
+*
+* SAMPLE: Sample usage of LZSS Library
+* Copyright (C) 2004, 2006, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "lzss.h"
+#include "optlist.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+typedef enum
+{
+    ENCODE,
+    DECODE
+} MODES;
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+char *RemovePath(char *fullPath);
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : main
+*   Description: This is the main function for this program, it validates
+*                the command line input and, if valid, it will either
+*                encode a file using the LZSS algorithm or decode a
+*                file encoded with the LZSS algorithm.
+*   Parameters : argc - number of parameters
+*                argv - parameter list
+*   Effects    : Encodes/Decodes input file
+*   Returned   : EXIT_SUCCESS for success, otherwise EXIT_FAILURE.
+****************************************************************************/
+int main(int argc, char *argv[])
+{
+    option_t *optList, *thisOpt;
+    FILE *fpIn, *fpOut;      /* pointer to open input & output files */
+    MODES mode;
+
+    /* initialize data */
+    fpIn = NULL;
+    fpOut = NULL;
+    mode = ENCODE;
+
+    /* parse command line */
+    optList = GetOptList(argc, argv, "cdi:o:h?");
+    thisOpt = optList;
+
+    while (thisOpt != NULL)
+    {
+        switch(thisOpt->option)
+        {
+            case 'c':       /* compression mode */
+                mode = ENCODE;
+                break;
+
+            case 'd':       /* decompression mode */
+                mode = DECODE;
+                break;
+
+            case 'i':       /* input file name */
+                if (fpIn != NULL)
+                {
+                    fprintf(stderr, "Multiple input files not allowed.\n");
+                    fclose(fpIn);
+
+                    if (fpOut != NULL)
+                    {
+                        fclose(fpOut);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                /* open input file as binary */
+                fpIn = fopen(thisOpt->argument, "rb");
+                if (fpIn == NULL)
+                {
+                    perror("Opening input file");
+
+                    if (fpOut != NULL)
+                    {
+                        fclose(fpOut);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                break;
+
+            case 'o':       /* output file name */
+                if (fpOut != NULL)
+                {
+                    fprintf(stderr, "Multiple output files not allowed.\n");
+                    fclose(fpOut);
+
+                    if (fpIn != NULL)
+                    {
+                        fclose(fpIn);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                /* open output file as binary */
+                fpOut = fopen(thisOpt->argument, "wb");
+                if (fpOut == NULL)
+                {
+                    perror("Opening output file");
+
+                    if (fpIn != NULL)
+                    {
+                        fclose(fpIn);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                break;
+
+            case 'h':
+            case '?':
+                printf("Usage: %s <options>\n\n", RemovePath(argv[0]));
+                printf("options:\n");
+                printf("  -c : Encode input file to output file.\n");
+                printf("  -d : Decode input file to output file.\n");
+                printf("  -i <filename> : Name of input file.\n");
+                printf("  -o <filename> : Name of output file.\n");
+                printf("  -h | ?  : Print out command line options.\n\n");
+                printf("Default: %s -c -i stdin -o stdout\n",
+                    RemovePath(argv[0]));
+
+                FreeOptList(optList);
+                return(EXIT_SUCCESS);
+        }
+
+        optList = thisOpt->next;
+        free(thisOpt);
+        thisOpt = optList;
+    }
+
+
+    /* use stdin/out if no files are provided */
+    if (fpIn == NULL)
+    {
+        fpIn = stdin;
+    }
+
+
+    if (fpOut == NULL)
+    {
+        fpOut = stdout;
+    }
+
+    /* we have valid parameters encode or decode */
+    if (mode == ENCODE)
+    {
+        EncodeLZSSByFile(fpIn, fpOut);
+    }
+    else
+    {
+        DecodeLZSSByFile(fpIn, fpOut);
+    }
+
+    /* remember to close files */
+    fclose(fpIn);
+    fclose(fpOut);
+    return EXIT_SUCCESS;
+}
+
+/***************************************************************************
+*   Function   : RemovePath
+*   Description: This is function accepts a pointer to the name of a file
+*                along with path information and returns a pointer to the
+*                character that is not part of the path.
+*   Parameters : fullPath - pointer to an array of characters containing
+*                           a file name and possible path modifiers.
+*   Effects    : None
+*   Returned   : Returns a pointer to the first character after any path
+*                information.
+***************************************************************************/
+char *RemovePath(char *fullPath)
+{
+    int i;
+    char *start, *tmp;                          /* start of file name */
+    const char delim[3] = {'\\', '/', ':'};     /* path deliminators */
+
+    start = fullPath;
+
+    /* find the first character after all file path delimiters */
+    for (i = 0; i < 3; i++)
+    {
+        tmp = strrchr(start, delim[i]);
+
+        if (tmp != NULL)
+        {
+            start = tmp + 1;
+        }
+    }
+
+    return start;
+}
diff --git a/test/compression/lzvars.c b/test/compression/lzvars.c
new file mode 100644
index 000000000..cc82a6fbb
--- /dev/null
+++ b/test/compression/lzvars.c
@@ -0,0 +1,76 @@
+/***************************************************************************
+*            Lempel, Ziv, Storer, and Szymanski Global Variables
+*
+*   File    : lzvars.c
+*   Purpose : Declare global variables used lzss coding
+*             compress/decompress files.
+*   Author  : Michael Dipperstein
+*   Date    : November 07, 2003
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzvars.c,v 1.3 2007/09/20 04:34:25 michael Exp $
+*   $Log: lzvars.c,v $
+*   Revision 1.3  2007/09/20 04:34:25  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.2  2006/12/26 04:09:09  michael
+*   Updated e-mail address and minor text clean-up.
+*
+*   Revision 1.1  2004/11/08 05:54:18  michael
+*   1. Split encode and decode routines for smarter linking
+*   2. Renamed lzsample.c sample.c to match my other samples
+*   3. Makefile now builds code as libraries for better LGPL compliance.
+*
+*
+****************************************************************************
+*
+* LZvars: Global variables used in ANSI C LZSS Encoding/Decoding Routines
+* Copyright (C) 2003, 2004, 2006, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzss library.
+*
+* The lzss library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzss library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include "lzlocal.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+/* cyclic buffer sliding window of already read characters */
+unsigned char slidingWindow[WINDOW_SIZE];
+unsigned char uncodedLookahead[MAX_CODED];
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
diff --git a/test/compression/lzw.h b/test/compression/lzw.h
new file mode 100644
index 000000000..89255f060
--- /dev/null
+++ b/test/compression/lzw.h
@@ -0,0 +1,70 @@
+/***************************************************************************
+*          Header for Lempel-Ziv-Welch Encoding and Decoding Library
+*
+*   File    : arcode.h
+*   Purpose : Provides prototypes for functions that use Lempel-Ziv-Welch
+*             coding to encode/decode files.
+*   Author  : Michael Dipperstein
+*   Date    : January 30, 2004
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzw.h,v 1.3 2007/09/29 01:28:09 michael Exp $
+*   $Log: lzw.h,v $
+*   Revision 1.3  2007/09/29 01:28:09  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.2  2005/03/27 21:12:03  michael
+*   Update e-mail address in copyright notice.
+*
+*   Revision 1.1.1.1  2005/02/21 03:35:34  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* LZW: An ANSI C Lempel-Ziv-Welch Encoding/Decoding Routines
+* Copyright (C) 2005, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzw library.
+*
+* The lzw library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzw library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+#ifndef _LZW_H_
+#define _LZW_H_
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#ifndef FALSE
+#define FALSE       0
+#endif
+
+#ifndef TRUE
+#define TRUE        1
+#endif
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+ /* encode inFile */
+int LZWEncodeFile(char *inFile, char *outFile);
+
+/* decode inFile*/
+int LZWDecodeFile(char *inFile, char *outFile);
+
+#endif  /* ndef _ARCODE_H_ */
diff --git a/test/compression/lzwdecode.c b/test/compression/lzwdecode.c
new file mode 100644
index 000000000..1f8b5eca1
--- /dev/null
+++ b/test/compression/lzwdecode.c
@@ -0,0 +1,306 @@
+/***************************************************************************
+*                    Lempel-Ziv-Welch Decoding Functions
+*
+*   File    : lzwdecode.c
+*   Purpose : Provides a function for decoding Lempel-Ziv-Welch encoded
+*             file streams
+*   Author  : Michael Dipperstein
+*   Date    : January 30, 2005
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzwdecode.c,v 1.2 2007/09/29 01:28:09 michael Exp $
+*   $Log: lzwdecode.c,v $
+*   Revision 1.2  2007/09/29 01:28:09  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.1  2005/04/09 03:11:22  michael
+*   Separated encode and decode routines into two different files in order to
+*   make future enhancements easier.
+*
+*   Revision 1.4  2005/03/27 15:56:47  michael
+*   Use variable length code words.
+*   Include new e-mail addres.
+*
+*   Revision 1.3  2005/03/10 14:26:58  michael
+*   User variable names that match discription in web page.
+*
+*   Revision 1.2  2005/03/10 05:44:02  michael
+*   Minimize the size of the dictionary.
+*
+*   Revision 1.1.1.1  2005/02/21 03:35:34  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* LZW: An ANSI C Lempel-Ziv-Welch Encoding/Decoding Routines
+* Copyright (C) 2005, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzw library.
+*
+* The lzw library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzw library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <string.h>
+#include "lzw.h"
+#include "bitfile.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+typedef struct
+{
+    unsigned char suffixChar;   /* last char in encoded string */
+    unsigned int prefixCode;    /* code for remaining chars in string */
+} decode_dictionary_t;
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+#define EMPTY           -1
+
+#define MIN_CODE_LEN    9                   /* min # bits in a code word */
+#define MAX_CODE_LEN    15                  /* max # bits in a code word */
+
+#define FIRST_CODE      (1 << CHAR_BIT)     /* value of 1st string code */
+
+#define MAX_CODES       (1 << MAX_CODE_LEN)
+
+#define DICT_SIZE       (MAX_CODES - FIRST_CODE)
+
+#if (MIN_CODE_LEN <= CHAR_BIT)
+#error Code words must be larger than 1 character
+#endif
+
+#if (MAX_CODE_LEN > (2 * CHAR_BIT))
+#error Code words larger than 2 characters require a re-write of GetCodeWord\
+ and PutCodeWord
+#endif
+
+#if ((MAX_CODES - 1) > INT_MAX)
+#error There cannot be more codes than can fit in an integer
+#endif
+
+/***************************************************************************
+*                                  MACROS
+***************************************************************************/
+#define CODE_MS_BITS(BITS)      ((BITS) - CHAR_BIT)
+#define MS_BITS_MASK(BITS)      (UCHAR_MAX << (CHAR_BIT - CODE_MS_BITS(BITS)))
+#define CURRENT_MAX_CODES(BITS)     (1 << (BITS))
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+
+/* dictionary of string the code word is the dictionary index */
+decode_dictionary_t dictionary[DICT_SIZE];
+char currentCodeLen;
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+unsigned char DecodeRecursive(unsigned int code, FILE *fpOut);
+
+/* read encoded data */
+int GetCodeWord(bit_file_t *bfpIn);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/***************************************************************************
+*   Function   : LZWDecodeFile
+*   Description: This routine reads an input file 1 encoded string at a
+*                time and decodes it using the LZW algorithm.
+*   Parameters : inFile - Name of file to decode
+*                outFile - Name of file to write decoded output to
+*   Effects    : File is decoded using the LZW algorithm with CODE_LEN
+*                codes.
+*   Returned   : TRUE for success, otherwise FALSE.
+***************************************************************************/
+int LZWDecodeFile(char *inFile, char *outFile)
+{
+    bit_file_t *bfpIn;                  /* encoded input */
+    FILE *fpOut;                        /* decoded output */
+
+    unsigned int nextCode;              /* value of next code */
+    unsigned int lastCode;              /* last decoded code word */
+    unsigned int code;                  /* code word to decode */
+    unsigned char c;                    /* last decoded character */
+
+    /* open input and output files */
+    if (NULL == (bfpIn = BitFileOpen(inFile, BF_READ)))
+    {
+        perror(inFile);
+        return FALSE;
+    }
+
+    if (NULL == outFile)
+    {
+        fpOut = stdout;
+    }
+    else
+    {
+        if (NULL == (fpOut = fopen(outFile, "wb")))
+        {
+            BitFileClose(bfpIn);
+            perror(outFile);
+            return FALSE;
+        }
+    }
+
+    /* start with 9 bit code words */
+    currentCodeLen = 9;
+
+    /* initialize for decoding */
+    nextCode = FIRST_CODE;  /* code for next (first) string */
+
+    /* first code from file must be a character.  use it for initial values */
+    lastCode = GetCodeWord(bfpIn);
+    c = lastCode;
+    fputc(lastCode, fpOut);
+
+    /* decode rest of file */
+    while ((code = GetCodeWord(bfpIn)) != EOF)
+    {
+
+        /* look for code length increase marker */
+        if (((CURRENT_MAX_CODES(currentCodeLen) - 1) == code) &&
+            (currentCodeLen < MAX_CODE_LEN))
+        {
+            currentCodeLen++;
+            code = GetCodeWord(bfpIn);
+        }
+
+        if (code < nextCode)
+        {
+            /* we have a known code.  decode it */
+            c = DecodeRecursive(code, fpOut);
+        }
+        else
+        {
+            /***************************************************************
+            * We got a code that's not in our dictionary.  This must be due
+            * to the string + char + string + char + string exception.
+            * Build the decoded string using the last character + the
+            * string from the last code.
+            ***************************************************************/
+            unsigned char tmp;
+
+            tmp = c;
+            c = DecodeRecursive(lastCode, fpOut);
+            fputc(tmp, fpOut);
+        }
+
+        /* if room, add new code to the dictionary */
+        if (nextCode < MAX_CODES)
+        {
+            dictionary[nextCode - FIRST_CODE].prefixCode = lastCode;
+            dictionary[nextCode - FIRST_CODE].suffixChar = c;
+            nextCode++;
+        }
+
+        /* save character and code for use in unknown code word case */
+        lastCode = code;
+    }
+
+    fclose(fpOut);
+    BitFileClose(bfpIn);
+    return TRUE;
+}
+
+/***************************************************************************
+*   Function   : DecodeRecursive
+*   Description: This function uses the dictionary to decode a code word
+*                into the string it represents and write it to the output
+*                file.  The string is actually built in reverse order and
+*                recursion is used to write it out in the correct order.
+*   Parameters : code - the code word to decode
+*                fpOut - the file that the decoded code word is written to
+*   Effects    : Decoded code word is written to a file
+*   Returned   : The first character in the decoded string
+***************************************************************************/
+unsigned char DecodeRecursive(unsigned int code, FILE *fpOut)
+{
+    unsigned char c;
+    unsigned char firstChar;
+
+    if (code >= FIRST_CODE)
+    {
+        /* code word is string + c */
+        c = dictionary[code - FIRST_CODE].suffixChar;
+        code = dictionary[code - FIRST_CODE].prefixCode;
+
+        /* evaluate new code word for remaining string */
+        firstChar = DecodeRecursive(code, fpOut);
+    }
+    else
+    {
+        /* code word is just c */
+        c = code;
+        firstChar = code;
+    }
+
+    fputc(c, fpOut);
+    
+    return firstChar;
+}
+
+/***************************************************************************
+*   Function   : GetCodeWord
+*   Description: This function reads and returns a code word from an
+*                encoded file.  In order to deal with endian issue the
+*                code word is read least significant byte followed by the
+*                remaining bits.
+*   Parameters : bfpIn - bit file containing the encoded data
+*   Effects    : code word is read from encoded input
+*   Returned   : The next code word in the encoded file.  EOF if the end
+*                of file has been reached.
+*
+*   NOTE: If the code word contains more than 16 bits, this routine should
+*         be modified to read in all the bytes from least significant to
+*         most significant followed by any left over bits.
+***************************************************************************/
+int GetCodeWord(bit_file_t *bfpIn)
+{
+    unsigned char byte;
+    int code;
+
+    /* get LS character */
+    if ((code = BitFileGetChar(bfpIn)) == EOF)
+    {
+        return EOF;
+    }
+
+
+    /* get remaining bits */
+    if (BitFileGetBits(bfpIn, &byte, CODE_MS_BITS(currentCodeLen)) == EOF)
+    {
+        return EOF;
+    }
+
+    code |= ((int)byte) << CODE_MS_BITS(currentCodeLen);
+
+    return code;
+}
diff --git a/test/compression/lzwencode.c b/test/compression/lzwencode.c
new file mode 100644
index 000000000..efd4a1295
--- /dev/null
+++ b/test/compression/lzwencode.c
@@ -0,0 +1,465 @@
+/***************************************************************************
+*                    Lempel-Ziv-Welch Encoding Functions
+*
+*   File    : lzwencode.c
+*   Purpose : Provides a function for Lempel-Ziv-Welch encoding of file
+*             streams
+*   Author  : Michael Dipperstein
+*   Date    : January 30, 2005
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: lzwencode.c,v 1.3 2007/09/29 01:28:09 michael Exp $
+*   $Log: lzwencode.c,v $
+*   Revision 1.3  2007/09/29 01:28:09  michael
+*   Changes required for LGPL v3.
+*
+*   Revision 1.2  2005/04/21 04:26:57  michael
+*   Encoding dictionary is now built using a binary tree.
+*
+*   Revision 1.1  2005/04/09 03:11:22  michael
+*   Separated encode and decode routines into two different files in order to
+*   make future enhancements easier.
+*
+*   Revision 1.4  2005/03/27 15:56:47  michael
+*   Use variable length code words.
+*   Include new e-mail addres.
+*
+*   Revision 1.3  2005/03/10 14:26:58  michael
+*   User variable names that match discription in web page.
+*
+*   Revision 1.2  2005/03/10 05:44:02  michael
+*   Minimize the size of the dictionary.
+*
+*   Revision 1.1.1.1  2005/02/21 03:35:34  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* LZW: An ANSI C Lempel-Ziv-Welch Encoding/Decoding Routines
+* Copyright (C) 2005, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzw library.
+*
+* The lzw library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzw library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+#include <string.h>
+#include "lzw.h"
+#include "bitfile.h"
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+/* node in dictionary tree */
+typedef struct dict_node_t
+{
+    unsigned int codeWord;      /* code word for this entry */
+    unsigned char suffixChar;   /* last char in encoded string */
+    unsigned int prefixCode;    /* code for remaining chars in string */
+
+    /* pointer to child nodes */
+    struct dict_node_t *left;   /* child with < key */
+    struct dict_node_t *right;  /* child with >= key */
+} dict_node_t;
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+#define MIN_CODE_LEN    9                   /* min # bits in a code word */
+#define MAX_CODE_LEN    15                  /* max # bits in a code word */
+
+#define FIRST_CODE      (1 << CHAR_BIT)     /* value of 1st string code */
+
+#define MAX_CODES       (1 << MAX_CODE_LEN)
+
+#if (MIN_CODE_LEN <= CHAR_BIT)
+#error Code words must be larger than 1 character
+#endif
+
+#if (MAX_CODE_LEN > (2 * CHAR_BIT))
+#error Code words larger than 2 characters require a re-write of GetCodeWord\
+ and PutCodeWord
+#endif
+
+#if ((MAX_CODES - 1) > INT_MAX)
+#error There cannot be more codes than can fit in an integer
+#endif
+
+/***************************************************************************
+*                                  MACROS
+***************************************************************************/
+#define CODE_MS_BITS(BITS)      ((BITS) - CHAR_BIT)
+#define MS_BITS_MASK(BITS)      (UCHAR_MAX << (CHAR_BIT - CODE_MS_BITS(BITS)))
+#define CURRENT_MAX_CODES(BITS)     (1 << (BITS))
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+
+/* dictionary of string codes (encode process uses a hash key) */
+/* XL hack: reuse that of lzwdecode.c */
+extern char currentCodeLen;
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+
+/* dictionary tree node create/free */
+dict_node_t *MakeNode(unsigned int codeWord,
+    unsigned int prefixCode, unsigned char suffixChar);
+void FreeTree(dict_node_t *node);
+
+/* searches tree for matching dictionary entry */
+dict_node_t *FindDictionaryEntry(dict_node_t *root, int prefixCode,
+    unsigned char c);
+
+/* makes key from prefix code and character */
+unsigned int MakeKey(unsigned int prefixCode, unsigned char suffixChar);
+
+/* write encoded data */
+void PutCodeWord(int code, bit_file_t *bfpOut);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/***************************************************************************
+*   Function   : LZWEncodeFile
+*   Description: This routine reads an input file 1 character at a time and
+*                writes out an LZW encoded version of that file.
+*   Parameters : inFile - Name of file to encode
+*                outFile - Name of file to write encoded output to
+*   Effects    : File is encoded using the LZW algorithm with CODE_LEN
+*                codes.
+*   Returned   : TRUE for success, otherwise FALSE.
+***************************************************************************/
+int LZWEncodeFile(char *inFile, char *outFile)
+{
+    FILE *fpIn;                         /* uncoded input */
+    bit_file_t *bfpOut;                 /* encoded output */
+
+    int code;                           /* code for current string */
+    unsigned int nextCode;              /* next available code index */
+    int c;                              /* character to add to string */
+
+    dict_node_t *dictRoot;              /* root of dictionary tree */
+    dict_node_t *node;                  /* node of dictionary tree */
+
+    /* open input and output files */
+    if (NULL == (fpIn = fopen(inFile, "rb")))
+    {
+        perror(inFile);
+        return FALSE;
+    }
+
+    if (NULL == outFile)
+    {
+        bfpOut = MakeBitFile(stdout, BF_WRITE);
+    }
+    else
+    {
+        if (NULL == (bfpOut = BitFileOpen(outFile, BF_WRITE)))
+        {
+            fclose(fpIn);
+            perror(outFile);
+            return FALSE;
+        }
+    }
+
+    /* initialize dictionary as empty */
+    dictRoot = NULL;
+
+    /* start with 9 bit code words */
+    currentCodeLen = 9;
+
+    nextCode = FIRST_CODE;  /* code for next (first) string */
+
+    /* now start the actual encoding process */
+
+    code = fgetc(fpIn);     /* start with code string = first character */
+
+    /* create a tree root from 1st 2 character string */
+    if ((c = fgetc(fpIn)) != EOF)
+    {
+        /* special case for NULL root */
+        dictRoot = MakeNode(nextCode, code, c);
+        nextCode++;
+
+        /* write code for 1st char */
+        PutCodeWord(code, bfpOut);
+
+        /* new code is just 2nd char */
+        code = c;
+    }
+
+    /* now encode normally */
+    while ((c = fgetc(fpIn)) != EOF)
+    {
+        /* look for code + c in the dictionary */
+        node = FindDictionaryEntry(dictRoot, code, c);
+
+        if ((node->prefixCode == code) &&
+            (node->suffixChar == c))
+        {
+            /* code + c is in the dictionary, make it's code the new code */
+            code = node->codeWord;
+        }
+        else
+        {
+            /* code + c is not in the dictionary, add it if there's room */
+            if (nextCode < MAX_CODES)
+            {
+                dict_node_t *tmp;
+
+                tmp = MakeNode(nextCode, code, c);
+                nextCode++;
+
+                if (MakeKey(code, c) <
+                    MakeKey(node->prefixCode, node->suffixChar))
+                {
+                    node->left = tmp;
+                }
+                else
+                {
+                    node->right = tmp;
+                }
+            }
+
+            /* are we using enough bits to write out this code word? */
+            if ((code >= (CURRENT_MAX_CODES(currentCodeLen) - 1)) &&
+                (currentCodeLen < MAX_CODE_LEN))
+            {
+                /* mark need for bigger code word with all ones */
+                PutCodeWord((CURRENT_MAX_CODES(currentCodeLen) - 1), bfpOut);
+                currentCodeLen++;
+            }
+
+            /* write out code for the string before c was added */
+            PutCodeWord(code, bfpOut);
+
+            /* new code is just c */
+            code = c;
+        }
+    }
+
+    /* no more input.  write out last of the code. */
+    PutCodeWord(code, bfpOut);
+
+    fclose(fpIn);
+    BitFileClose(bfpOut);
+
+    /* free the dictionary */
+    if (dictRoot != NULL)
+    {
+        FreeTree(dictRoot);
+    }
+
+    return TRUE;
+}
+
+/***************************************************************************
+*   Function   : MakeKey
+*   Description: This routine creates a simple key from a prefix code and
+*                an appended character.  The key may be used to establish
+*                an order when building/searching a dictionary tree.
+*   Parameters : prefixCode - code for all but the last character of a
+*                             string.
+*                suffixChar - the last character of a string
+*   Effects    : None
+*   Returned   : Key built from string represented as a prefix + char.  Key
+*                format is {ms nibble of c} + prefix + {ls nibble of c}
+***************************************************************************/
+unsigned int MakeKey(unsigned int prefixCode, unsigned char suffixChar)
+{
+    unsigned int key;
+
+    /* position ms nibble */
+    key = suffixChar & 0xF0;
+    key <<= MAX_CODE_LEN;
+
+    /* include prefix code */
+    key |= (prefixCode << 4);
+
+    /* inclulde ls nibble */
+    key |= (suffixChar & 0x0F);
+
+    return key;
+}
+
+/***************************************************************************
+*   Function   : MakeNode
+*   Description: This routine creates and initializes a dictionary entry
+*                for a string and the code word that encodes it.
+*   Parameters : codeWord - code word used to encode the string prefixCode +
+*                           suffixChar
+*                prefixCode - code for all but the last character of a
+*                             string.
+*                suffixChar - the last character of a string
+*   Effects    : Node is allocated for new dictionary entry
+*   Returned   : Pointer to newly allocated node
+***************************************************************************/
+dict_node_t *MakeNode(unsigned int codeWord,
+    unsigned int prefixCode, unsigned char suffixChar)
+{
+    dict_node_t *node;
+
+    node = malloc(sizeof(dict_node_t));
+
+    if (NULL != node)
+    {
+        node->codeWord = codeWord;
+        node->prefixCode = prefixCode;
+        node->suffixChar = suffixChar;
+
+        node->left = NULL;
+        node->right = NULL;
+    }
+    else
+    {
+        perror("allocating dictionary node");
+        exit(EXIT_FAILURE);
+    }
+
+    return node;
+}
+
+/***************************************************************************
+*   Function   : FreeTree
+*   Description: This routine will free all nodes of a tree rooted at the
+*                node past as a parameter.
+*   Parameters : node - root of tree to free
+*   Effects    : frees allocated tree node from initial parameter down.
+*   Returned   : none
+***************************************************************************/
+void FreeTree(dict_node_t *node)
+{
+    /* free left branch */
+    if (node->left != NULL)
+    {
+        FreeTree(node->left);
+    }
+
+    /* free right branch */
+    if (node->right != NULL)
+    {
+        FreeTree(node->right);
+    }
+
+    /* free root */
+    free(node);
+}
+
+/***************************************************************************
+*   Function   : FindDictionaryEntry
+*   Description: This routine searches the dictionary tree for an entry
+*                with a matching string (prefix code + suffix character).
+*                If one isn't found, the parent node for that string is
+*                returned.
+*   Parameters : prefixCode - code for the prefix of string
+*                c - last character in string
+*   Effects    : None
+*   Returned   : If string is in dictionary, pointer to node containing
+*                string, otherwise pointer to suitable parent node.  NULL
+*                is returned for an empty tree.
+***************************************************************************/
+dict_node_t *FindDictionaryEntry(dict_node_t *root, int prefixCode,
+    unsigned char c)
+{
+    unsigned int searchKey, key;
+
+    if (NULL == root)
+    {
+        return NULL;
+    }
+
+    searchKey = MakeKey(prefixCode, c);     /* key of string to find */
+
+    while (1)
+    {
+        /* key of current node */
+        key = MakeKey(root->prefixCode, root->suffixChar);
+
+        if (key == searchKey)
+        {
+            /* current node contains string */
+            return root;
+        }
+        else if (searchKey < key)
+        {
+            if (NULL != root->left)
+            {
+                /* check left branch for string */
+                root = root->left;
+            }
+            else
+            {
+                /* string isn't in tree, it can be added as a left child */
+                return root;
+            }
+        }
+        else
+        {
+            if (NULL != root->right)
+            {
+                /* check right branch for string */
+                root = root->right;
+            }
+            else
+            {
+                /* string isn't in tree, it can be added as a right child */
+                return root;
+            }
+        }
+    }
+}
+
+/***************************************************************************
+*   Function   : PutCodeWord
+*   Description: This function writes a code word from to an encoded file.
+*                In order to deal with endian issue the code word is
+*                written least significant byte followed by the remaining
+*                bits.
+*   Parameters : bfpOut - bit file containing the encoded data
+*                code - code word to add to the encoded data
+*   Effects    : code word is written to the encoded output
+*   Returned   : None
+*
+*   NOTE: If the code word contains more than 16 bits, this routine should
+*         be modified to write out all the bytes from least significant to
+*         most significant followed by any left over bits.
+***************************************************************************/
+void PutCodeWord(int code, bit_file_t *bfpOut)
+{
+    unsigned char byte;
+
+    /* write LS character */
+    byte = code & 0xFF;
+    BitFilePutChar(byte, bfpOut);
+
+    /* write remaining bits */
+    byte = (code >> CODE_MS_BITS(currentCodeLen)) &
+        MS_BITS_MASK(currentCodeLen);
+    BitFilePutBits(bfpOut, &byte, CODE_MS_BITS(currentCodeLen));
+}
diff --git a/test/compression/lzwmain.c b/test/compression/lzwmain.c
new file mode 100644
index 000000000..9f5a8b867
--- /dev/null
+++ b/test/compression/lzwmain.c
@@ -0,0 +1,260 @@
+/***************************************************************************
+*                 Sample Program Using LZW Encoding Library
+*
+*   File    : sample.c
+*   Purpose : Demonstrate usage of LZW encoding library
+*   Author  : Michael Dipperstein
+*   Date    : January 30, 2005
+*
+****************************************************************************
+*   UPDATES
+*
+*   $Id: sample.c,v 1.4 2007/09/29 01:28:27 michael Exp $
+*   $Log: sample.c,v $
+*   Revision 1.4  2007/09/29 01:28:27  michael
+*   Replace getopt with optlist.
+*   Changes required for LGPL v3.
+*
+*   Revision 1.3  2005/03/27 21:10:22  michael
+*   Update e-mail address in copyright notice.
+*
+*   Revision 1.2  2005/02/21 03:52:06  michael
+*   Correct comments in headers.
+*
+*   Revision 1.1.1.1  2005/02/21 03:35:34  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* SAMPLE: Sample usage of Lempel-Ziv-Welch Encoding Library
+* Copyright (C) 2005, 2007 by
+* Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the lzw library.
+*
+* The lzw library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public License as
+* published by the Free Software Foundation; either version 3 of the
+* License, or (at your option) any later version.
+*
+* The lzw library is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "optlist.h"
+#include "lzw.h"
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+char *RemovePath(char *fullPath);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : main
+*   Description: This is the main function for this program, it validates
+*                the command line input and, if valid, it will call
+*                functions to encode or decode a file using the lzw
+*                coding algorithm.
+*   Parameters : argc - number of parameters
+*                argv - parameter list
+*   Effects    : Encodes/Decodes input file
+*   Returned   : EXIT_SUCCESS for success, otherwise EXIT_FAILURE.
+****************************************************************************/
+int main(int argc, char *argv[])
+{
+    option_t *optList, *thisOpt;
+    char *inFile, *outFile; /* name of input & output files */
+    char encode;            /* encode/decode */
+
+    /* initialize data */
+    inFile = NULL;
+    outFile = NULL;
+    encode = TRUE;
+
+    /* parse command line */
+    optList = GetOptList(argc, argv, "cdi:o:h?");
+    thisOpt = optList;
+
+    while (thisOpt != NULL)
+    {
+        switch(thisOpt->option)
+        {
+            case 'c':       /* compression mode */
+                encode = TRUE;
+                break;
+
+            case 'd':       /* decompression mode */
+                encode = FALSE;
+                break;
+
+            case 'i':       /* input file name */
+                if (inFile != NULL)
+                {
+                    fprintf(stderr, "Multiple input files not allowed.\n");
+                    free(inFile);
+
+                    if (outFile != NULL)
+                    {
+                        free(outFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                else if ((inFile =
+                    (char *)malloc(strlen(thisOpt->argument) + 1)) == NULL)
+                {
+                    perror("Memory allocation");
+
+                    if (outFile != NULL)
+                    {
+                        free(outFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                strcpy(inFile, thisOpt->argument);
+                break;
+
+            case 'o':       /* output file name */
+                if (outFile != NULL)
+                {
+                    fprintf(stderr, "Multiple output files not allowed.\n");
+                    free(outFile);
+
+                    if (inFile != NULL)
+                    {
+                        free(inFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+                else if ((outFile =
+                    (char *)malloc(strlen(thisOpt->argument) + 1)) == NULL)
+                {
+                    perror("Memory allocation");
+
+                    if (inFile != NULL)
+                    {
+                        free(inFile);
+                    }
+
+                    FreeOptList(optList);
+                    exit(EXIT_FAILURE);
+                }
+
+                strcpy(outFile, thisOpt->argument);
+                break;
+
+            case 'h':
+            case '?':
+                printf("Usage: %s <options>\n\n", RemovePath(argv[0]));
+                printf("options:\n");
+                printf("  -c : Encode input file to output file.\n");
+                printf("  -d : Decode input file to output file.\n");
+                printf("  -i <filename> : Name of input file.\n");
+                printf("  -o <filename> : Name of output file.\n");
+                printf("  -h | ?  : Print out command line options.\n\n");
+                printf("Default: %s -c\n", RemovePath(argv[0]));
+
+                FreeOptList(optList);
+                return(EXIT_SUCCESS);
+        }
+
+        optList = thisOpt->next;
+        free(thisOpt);
+        thisOpt = optList;
+    }
+
+    /* validate command line */
+    if (inFile == NULL)
+    {
+        fprintf(stderr, "Input file must be provided\n");
+        fprintf(stderr, "Enter \"%s -?\" for help.\n", RemovePath(argv[0]));
+
+        if (outFile != NULL)
+        {
+            free(outFile);
+        }
+
+        exit (EXIT_FAILURE);
+    }
+    else if (outFile == NULL)
+    {
+        fprintf(stderr, "Output file must be provided\n");
+        fprintf(stderr, "Enter \"%s -?\" for help.\n", RemovePath(argv[0]));
+
+        if (inFile != NULL)
+        {
+            free(inFile);
+        }
+
+        exit (EXIT_FAILURE);
+    }
+
+    /* we have valid parameters encode or decode */
+    if (encode)
+    {
+        LZWEncodeFile(inFile, outFile);
+    }
+    else
+    {
+        LZWDecodeFile(inFile, outFile);
+    }
+
+    free(inFile);
+    free(outFile);
+    return EXIT_SUCCESS;
+}
+
+/****************************************************************************
+*   Function   : RemovePath
+*   Description: This is function accepts a pointer to the name of a file
+*                along with path information and returns a pointer to the
+*                character that is not part of the path.
+*   Parameters : fullPath - pointer to an array of characters containing
+*                           a file name and possible path modifiers.
+*   Effects    : None
+*   Returned   : Returns a pointer to the first character after any path
+*                information.
+****************************************************************************/
+char *RemovePath(char *fullPath)
+{
+    int i;
+    char *start, *tmp;                          /* start of file name */
+    const char delim[3] = {'\\', '/', ':'};     /* path deliminators */
+
+    start = fullPath;
+
+    /* find the first character after all file path delimiters */
+    for (i = 0; i < 3; i++)
+    {
+        tmp = strrchr(start, delim[i]);
+
+        if (tmp != NULL)
+        {
+            start = tmp + 1;
+        }
+    }
+
+    return start;
+}
diff --git a/test/compression/optlist.c b/test/compression/optlist.c
new file mode 100644
index 000000000..d1fc013f8
--- /dev/null
+++ b/test/compression/optlist.c
@@ -0,0 +1,228 @@
+/***************************************************************************
+*                       Command Line Option Parser
+*
+*   File    : optlist.c
+*   Purpose : Provide getopt style command line option parsing
+*   Author  : Michael Dipperstein
+*   Date    : August 1, 2007
+*
+****************************************************************************
+*   HISTORY
+*
+*   $Id: optlist.c,v 1.1.1.2 2007/09/04 04:45:42 michael Exp $
+*   $Log: optlist.c,v $
+*   Revision 1.1.1.2  2007/09/04 04:45:42  michael
+*   Added FreeOptList.
+*
+*   Revision 1.1.1.1  2007/08/07 05:01:48  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* OptList: A command line option parsing library
+* Copyright (C) 2007 by Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the OptList library.
+*
+* OptList is free software; you can redistribute it and/or modify it
+* under the terms of the GNU Lesser General Public License as published by
+* the Free Software Foundation; either version 3 of the License, or (at
+* your option) any later version.
+*
+* OptList is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+#include "optlist.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+
+/***************************************************************************
+*                            GLOBAL VARIABLES
+***************************************************************************/
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+option_t *MakeOpt(const char option, char *const argument, const int index);
+
+/***************************************************************************
+*                                FUNCTIONS
+***************************************************************************/
+
+/****************************************************************************
+*   Function   : GetOptList
+*   Description: This function is similar to the POSIX function getopt.  All
+*                options and their corresponding arguments are returned in a
+*                linked list.  This function should only be called once per
+*                an option list and it does not modify argv or argc.
+*   Parameters : argc - the number of command line arguments (including the
+*                       name of the executable)
+*                argv - pointer to the open binary file to write encoded
+*                       output
+*                options - getopt style option list.  A NULL terminated
+*                          string of single character options.  Follow an
+*                          option with a colon to indicate that it requires
+*                          an argument.
+*   Effects    : Creates a link list of command line options and their
+*                arguments.
+*   Returned   : option_t type value where the option and arguement fields
+*                contain the next option symbol and its argument (if any).
+*                The argument field will be set to NULL if the option is
+*                specified as having no arguments or no arguments are found.
+*                The option field will be set to PO_NO_OPT if no more
+*                options are found.
+*
+*   NOTE: The caller is responsible for freeing up the option list when it
+*         is no longer needed.
+****************************************************************************/
+option_t *GetOptList(const int argc, char *const argv[], char *const options)
+{
+    int nextArg;
+    option_t *head, *tail;
+    int optIndex;
+
+    /* start with first argument and nothing found */
+    nextArg = 1;
+    head = NULL;
+    tail = NULL;
+
+    /* loop through all of the command line arguments */
+    while (nextArg < argc)
+    {
+        if ((strlen(argv[nextArg]) > 1) && ('-' == argv[nextArg][0]))
+        {
+            /* possible option */
+            optIndex = 0;
+
+            /* attempt to find a matching option */
+            while ((options[optIndex] != '\0') &&
+                (options[optIndex] != argv[nextArg][1]))
+            {
+                do
+                {
+                    optIndex++;
+                }
+                while ((options[optIndex] != '\0') &&
+                    (':' == options[optIndex]));
+            }
+
+            if (options[optIndex] == argv[nextArg][1])
+            {
+                /* we found the matching option */
+                if (NULL == head)
+                {
+                    head = MakeOpt(options[optIndex], NULL, OL_NOINDEX);
+                    tail = head;
+                }
+                else
+                {
+                    tail->next = MakeOpt(options[optIndex], NULL, OL_NOINDEX);
+                    tail = tail->next;
+                }
+
+                if (':' == options[optIndex + 1])
+                {
+                    /* the option found should have a text arguement */
+                    if (strlen(argv[nextArg]) > 2)
+                    {
+                        /* no space between argument and option */
+                        tail->argument = &(argv[nextArg][2]);
+                        tail->argIndex = nextArg;
+                    }
+                    else if (nextArg < argc)
+                    {
+                        /* there must be space between the argument option */
+                        nextArg++;
+                        tail->argument = argv[nextArg];
+                        tail->argIndex = nextArg;
+                    }
+                }
+            }
+        }
+
+        nextArg++;
+    }
+
+    return head;
+}
+
+/****************************************************************************
+*   Function   : MakeOpt
+*   Description: This function uses malloc to allocate space for an option_t
+*                type structure and initailizes the structure with the
+*                values passed as a parameter.
+*   Parameters : option - this option character
+*                argument - pointer string containg the argument for option.
+*                           Use NULL for no argument
+*                index - argv[index] contains argument us OL_NOINDEX for
+*                        no argument
+*   Effects    : A new option_t type variable is created on the heap.
+*   Returned   : Pointer to newly created and initialized option_t type
+*                structure.  NULL if space for structure can't be allocated.
+****************************************************************************/
+option_t *MakeOpt(const char option, char *const argument, const int index)
+{
+    option_t *opt;
+
+    opt = malloc(sizeof(option_t));
+
+    if (opt != NULL)
+    {
+        opt->option = option;
+        opt->argument = argument;
+        opt->argIndex = index;
+        opt->next = NULL;
+    }
+    else
+    {
+        perror("Failed to Allocate option_t");
+    }
+
+    return opt;
+}
+
+/****************************************************************************
+*   Function   : FreeOptList
+*   Description: This function will free all the elements in an option_t
+*                type linked list starting from the node passed as a
+*                parameter.
+*   Parameters : list - head of linked list to be freed
+*   Effects    : All elements of the linked list pointed to by list will
+*                be freed and list will be set to NULL.
+*   Returned   : None
+****************************************************************************/
+void FreeOptList(option_t *list)
+{
+    option_t *head, *next;
+
+    head = list;
+    list = NULL;
+
+    while (head != NULL)
+    {
+        next = head->next;
+        free(head);
+        head = next;
+    }
+
+    return;
+}
diff --git a/test/compression/optlist.h b/test/compression/optlist.h
new file mode 100644
index 000000000..046126336
--- /dev/null
+++ b/test/compression/optlist.h
@@ -0,0 +1,74 @@
+/***************************************************************************
+*                       Command Line Option Parser
+*
+*   File    : optlist.h
+*   Purpose : Header for getopt style command line option parsing
+*   Author  : Michael Dipperstein
+*   Date    : August 1, 2007
+*
+****************************************************************************
+*   HISTORY
+*
+*   $Id: optlist.h,v 1.1.1.2 2007/09/04 04:45:42 michael Exp $
+*   $Log: optlist.h,v $
+*   Revision 1.1.1.2  2007/09/04 04:45:42  michael
+*   Added FreeOptList.
+*
+*   Revision 1.1.1.1  2007/08/07 05:01:48  michael
+*   Initial Release
+*
+****************************************************************************
+*
+* OptList: A command line option parsing library
+* Copyright (C) 2007 by Michael Dipperstein (mdipper@alumni.engr.ucsb.edu)
+*
+* This file is part of the OptList library.
+*
+* OptList is free software; you can redistribute it and/or modify it
+* under the terms of the GNU Lesser General Public License as published by
+* the Free Software Foundation; either version 3 of the License, or (at
+* your option) any later version.
+*
+* OptList is distributed in the hope that it will be useful, but
+* WITHOUT ANY WARRANTY; without even the implied warranty of
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser
+* General Public License for more details.
+*
+* You should have received a copy of the GNU Lesser General Public License
+* along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*
+***************************************************************************/
+#ifndef OPTLIST_H
+#define OPTLIST_H
+
+/***************************************************************************
+*                             INCLUDED FILES
+***************************************************************************/
+
+/***************************************************************************
+*                                 MACROS
+***************************************************************************/
+
+/***************************************************************************
+*                                CONSTANTS
+***************************************************************************/
+#define    OL_NOINDEX    -1        /* this option has no arguement */
+
+/***************************************************************************
+*                            TYPE DEFINITIONS
+***************************************************************************/
+typedef struct option_t
+{
+    char option;
+    char *argument;
+    int argIndex;
+    struct option_t *next;
+} option_t;
+
+/***************************************************************************
+*                               PROTOTYPES
+***************************************************************************/
+option_t *GetOptList(int argc, char *const argv[], char *const options);
+void FreeOptList(option_t *list);
+
+#endif  /* ndef OPTLIST_H */
diff --git a/test/raytracer/.depend b/test/raytracer/.depend
new file mode 100644
index 000000000..c31fbcdc1
--- /dev/null
+++ b/test/raytracer/.depend
@@ -0,0 +1,20 @@
+arrays.o: arrays.c config.h arrays.h
+eval.o: eval.c config.h arrays.h gml.h point.h vector.h eval.h object.h \
+  light.h render.h
+gmllexer.o: gmllexer.c config.h gmllexer.h gml.h
+gmlparser.o: gmlparser.c config.h arrays.h gml.h gmllexer.h gmlparser.h
+intersect.o: intersect.c config.h arrays.h point.h vector.h eval.h \
+  object.h matrix.h intersect.h
+light.o: light.c config.h point.h vector.h eval.h object.h intersect.h \
+  light.h
+main.o: main.c config.h arrays.h gmllexer.h gmlparser.h eval.h
+matrix.o: matrix.c config.h point.h vector.h matrix.h
+memory.o: memory.c config.h
+object.o: object.c config.h arrays.h point.h vector.h eval.h object.h \
+  matrix.h
+render.o: render.c config.h point.h vector.h matrix.h eval.h object.h \
+  light.h intersect.h surface.h simplify.h render.h
+simplify.o: simplify.c config.h point.h vector.h arrays.h matrix.h eval.h \
+  object.h simplify.h
+surface.o: surface.c config.h point.h vector.h eval.h object.h surface.h
+vector.o: vector.c config.h point.h vector.h
diff --git a/test/raytracer/Makefile b/test/raytracer/Makefile
new file mode 100644
index 000000000..cd442c586
--- /dev/null
+++ b/test/raytracer/Makefile
@@ -0,0 +1,43 @@
+CC=../../ccomp 
+CFLAGS=-U__GNUC__ -stdlib ../../runtime -dclight -dasm
+LIBS=
+
+OBJS=memory.o gmllexer.o gmlparser.o eval.o \
+  arrays.o vector.o matrix.o object.o intersect.o surface.o light.o \
+  simplify.o render.o main.o
+
+render: $(OBJS)
+	$(CC) $(CFLAGS) -o render $(OBJS) $(LIBS)
+
+clean:
+	rm -f *.o *.light.c *.s *.ppm render
+
+include .depend
+
+depend:
+	gcc -MM *.c > .depend
+
+gcc0:
+	$(MAKE) clean
+	$(MAKE) CC=gcc CFLAGS="-O0 -Wall"
+	mv render render.gcc0
+	$(MAKE) clean
+
+gcc1:
+	$(MAKE) clean
+	$(MAKE) CC=gcc CFLAGS="-O1 -Wall"
+	mv render render.gcc1
+	$(MAKE) clean
+
+gcc2:
+	$(MAKE) clean
+	$(MAKE) CC=gcc CFLAGS="-O2 -Wall"
+	mv render render.gcc2
+	$(MAKE) clean
+
+gcc3:
+	$(MAKE) clean
+	$(MAKE) CC=gcc CFLAGS="-O3 -Wall"
+	mv render render.gcc3
+	$(MAKE) clean
+
diff --git a/test/raytracer/arrays.c b/test/raytracer/arrays.c
new file mode 100644
index 000000000..43508ee66
--- /dev/null
+++ b/test/raytracer/arrays.c
@@ -0,0 +1,30 @@
+#include "config.h"
+#include "arrays.h"
+
+struct array * alloc_array(int eltsize, int initialsize)
+{
+  struct array * a = arena_alloc(sizeof(struct array));
+  a->size = 0;
+  a->capacity = initialsize;
+  a->data = arena_alloc(eltsize * initialsize);
+  return a;
+}
+
+void grow_array(int eltsize, struct array * arr)
+{
+  void * newdata = arena_alloc(arr->capacity * 2 * eltsize);
+  memcpy(newdata, arr->data, arr->size * eltsize);
+  arr->data = newdata;
+  arr->capacity *= 2;
+}
+
+struct array * copy_array(int eltsize, struct array * arr, int extrasize)
+{
+  struct array * a = arena_alloc(sizeof(struct array));
+  a->size = arr->size;
+  a->capacity = arr->size + extrasize;
+  a->data = arena_alloc(eltsize * a->capacity);
+  memcpy(a->data, arr->data, eltsize * a->size);
+  return a;
+}
+
diff --git a/test/raytracer/arrays.h b/test/raytracer/arrays.h
new file mode 100644
index 000000000..47634f62d
--- /dev/null
+++ b/test/raytracer/arrays.h
@@ -0,0 +1,31 @@
+/* Resizable arrays */
+
+struct array {
+  int size;                     /* Number of elements in the array */
+  int capacity;                 /* Max number of elements before resizing */
+  void * data;                  /* Actual data -- variable sized */
+};
+
+struct array * alloc_array(int eltsize, int initialsize);
+
+void grow_array(int eltsize, struct array * arr);
+
+struct array * copy_array(int eltsize, struct array * arr, int extrasize);
+
+#define new_array(ty,sz) alloc_array(sizeof(ty), sz)
+
+#define data_array(ty,arr) ((ty *) (arr)->data)
+
+#define get_array(ty,arr,idx) (data_array(ty,arr)[idx])
+
+#define get_array_large(dst,ty,arr,idx) \
+  ASSIGN(dst, data_array(ty,arr)[idx])
+
+#define set_array(ty,arr,idx,newval) (data_array(ty,arr)[idx] = (newval))
+
+#define set_array_large(ty,arr,idx,newval) \
+  ASSIGN(data_array(ty,arr)[idx], newval)
+
+#define extend_array(ty,arr) \
+  if ((arr)->size >= (arr)->capacity) grow_array(sizeof(ty), (arr)); \
+  (arr)->size++
diff --git a/test/raytracer/config.h b/test/raytracer/config.h
new file mode 100644
index 000000000..ed5b8478c
--- /dev/null
+++ b/test/raytracer/config.h
@@ -0,0 +1,27 @@
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+
+#ifdef SINGLE_PRECISION
+typedef float flt;
+#else
+typedef double flt;
+#endif
+
+#if 0
+extern void abord(void);
+
+#define assert(cond) \
+  if (!(cond)) { fprintf(stderr, "%s:%u: failed assertion\n", __FILE__, __LINE__); abort(); }
+
+#endif
+
+#define ASSIGN(lv,rv) memcpy(&(lv), &(rv), sizeof(rv))
+
+void arena_init(void);
+void arena_clear(void);
+void * arena_alloc(int size);
+
diff --git a/test/raytracer/eval.c b/test/raytracer/eval.c
new file mode 100644
index 000000000..4d3f3dcc4
--- /dev/null
+++ b/test/raytracer/eval.c
@@ -0,0 +1,672 @@
+/* Evaluator for GML */
+
+#include "config.h"
+#include "arrays.h"
+#include "gml.h"
+#include "point.h"
+#include "vector.h"
+#include "eval.h"
+#include "object.h"
+#include "light.h"
+#include "render.h"
+
+struct value {
+  enum { B, I, R, S, Arr, Clos, Point, Obj, Light } tag;
+  union {
+    int i;                      /* B, I */
+    flt r;                   /* R */
+    char * s;                   /* S */
+    struct array * arr;         /* Arr */
+    struct closure clos;        /* Clos */
+    struct point * point;       /* Point */
+    struct object * obj;        /* Obj */
+    struct light * light;       /* Light */
+  } u;
+};
+
+struct binding {
+  char * name;
+  int mutable;
+  struct value val;
+};
+
+/* Lookup an identifier in an environment */
+
+static int lookup(struct array * env, char * name, /*out*/ struct value * res)
+{
+  int i;
+  for (i = env->size - 1; i >= 0; i--) {
+    struct binding * b = &get_array(struct binding, env, i);
+    if (name == b->name) {
+      ASSIGN(*res, b->val);
+      return 1;
+    }
+  }
+  return 0;
+}
+
+/* Assign an identifier in an environment */
+
+static void assign(struct array * env, char * name, struct value * newval)
+{
+  int i;
+  struct binding * b;
+  for (i = env->size - 1; i >= 0; i--) {
+    b = &get_array(struct binding, env, i);
+    if (! b->mutable) break;
+    if (name == b->name) {
+      ASSIGN(b->val, *newval);
+      return;
+    }
+  }
+  extend_array(struct binding, env);
+  b = &get_array(struct binding, env, env->size - 1);
+  b->name = name;
+  b->mutable = 1;
+  ASSIGN(b->val, *newval);
+}
+
+/* Take an immutable copy of an environment */
+
+static struct array * snapshot_env(struct array * env)
+{
+  int i;
+  struct array * nenv = copy_array(sizeof(struct binding), env, 10);
+  for (i = 0; i < nenv->size; i++)
+    get_array(struct binding, nenv, i).mutable = 0;
+  return nenv;
+}
+
+/* Utility math functions */
+
+static inline flt degrees_to_radians(flt d)
+{ return d * (M_PI / 180.0); }
+
+static inline flt radians_to_degrees(flt d)
+{ return d * (180.0 / M_PI); }
+
+static inline flt deg_cos(flt a)
+{ return cos(degrees_to_radians(a)); }
+
+static inline flt deg_sin(flt a)
+{ return sin(degrees_to_radians(a)); }
+
+static inline flt deg_acos(flt a)
+{ return radians_to_degrees(acos(a)); }
+
+static inline flt deg_asin(flt a)
+{ return radians_to_degrees(asin(a)); }
+
+static inline flt clampf(flt r)
+{
+  if (r < 0.0) return 0.0;
+  if (r > 1.0) return 1.0;
+  return r;
+}
+
+/* For error printing */
+
+static char * operator_names [] = {
+  "<Identifier>", "<Binder>", "<Boolean>", "<Integer>", "<Real>",
+  "<String>", "<Array>", "<Function>", "acos", "addi", "addf", "apply",
+  "asin", "clampf", "cone", "cos", "cube", "cylinder", "difference",
+  "divi", "divf", "eqi", "eqf", "floor", "frac", "get", "getx", "gety",
+  "getz", "if", "intersect", "length", "lessi", "lessf", "light",
+  "modi", "muli", "mulf", "negi", "negf", "plane", "point",
+  "pointlight", "real", "render", "rotatex", "rotatey", "rotatez",
+  "scale", "sin", "sphere", "spotlight", "sqrt", "subi", "subf",
+  "translate", "union", "uscale", "print",
+};
+
+/* Convert a GML array of lights into a C array of lights */
+
+static struct light ** light_array(struct array * arr)
+{
+  int sz = arr->size;
+  struct light ** l = arena_alloc((sz + 1) * sizeof(struct light *));
+  int i;
+  for (i = 0; i < sz; i++) {
+    struct value * v = &get_array(struct value, arr, i);
+    if (v->tag != Light) {
+      fprintf(stderr, "Light expected in array argument to `render'\n");
+      exit(2);
+    }
+    l[i] = v->u.light;
+  }
+  l[sz] = NULL;
+  return l;
+}
+
+/* Pretty-print a value */
+
+static void print_value(struct value * s)
+{
+  switch (s->tag) {
+  case B:
+    printf("%s\n", s->u.i ? "true" : "false"); break;
+  case I:
+    printf("%d\n", s->u.i); break;
+  case R:
+    printf("%e\n", s->u.r); break;
+  case S:
+    printf("\"%s\"\n", s->u.s); break;
+  case Arr:
+    printf("array\n"); break;
+  case Clos:
+    printf("closure\n"); break;
+  case Point:
+    printf("point %e %e %e\n", 
+           s->u.point->x, s->u.point->y, s->u.point->z);
+    break;
+  case Obj:
+    printf("object\n"); break;
+  case Light:
+    printf("light\n"); break;
+  }
+}
+
+/* The evaluation stack */
+
+#define MAIN_STACK_SIZE 10000
+#define SURFACE_STACK_SIZE 1000
+
+static struct value main_stack[MAIN_STACK_SIZE];
+static struct value surface_stack[SURFACE_STACK_SIZE];
+
+/* Error handling functions */
+#define ERRORFUN(name, msg) \
+  static void name(void) { fprintf(stderr, msg); exit(2); }
+
+ERRORFUN(stack_overflow, "Stack overflow\n")
+ERRORFUN(stack_underflow, "Stack underflow\n")
+ERRORFUN(type_error, "Type error\n")
+ERRORFUN(division_by_zero, "Division by zero\n")
+ERRORFUN(bound_error, "Out-of-bound array access\n")
+ERRORFUN(negative_sqrt, "Square root of negative number\n")
+
+/* Macros for stack checking and type checking */
+
+#define push() if (--sp < bos) stack_overflow()
+#define can_pop(n) if (sp + (n) > tos) stack_underflow()
+#define check(n,ty) if (sp[n].tag != ty) type_error()
+
+/* Execute the given token list in the given environment */
+
+static struct value * execute_list(struct array * code,
+                                   struct array * env,
+                                   struct value * bos,
+                                   struct value * tos,
+                                   struct value * sp)
+{
+  int i;
+  struct tok * t;
+
+  for (i = 0; i < code->size; i++) {
+    t = &get_array(struct tok, code, i);
+    switch (t->tag) {
+    case Identifier:
+      push();
+      if (! lookup(env, t->u.s, sp)) {
+        fprintf(stderr, "Unbound identifier %s\n", t->u.s);
+        exit(2);
+      }
+      break;
+    case Binder:
+      can_pop(1);
+      assign(env, t->u.s, sp);
+      sp++;
+      break;
+    case Boolean:
+      push();
+      sp[0].tag = B; sp[0].u.i = t->u.i;
+      break;
+    case Integer:
+      push();
+      sp[0].tag = I; sp[0].u.i = t->u.i;
+      break;
+    case Real:
+      push();
+      sp[0].tag = R; sp[0].u.r = t->u.d;
+      break;
+    case String:
+      push();
+      sp[0].tag = S; sp[0].u.s = t->u.s;
+      break;
+    case Array: {
+      struct value * esp = execute_list(t->u.a, env, bos, sp, sp);
+      int sz = sp - esp;
+      struct array * a = new_array(struct value, sz);
+      int j;
+      a->size = sz;
+      for (j = 0; j < sz; j++) set_array_large(struct value, a, j, sp[-1-j]);
+      push();
+      sp[0].tag = Arr; sp[0].u.arr = a;
+      break;
+    }
+    case Function:
+      push();
+      sp[0].tag = Clos;
+      sp[0].u.clos.code = t->u.a;
+      sp[0].u.clos.env = snapshot_env(env);
+      break;
+    case Op_acos:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = deg_acos(sp[0].u.r);
+      break;
+    case Op_addi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      sp[1].u.i = sp[1].u.i + sp[0].u.i;
+      sp += 1;
+      break;
+    case Op_addf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].u.r = sp[1].u.r + sp[0].u.r;
+      sp += 1;
+      break;
+    case Op_apply:
+      can_pop(1);
+      check(0, Clos);
+      sp = execute_list(sp[0].u.clos.code, sp[0].u.clos.env, bos, tos, sp + 1);
+      break;
+    case Op_asin:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = deg_asin(sp[0].u.r);
+      break;
+    case Op_clampf:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = clampf(sp[0].u.r);
+      break;
+    case Op_cone:
+      can_pop(1);
+      check(0, Clos);
+      sp[0].tag = Obj;
+      sp[0].u.obj = cone(&sp[0].u.clos);
+      break;
+    case Op_cos:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = deg_cos(sp[0].u.r);
+      break;
+    case Op_cube:
+      can_pop(1);
+      check(0, Clos);
+      sp[0].tag = Obj;
+      sp[0].u.obj = cube(&sp[0].u.clos);
+      break;
+    case Op_cylinder:
+      can_pop(1);
+      check(0, Clos);
+      sp[0].tag = Obj;
+      sp[0].u.obj = cylinder(&sp[0].u.clos);
+      break;
+    case Op_difference:
+      can_pop(2);
+      check(1, Obj);
+      check(0, Obj);
+      sp[1].u.obj = odifference(sp[1].u.obj, sp[0].u.obj);
+      sp += 1;
+      break;
+    case Op_divi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      if (sp[0].u.i == 0) division_by_zero();
+      sp[1].u.i = sp[1].u.i / sp[0].u.i;
+      sp += 1;
+      break;
+    case Op_divf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].u.r = sp[1].u.r / sp[0].u.r;
+      sp += 1;
+      break;
+    case Op_eqi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      sp[1].tag = B;
+      sp[1].u.i = (sp[1].u.i == sp[0].u.i);
+      sp += 1;
+      break;
+    case Op_eqf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].tag = B;
+      sp[1].u.i = (sp[1].u.r == sp[0].u.r);
+      sp += 1;
+      break;
+    case Op_floor:
+      can_pop(1);
+      check(0, R);
+      sp[0].tag = I;
+      sp[0].u.i = (int) floor(sp[0].u.r);
+      break;
+    case Op_frac: {
+      double rem;
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = modf(sp[0].u.r, &rem);
+      break;
+    }
+    case Op_get: {
+      struct array * a;
+      int idx;
+      can_pop(2);
+      check(1, Arr);
+      check(0, I);
+      a = sp[1].u.arr;
+      idx = sp[0].u.i;
+      if (idx < 0 || idx >= a->size) bound_error();
+      get_array_large(sp[1], struct value, a, idx);
+      sp++;
+      break;
+    }
+    case Op_getx:
+      can_pop(1);
+      check(0, Point);
+      sp[0].tag = R;
+      sp[0].u.r = sp[0].u.point->x;
+      break;
+    case Op_gety:
+      can_pop(1);
+      check(0, Point);
+      sp[0].tag = R;
+      sp[0].u.r = sp[0].u.point->y;
+      break;
+    case Op_getz:
+      can_pop(1);
+      check(0, Point);
+      sp[0].tag = R;
+      sp[0].u.r = sp[0].u.point->z;
+      break;
+    case Op_if:
+      can_pop(3);
+      check(2, B);
+      check(1, Clos);
+      check(0, Clos);
+      if (sp[2].u.i)
+        sp = execute_list(sp[1].u.clos.code, sp[1].u.clos.env, bos, tos, sp + 3);
+      else
+        sp = execute_list(sp[0].u.clos.code, sp[0].u.clos.env, bos, tos, sp + 3);
+      break;
+    case Op_intersect:
+      can_pop(2);
+      check(1, Obj);
+      check(0, Obj);
+      sp[1].u.obj = ointersect(sp[1].u.obj, sp[0].u.obj);
+      sp += 1;
+      break;
+    case Op_length:
+      can_pop(1);
+      check(0, Arr);
+      sp[0].tag = I;
+      sp[0].u.i = sp[0].u.arr->size;
+      break;
+    case Op_lessi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      sp[1].tag = B;
+      sp[1].u.i = (sp[1].u.i < sp[0].u.i);
+      sp += 1;
+      break;
+    case Op_lessf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].tag = B;
+      sp[1].u.i = (sp[1].u.r < sp[0].u.r);
+      sp += 1;
+      break;
+    case Op_light:
+      can_pop(2);
+      check(1, Point);
+      check(0, Point);
+      sp[1].tag = Light;
+      sp[1].u.light = dirlight(sp[1].u.point, sp[0].u.point);
+      sp += 1;
+      break;
+    case Op_modi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      if (sp[0].u.i == 0) division_by_zero();
+      sp[1].u.i = sp[1].u.i % sp[0].u.i;
+      sp += 1;
+      break;
+    case Op_muli:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      sp[1].u.i = sp[1].u.i * sp[0].u.i;
+      sp += 1;
+      break;
+    case Op_mulf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].u.r = sp[1].u.r * sp[0].u.r;
+      sp += 1;
+      break;
+    case Op_negi:
+      can_pop(1);
+      check(0, I);
+      sp[0].u.i = - sp[0].u.i;
+      break;
+    case Op_negf:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = - sp[0].u.r;
+      break;
+    case Op_plane:
+      can_pop(1);
+      check(0, Clos);
+      sp[0].tag = Obj;
+      sp[0].u.obj = plane(&sp[0].u.clos);
+      break;
+    case Op_point: {
+      struct point * p;
+      can_pop(3);
+      check(2, R);
+      check(1, R);
+      check(0, R);
+      p = arena_alloc(sizeof(struct point));
+      p->x = sp[2].u.r;
+      p->y = sp[1].u.r;
+      p->z = sp[0].u.r;
+      sp += 2;
+      sp[0].tag = Point;
+      sp[0].u.point = p;
+      break;
+    }
+    case Op_pointlight:
+      can_pop(2);
+      check(1, Point);
+      check(0, Point);
+      sp[1].tag = Light;
+      sp[1].u.light = pointlight(sp[1].u.point, sp[0].u.point);
+      sp += 1;
+      break;
+    case Op_real:
+      can_pop(1);
+      check(0, I);
+      sp[0].tag = R;
+      sp[0].u.r = (flt) sp[0].u.i;
+      break;
+    case Op_render:
+      can_pop(8);
+      check(0, S);
+      check(1, I);
+      check(2, I);
+      check(3, R);
+      check(4, I);
+      check(5, Obj);
+      check(6, Arr);
+      check(7, Point);
+      render(sp[7].u.point,
+             sp[6].u.arr->size,
+             light_array(sp[6].u.arr),
+             sp[5].u.obj,
+             sp[4].u.i,
+             degrees_to_radians(sp[3].u.r),
+             sp[2].u.i,
+             sp[1].u.i,
+             sp[0].u.s);
+      sp += 8;
+      break;
+    case Op_rotatex:
+      can_pop(2);
+      check(1, Obj);
+      check(0, R);
+      sp[1].u.obj = orotatex(sp[1].u.obj, degrees_to_radians(sp[0].u.r));
+      sp += 1;
+      break;
+    case Op_rotatey:
+      can_pop(2);
+      check(1, Obj);
+      check(0, R);
+      sp[1].u.obj = orotatey(sp[1].u.obj, degrees_to_radians(sp[0].u.r));
+      sp += 1;
+      break;
+    case Op_rotatez:
+      can_pop(2);
+      check(1, Obj);
+      check(0, R);
+      sp[1].u.obj = orotatez(sp[1].u.obj, degrees_to_radians(sp[0].u.r));
+      sp += 1;
+      break;
+    case Op_scale:
+      can_pop(4);
+      check(3, Obj);
+      check(2, R);
+      check(1, R);
+      check(0, R);
+      sp[3].u.obj = oscale(sp[3].u.obj, sp[2].u.r, sp[1].u.r, sp[0].u.r);
+      sp += 3;
+      break;
+    case Op_sin:
+      can_pop(1);
+      check(0, R);
+      sp[0].u.r = deg_sin(sp[0].u.r);
+      break;
+    case Op_sphere:
+      can_pop(1);
+      check(0, Clos);
+      sp[0].tag = Obj;
+      sp[0].u.obj = sphere(&sp[0].u.clos);
+      break;
+    case Op_spotlight:
+      can_pop(5);
+      check(4, Point);
+      check(3, Point);
+      check(2, Point);
+      check(1, R);
+      check(0, R);
+      sp[4].tag = Light;
+      sp[4].u.light = spotlight(sp[4].u.point, sp[3].u.point, sp[2].u.point, 
+                                degrees_to_radians(sp[1].u.r), sp[0].u.r);
+      sp += 4;
+      break;
+    case Op_sqrt:
+      can_pop(1);
+      check(0, R);
+      if (sp[0].u.r < 0) negative_sqrt();
+      sp[0].u.r = sqrt(sp[0].u.r);
+      break;
+    case Op_subi:
+      can_pop(2);
+      check(1, I);
+      check(0, I);
+      sp[1].u.i = sp[1].u.i - sp[0].u.i;
+      sp += 1;
+      break;
+    case Op_subf:
+      can_pop(2);
+      check(1, R);
+      check(0, R);
+      sp[1].u.r = sp[1].u.r - sp[0].u.r;
+      sp += 1;
+      break;
+    case Op_translate:
+      can_pop(4);
+      check(3, Obj);
+      check(2, R);
+      check(1, R);
+      check(0, R);
+      sp[3].u.obj = otranslate(sp[3].u.obj, sp[2].u.r, sp[1].u.r, sp[0].u.r);
+      sp += 3;
+      break;
+    case Op_union:
+      can_pop(2);
+      check(1, Obj);
+      check(0, Obj);
+      sp[1].u.obj = ounion(sp[1].u.obj, sp[0].u.obj);
+      sp += 1;
+      break;
+    case Op_uscale:
+      can_pop(2);
+      check(1, Obj);
+      check(0, R);
+      sp[1].u.obj = ouscale(sp[1].u.obj, sp[0].u.r);
+      sp += 1;
+      break;
+    case Op_print:
+      can_pop(1);
+      print_value(sp);
+      sp += 1;
+      break;
+    }
+  }
+  return sp;
+}
+
+/* Evaluate a surface function */
+
+void surface_function(struct closure * clos, int face, flt u, flt v,
+                      /*out*/ struct surface_characteristics * sc)
+{
+  struct value * sp;
+  sp = surface_stack + SURFACE_STACK_SIZE - 3;
+  sp[2].tag = I;
+  sp[2].u.i = face;
+  sp[1].tag = R;
+  sp[1].u.i = u;
+  sp[0].tag = R;
+  sp[0].u.i = v;
+  sp =
+    execute_list(clos->code, clos->env, surface_stack,
+                 surface_stack + SURFACE_STACK_SIZE, sp);
+  if (sp != surface_stack + SURFACE_STACK_SIZE - 4 ||
+      sp[0].tag != R ||
+      sp[1].tag != R ||
+      sp[2].tag != R ||
+      sp[3].tag != Point) {
+    fprintf(stderr, "Wrong result for surface function\n");
+    exit(2);
+  }
+  sc->x = sp[3].u.point->x;
+  sc->y = sp[3].u.point->y;
+  sc->z = sp[3].u.point->z;
+  sc->kd = sp[2].u.r;
+  sc->ks = sp[1].u.r;
+  sc->phong = sp[0].u.r;
+}
+
+/* Execute the main program */
+
+void execute_program(struct array * toklist)
+{
+  execute_list(toklist, new_array(struct binding, 50),
+               main_stack,
+               main_stack + MAIN_STACK_SIZE,
+               main_stack + MAIN_STACK_SIZE);
+}
diff --git a/test/raytracer/eval.h b/test/raytracer/eval.h
new file mode 100644
index 000000000..610b5d967
--- /dev/null
+++ b/test/raytracer/eval.h
@@ -0,0 +1,18 @@
+/* Evaluator for GML */
+
+struct closure {
+  struct array * code;
+  struct array * env;
+};
+
+struct surface_characteristics {
+  flt x, y, z;
+  flt kd;
+  flt ks;
+  flt phong;
+};
+
+void execute_program(struct array * toklist);
+
+void surface_function(struct closure * clos, int face, flt u, flt v,
+                      /*out*/ struct surface_characteristics * sc);
diff --git a/test/raytracer/gml.h b/test/raytracer/gml.h
new file mode 100644
index 000000000..20ca3e94d
--- /dev/null
+++ b/test/raytracer/gml.h
@@ -0,0 +1,73 @@
+/* The GML abstract syntax tree */
+
+enum operation {
+    Identifier,
+    Binder,
+    Boolean,
+    Integer,
+    Real,
+    String,
+    Array,
+    Function,
+    Op_acos,
+    Op_addi,
+    Op_addf,
+    Op_apply,
+    Op_asin,
+    Op_clampf,
+    Op_cone,
+    Op_cos,
+    Op_cube,
+    Op_cylinder,
+    Op_difference,
+    Op_divi,
+    Op_divf,
+    Op_eqi,
+    Op_eqf,
+    Op_floor,
+    Op_frac,
+    Op_get,
+    Op_getx,
+    Op_gety,
+    Op_getz,
+    Op_if,
+    Op_intersect,
+    Op_length,
+    Op_lessi,
+    Op_lessf,
+    Op_light,
+    Op_modi,
+    Op_muli,
+    Op_mulf,
+    Op_negi,
+    Op_negf,
+    Op_plane,
+    Op_point,
+    Op_pointlight,
+    Op_real,
+    Op_render,
+    Op_rotatex,
+    Op_rotatey,
+    Op_rotatez,
+    Op_scale,
+    Op_sin,
+    Op_sphere,
+    Op_spotlight,
+    Op_sqrt,
+    Op_subi,
+    Op_subf,
+    Op_translate,
+    Op_union,
+    Op_uscale,
+    Op_print
+};
+
+struct tok {
+  enum operation tag;
+  union {
+    char * s;                   /* Identifier, Binder, String */
+    int i;                      /* Boolean, Integer */
+    flt d;                      /* Real */
+    struct array * a;           /* Array, Function */
+  } u;
+};
diff --git a/test/raytracer/gmllexer.c b/test/raytracer/gmllexer.c
new file mode 100644
index 000000000..332126b29
--- /dev/null
+++ b/test/raytracer/gmllexer.c
@@ -0,0 +1,297 @@
+/* Lexer for GML */
+
+#include "config.h"
+#include "gmllexer.h"
+#include "gml.h"
+
+#define isdigit(c) (c >= '0' && c <= '9')
+#define isalpha(c) ((c >= 'A' && c <= 'Z') || (c >= 'a' && c <= 'z'))
+#define isalnum(c) (isdigit(c) || isalpha(c))
+#define isprint(c) (c >= ' ' && c <= 126)
+
+struct lexeme current_lexeme;
+
+struct bucket {
+  int kind;
+  struct bucket * next;
+  int op;
+  char string[1];
+};
+
+#define HASHTABLE_SIZE 256
+
+static struct bucket * hashtable[HASHTABLE_SIZE] = { NULL, /*all nulls*/};
+
+#define BUFFER_SIZE 1024
+
+static char buffer[BUFFER_SIZE];
+
+#define STORE_BUFFER(pos,c) if (pos < sizeof(buffer) - 1) buffer[pos++] = c
+
+static inline unsigned int hash_ident(char * s)
+{
+  unsigned int h;
+  for (h = 0; *s != 0; s++) h = (h << 4) + h + *s;
+  return h % HASHTABLE_SIZE;
+}
+
+static void get_ident(int firstchar)
+{
+  int c, pos;
+  unsigned int h;
+  struct bucket * b;
+
+  /* Read characters of the ident */
+  buffer[0] = firstchar;
+  pos = 1;
+  while (1) {
+    c = getchar();
+    if (! (isalnum(c) || c == '-' || c == '_')) break;
+    STORE_BUFFER(pos, c);
+  }
+  if (c != EOF) ungetc(c, stdin);
+  buffer[pos] = 0;
+  /* Hash the ident */
+  h = hash_ident(buffer);
+  /* Look it up in the hash table */
+  for (b = hashtable[h]; b != NULL; b = b->next) {
+    if (strcmp(b->string, buffer) == 0) {
+      /* Found: return hash table entry */
+      current_lexeme.kind = b->kind;
+      if (b->kind == IDENTIFIER)
+        current_lexeme.u.s = b->string;
+      else
+        current_lexeme.u.op = b->op;
+      return;
+    }
+  }
+  /* Not found: enter it */
+  b = malloc(sizeof(struct bucket) + pos);
+  b->kind = IDENTIFIER;
+  strcpy(b->string, buffer);
+  b->next = hashtable[h];
+  hashtable[h] = b;
+  current_lexeme.kind = IDENTIFIER;
+  current_lexeme.u.s = b->string;
+}
+
+static void get_binder(void)
+{
+  int c = getchar();
+  if (! isalpha(c)) {
+    fprintf(stderr, "Bad binder /%c...\n", c);
+    exit(2);
+  }
+  get_ident(c);
+  if (current_lexeme.kind != IDENTIFIER) {
+    fprintf(stderr, "Binder /%s rebinds reserved identifier\n",
+            current_lexeme.u.s);
+    exit(2);
+  }
+  current_lexeme.kind = BINDER;
+}
+
+static int get_number(int firstchar)
+{
+  int c, pos, is_real;
+
+  pos = 0;
+  is_real = 0;
+  c = firstchar;
+  /* Initial sign */
+  if (c == '-') {
+    STORE_BUFFER(pos, c);
+    c = getchar();
+    if (! isdigit(c)) return -1;
+  }
+  /* Decimal number */
+  do {
+    STORE_BUFFER(pos, c);
+    c = getchar();
+  } while (isdigit(c));
+  /* Period and fractional part */
+  if (c == '.') {
+    is_real = 1;
+    STORE_BUFFER(pos, c);
+    c = getchar();
+    if (! isdigit(c)) return -1;
+    do {
+      STORE_BUFFER(pos, c);
+      c = getchar();
+    } while (isdigit(c));
+  }
+  /* Exponent */
+  if (c == 'e' || c == 'E') {
+    is_real = 1;
+    STORE_BUFFER(pos, c);
+    c = getchar();
+    if (c == '-') {
+      STORE_BUFFER(pos, c);
+      c = getchar();
+    }
+    if (! isdigit(c)) return -1;
+    do {
+      STORE_BUFFER(pos, c);
+      c = getchar();
+    } while (isdigit(c));
+  }
+  if (c != EOF) ungetc(c, stdin);
+  /* Convert to token */
+  buffer[pos] = 0;
+  if (is_real) {
+    current_lexeme.kind = REAL;
+    current_lexeme.u.d = strtod(buffer, NULL);
+  } else {
+    current_lexeme.kind = INTEGER;
+    current_lexeme.u.i = atoi(buffer);
+  }
+  return 0;
+}
+
+static int get_string()
+{
+  int c, pos;
+  pos = 0;
+  while (1) {
+    c = getchar();
+    if (c == '"') break;
+    if (! isprint(c)) return -1;
+    STORE_BUFFER(pos, c);
+  }
+  buffer[pos] = 0;
+  current_lexeme.kind = STRING;
+  current_lexeme.u.s = strdup(buffer);
+  return 0;
+}
+
+void get_lexeme(void)
+{
+  int c;
+
+  if (current_lexeme.kind != NONE) return;
+  while (1) {
+    c = getchar();
+    switch (c) {
+    case EOF:
+      current_lexeme.kind = END_OF_FILE; break;
+    case ' ': case '\n': case '\t': case '\r': case 11:
+      continue;
+    case '%':
+      do { c = getchar(); } while (c != '\n' && c != EOF);
+      continue;
+    case '/':
+      get_binder(); break;
+    case '-': case '0': case '1': case '2': case '3': case '4':
+    case '5': case '6': case '7': case '8': case '9':
+      if (get_number(c) == -1) {
+        fprintf(stderr, "Bad number\n");
+        exit(2);
+      }
+      break;
+    case '"':
+      if (get_string() == -1) {
+        fprintf(stderr, "Bad string literal\n");
+        exit(2);
+      }
+      break;
+    case '{':
+      current_lexeme.kind = LBRACE; break;
+    case '}':
+      current_lexeme.kind = RBRACE; break;
+    case '[':
+      current_lexeme.kind = LBRACKET; break;
+    case ']':
+      current_lexeme.kind = RBRACKET; break;
+    default:
+      if (isalpha(c)) {
+        get_ident(c);
+      } else {
+        fprintf(stderr, "Illegal character `%c'\n", c);
+        exit(2);
+      }
+      break;
+    }
+    return;
+  }
+}
+
+void discard_lexeme(void)
+{
+  current_lexeme.kind = NONE;
+}
+
+/* Enter the operators in the hashtable */
+
+struct op_init { char * name; int kind; int op; };
+
+static struct op_init operators[] =
+   { { "true", BOOLEAN, 1 },
+     { "false", BOOLEAN, 0 },
+     { "acos", OPERATOR, Op_acos },
+     { "addi", OPERATOR, Op_addi },
+     { "addf", OPERATOR, Op_addf },
+     { "apply", OPERATOR, Op_apply },
+     { "asin", OPERATOR, Op_asin },
+     { "clampf", OPERATOR, Op_clampf },
+     { "cone", OPERATOR, Op_cone },
+     { "cos", OPERATOR, Op_cos },
+     { "cube", OPERATOR, Op_cube },
+     { "cylinder", OPERATOR, Op_cylinder },
+     { "difference", OPERATOR, Op_difference },
+     { "divi", OPERATOR, Op_divi },
+     { "divf", OPERATOR, Op_divf },
+     { "eqi", OPERATOR, Op_eqi },
+     { "eqf", OPERATOR, Op_eqf },
+     { "floor", OPERATOR, Op_floor },
+     { "frac", OPERATOR, Op_frac },
+     { "get", OPERATOR, Op_get },
+     { "getx", OPERATOR, Op_getx },
+     { "gety", OPERATOR, Op_gety },
+     { "getz", OPERATOR, Op_getz },
+     { "if", OPERATOR, Op_if },
+     { "intersect", OPERATOR, Op_intersect },
+     { "length", OPERATOR, Op_length },
+     { "lessi", OPERATOR, Op_lessi },
+     { "lessf", OPERATOR, Op_lessf },
+     { "light", OPERATOR, Op_light },
+     { "modi", OPERATOR, Op_modi },
+     { "muli", OPERATOR, Op_muli },
+     { "mulf", OPERATOR, Op_mulf },
+     { "negi", OPERATOR, Op_negi },
+     { "negf", OPERATOR, Op_negf },
+     { "plane", OPERATOR, Op_plane },
+     { "point", OPERATOR, Op_point },
+     { "pointlight", OPERATOR, Op_pointlight },
+     { "real", OPERATOR, Op_real },
+     { "render", OPERATOR, Op_render },
+     { "rotatex", OPERATOR, Op_rotatex },
+     { "rotatey", OPERATOR, Op_rotatey },
+     { "rotatez", OPERATOR, Op_rotatez },
+     { "scale", OPERATOR, Op_scale },
+     { "sin", OPERATOR, Op_sin },
+     { "sphere", OPERATOR, Op_sphere },
+     { "spotlight", OPERATOR, Op_spotlight },
+     { "sqrt", OPERATOR, Op_sqrt },
+     { "subi", OPERATOR, Op_subi },
+     { "subf", OPERATOR, Op_subf },
+     { "translate", OPERATOR, Op_translate },
+     { "union", OPERATOR, Op_union },
+     { "uscale", OPERATOR, Op_uscale },
+     { "print", OPERATOR, Op_print },
+   };
+
+void init_lexer(void)
+{
+  int i;
+  for (i = 0; i < sizeof(operators) / sizeof(struct op_init); i++) {
+    struct op_init * opi = &(operators[i]);
+    struct bucket * b = malloc(sizeof(struct bucket) + strlen(opi->name));
+    unsigned int h = hash_ident(opi->name);
+    b->kind = opi->kind;
+    strcpy(b->string, opi->name);
+    b->op = opi->op;
+    b->next = hashtable[h];
+    hashtable[h] = b;
+  }
+}
+
diff --git a/test/raytracer/gmllexer.h b/test/raytracer/gmllexer.h
new file mode 100644
index 000000000..9b8d630f1
--- /dev/null
+++ b/test/raytracer/gmllexer.h
@@ -0,0 +1,21 @@
+/* Lexer for GML */
+
+struct lexeme {
+  enum {
+    NONE,
+    OPERATOR, IDENTIFIER, BINDER, BOOLEAN, INTEGER, REAL, STRING,
+    LBRACE, RBRACE, LBRACKET, RBRACKET, END_OF_FILE
+  } kind;
+  union {
+    int op;                     /* for operators */
+    char * s;                   /* for identifiers, binders, strings */
+    int i;                      /* for integer and boolean literals */
+    flt d;                      /* for float literals */
+  } u;
+};
+
+extern struct lexeme current_lexeme;
+
+extern void get_lexeme(void);
+extern void discard_lexeme(void);
+extern void init_lexer(void);
diff --git a/test/raytracer/gmlparser.c b/test/raytracer/gmlparser.c
new file mode 100644
index 000000000..b0a6cbc19
--- /dev/null
+++ b/test/raytracer/gmlparser.c
@@ -0,0 +1,102 @@
+/* Parser for GML */
+
+#include "config.h"
+
+#include "arrays.h"
+#include "gml.h"
+#include "gmllexer.h"
+#include "gmlparser.h"
+
+static struct array * parse_tokenlist(void);
+
+static int parse_token(struct tok * t)
+{
+  get_lexeme();
+  switch (current_lexeme.kind) {
+  case OPERATOR:
+    discard_lexeme();
+    t->tag = current_lexeme.u.op;
+    break;
+  case IDENTIFIER:
+    discard_lexeme();
+    t->tag = Identifier;
+    t->u.s = current_lexeme.u.s;
+    break;
+  case BINDER:
+    discard_lexeme();
+    t->tag = Binder;
+    t->u.s = current_lexeme.u.s;
+    break;
+  case BOOLEAN:
+    discard_lexeme();
+    t->tag = Boolean;
+    t->u.i = current_lexeme.u.i;
+    break;
+  case INTEGER:
+    discard_lexeme();
+    t->tag = Integer;
+    t->u.i = current_lexeme.u.i;
+    break;
+  case REAL:
+    discard_lexeme();
+    t->tag = Real;
+    t->u.d = current_lexeme.u.d;
+    break;
+  case STRING:
+    discard_lexeme();
+    t->tag = String;
+    t->u.s = current_lexeme.u.s;
+    break;
+  case LBRACE:
+    discard_lexeme();
+    t->tag = Function;
+    t->u.a = parse_tokenlist();
+    get_lexeme();
+    if (current_lexeme.kind != RBRACE) {
+      fprintf(stderr, "} expected\n");
+      exit(2);
+    }
+    discard_lexeme();
+    break;
+  case LBRACKET:
+    discard_lexeme();
+    t->tag = Array;
+    t->u.a = parse_tokenlist();
+    get_lexeme();
+    if (current_lexeme.kind != RBRACKET) {
+      fprintf(stderr, "] expected\n");
+      exit(2);
+    }
+    discard_lexeme();
+    break;
+  default:
+    return 0;
+  }
+  return 1;
+}
+
+static struct array * parse_tokenlist(void)
+{
+  struct array * a = alloc_array(sizeof(struct tok), 10);
+  struct tok t;
+  int i = 0;
+  while (parse_token(&t)) {
+    extend_array(struct tok, a);
+    set_array_large(struct tok, a, i, t);
+    i++;
+  }
+  return a;
+}
+
+struct array * parse_program(void)
+{
+  struct array * a = parse_tokenlist();
+  get_lexeme();
+  if (current_lexeme.kind != END_OF_FILE) {
+    fprintf(stderr, "syntax error at end of program\n");
+    exit(2);
+  }
+  return a;
+}
+
+  
diff --git a/test/raytracer/gmlparser.h b/test/raytracer/gmlparser.h
new file mode 100644
index 000000000..6a31546c9
--- /dev/null
+++ b/test/raytracer/gmlparser.h
@@ -0,0 +1,3 @@
+/* Parser for GML */
+
+struct array * parse_program(void);
diff --git a/test/raytracer/intersect.c b/test/raytracer/intersect.c
new file mode 100644
index 000000000..c8a212603
--- /dev/null
+++ b/test/raytracer/intersect.c
@@ -0,0 +1,416 @@
+#include "config.h"
+#include "arrays.h"
+#include "point.h"
+#include "vector.h"
+#include "eval.h"
+#include "object.h"
+#include "matrix.h"
+#include "intersect.h"
+
+/* Operations on interval lists */
+
+#define POS_INFTY 1e300
+
+struct intervlist {
+  flt beg;
+  struct object * obeg;
+  flt end;
+  struct object * oend;
+  struct intervlist * next;
+};
+
+typedef struct intervlist * intervlist;
+
+static intervlist cons(flt b, struct object * ob,
+                       flt e, struct object * oe,
+                       intervlist next)
+{
+  intervlist l = arena_alloc(sizeof(struct intervlist));
+  l->beg = b;
+  l->obeg = ob;
+  l->end = e;
+  l->oend = oe;
+  l->next = next;
+  return l;
+}
+
+static intervlist conshead(intervlist i, intervlist next)
+{
+  intervlist l = arena_alloc(sizeof(struct intervlist));
+  l->beg = i->beg;
+  l->obeg = i->obeg;
+  l->end = i->end;
+  l->oend = i->oend;
+  l->next = next;
+  return l;
+}
+
+static intervlist union_intervals(intervlist l1, intervlist l2)
+{
+  flt beg;
+  struct object * obeg;
+
+  if (l1 == NULL) return l2;
+  if (l2 == NULL) return l1;
+  if (l1->end < l2->beg)
+    /* int1 strictly before int2 */
+    return conshead(l1, union_intervals(l1->next, l2));
+  if (l2->end < l1->beg)
+    /* int2 strictly before int1 */
+    return conshead(l2, union_intervals(l1, l2->next));
+  /* int1 and int2 overlap, merge */
+  if (l1->beg < l2->beg) {
+    beg = l1->beg;
+    obeg = l1->obeg;
+  } else {
+    beg = l2->beg;
+    obeg = l2->obeg;
+  }
+  if (l1->end < l2->end)
+    return union_intervals(l1->next,
+                           cons(beg, obeg, l2->end, l2->oend, l2->next));
+  else
+    return union_intervals(cons(beg, obeg, l1->end, l1->oend, l1->next),
+                           l2->next);
+}
+
+static intervlist intersect_intervals(intervlist l1, intervlist l2)
+{
+  flt beg;
+  struct object * obeg;
+
+  if (l1 == NULL) return NULL;
+  if (l2 == NULL) return NULL;
+  if (l1->end < l2->beg)
+    /* int1 strictly before int2 */
+    return intersect_intervals(l1->next, l2);
+  if (l2->end < l1->beg)
+    /* int2 strictly before int1 */
+    return intersect_intervals(l1, l2->next);
+  /* int1 and int2 overlap, add intersection */
+  if (l1->beg > l2->beg) {
+    beg = l1->beg;
+    obeg = l1->obeg;
+  } else {
+    beg = l2->beg;
+    obeg = l2->obeg;
+  }
+  if (l1->end < l2->end)
+    return cons(beg, obeg, l1->end, l1->oend, intersect_intervals(l1->next, l2));
+  else
+    return cons(beg, obeg, l2->end, l2->oend, intersect_intervals(l1, l2->next));
+}
+
+static intervlist difference_intervals(intervlist l1, intervlist l2)
+{
+  intervlist l;
+
+  if (l1 == NULL) return NULL;
+  if (l2 == NULL) return l1;
+  if (l1->end < l2->beg)
+    /* int1 strictly before int2, keep int1 */
+    return conshead(l1, difference_intervals(l1->next, l2));
+  if (l2->end < l1->beg)
+    /* int2 strictly before int1, throw int2 away */
+    return difference_intervals(l1, l2->next);
+  /* int and int2 overlap */
+  if (l1->end > l2->end)
+    /* int1 extends beyond int2, add segment [l2->end,l1->end] */
+    l = difference_intervals(cons(l2->end, l2->oend, l1->end, l1->oend, l1->next), l2->next);
+  else
+    /* int2 doesn't extend beyond int2, nothing to add */
+    l = difference_intervals(l1->next, l2);
+  if (l1->beg < l2->beg)
+    /* int1 starts before l2, add segment [l1->beg,l2->beg] */
+    l = cons(l1->beg, l1->obeg, l2->beg, l2->obeg, l);
+  return l;
+}
+
+/* Intersections between rays (p + tv) and basic objects (bobj) */
+
+static intervlist intersect_ray_plane(struct object * bobj,
+                                      struct point * p,
+                                      struct vector * v)
+{
+  if (p->y >= 0)
+    if (v->dy >= 0)
+      return NULL;
+    else
+      return cons(- p->y / v->dy, bobj, POS_INFTY, bobj, NULL);
+  else
+    if (v->dy >= 0)
+      return cons(0, bobj, - p->y / v->dy, bobj, NULL);
+    else
+      return cons(0, bobj, POS_INFTY, bobj, NULL);
+}
+
+static intervlist intersect_ray_sphere(struct object * bobj,
+                                       struct point * p,
+                                       struct vector * v)
+{
+  flt a, b, c, d, sq, t0, t1, tswap;
+
+  a = v->dx * v->dx + v->dy * v->dy + v->dz * v->dz;
+  b = p->x * v->dx + p->y * v->dy + p->z * v->dz;
+  c = p->x * p->x + p->y * p->y + p->z * p->z - 1.0;
+  d = b * b - a * c;
+  if (d <= 0) return NULL;
+  sq = sqrt(d);
+  /* Numerically stable solution of quadratic */
+  if (b >= 0) {
+    t0 = c / (- b - sq);
+    t1 = (- b - sq) / a;
+  } else {
+    t0 = c / (- b + sq);
+    t1 = (- b + sq) / a;
+  }
+  if (t0 >= t1) {
+    tswap = t0; t0 = t1; t1 = tswap;
+  }
+  if (t1 <= 0) {
+    assert (t0 <= 0);
+    return NULL;
+  }
+  if (t0 >= 0) {
+    assert (t1 >= 0);
+    return cons(t0, bobj, t1, bobj, NULL);
+  }
+  return cons(0, bobj, t1, bobj, NULL);
+}
+
+static intervlist intersect_ray_slice(struct object * bobj,
+                                      flt pc, flt vc)
+{
+  if (pc > 1.0) {
+    if (vc >= 0.0)
+      return NULL;
+    else
+      return cons((1.0 - pc) / vc, bobj, - pc / vc, bobj, NULL);
+  }
+  if (pc < 0.0) {
+    if (vc <= 0.0)
+      return NULL;
+    else
+      return cons(- pc / vc, bobj, (1.0 - pc) / vc, bobj, NULL);
+  }
+  if (vc == 0.0)
+    return cons(0.0, bobj, POS_INFTY, bobj, NULL);
+  if (vc > 0.0)
+    return cons(0.0, bobj, (1.0 - pc) / vc, bobj, NULL);
+  else
+    return cons(0.0, bobj, - pc / vc, bobj, NULL);
+}
+
+static intervlist intersect_ray_cube(struct object * bobj,
+                                     struct point * p,
+                                     struct vector * v)
+{
+  return intersect_intervals(intersect_ray_slice(bobj, p->x, v->dx),
+           intersect_intervals(intersect_ray_slice(bobj, p->y, v->dy),
+                               intersect_ray_slice(bobj, p->z, v->dz)));
+}
+
+static intervlist intersect_ray_infinite_cylinder(struct object * bobj,
+                                                  struct point * p,
+                                                  struct vector * v)
+{
+  flt a, b, c, d, sq, t0, t1, tswap;
+
+  a = v->dx * v->dx + v->dz * v->dz;
+  b = p->x * v->dx + p->z * v->dz;
+  c = p->x * p->x + p->z * p->z - 1.0;
+  d = b * b - a * c;
+  if (d <= 0.0) return NULL;
+  sq = sqrt(d);
+  if (b >= 0.0) {
+    t0 = c / (- b - sq);
+    t1 = (- b - sq) / a;
+  } else {
+    t0 = c / (- b + sq);
+    t1 = (- b + sq) / a;
+  }
+  if (t0 >= t1) { tswap = t0; t0 = t1; t1 = tswap; }
+  if (t1 <= 0.0) {
+    assert (t0 <= 0.0);
+    return NULL;
+  }
+  if (t0 >= 0.0) {
+    assert (t1 >= 0.0);
+    return cons(t0, bobj, t1, bobj, NULL);
+  }
+  return cons(0.0, bobj, t1, bobj, NULL);
+}
+
+static intervlist intersect_ray_cylinder(struct object * bobj,
+                                         struct point * p,
+                                         struct vector * v)
+{
+  return intersect_intervals(intersect_ray_infinite_cylinder(bobj, p, v),
+                             intersect_ray_slice(bobj, p->y, v->dy));
+}
+
+static intervlist intersect_ray_infinite_cone(struct object * bobj,
+                                              struct point * p,
+                                              struct vector * v)
+{
+  flt a, b, c, d, sq, t, t0, t1, tswap;
+
+  a = v->dx * v->dx - v->dy * v->dy + v->dz * v->dz;
+  b = p->x * v->dx - p->y * v->dy + p->z * v->dz;
+  c = p->x * p->x - p->y * p->y + p->z * p->z;
+  if (a == 0.0) {
+    if (b == 0.0) return NULL;
+    t = - 0.5 * c / b;
+    if (c < 0.0) {
+      if (t <= 0.0)
+        return cons(0.0, bobj, POS_INFTY, bobj, NULL);
+      else
+        return cons(0.0, bobj, t, bobj, NULL);
+    }
+    if (t <= 0.0)
+      return NULL;
+    else
+      return cons(t, bobj, POS_INFTY, bobj, NULL);
+  }
+  d = b * b - a * c;
+  if (d <= 0.0) return NULL;
+  sq = sqrt(d);
+  if (b >= 0.0) {
+    t0 = c / (- b - sq);
+    t1 = (- b - sq) / a;
+  } else {
+    t0 = c / (- b + sq);
+    t1 = (- b + sq) / a;
+  }
+  if (t0 >= t1) { tswap = t0; t0 = t1; t1 = tswap; }
+  if (t1 <= 0.0) {
+      assert (t0 <= 0.0);
+      if (c < 0.0)
+        return cons(0.0, bobj, POS_INFTY, bobj, NULL);
+      else
+        return NULL;
+  }
+  if (t0 >= 0.0) {
+    assert (t1 >= 0.0);
+    if (c < 0.0)
+      return cons(0.0, bobj, t0, bobj,
+                  cons (t1, bobj, POS_INFTY, bobj, NULL));
+    else
+      return cons(t0, bobj, t1, bobj, NULL);
+  }
+  if (c < 0.0)
+    return cons(0.0, bobj, t1, bobj, NULL);
+  else
+    return cons(t1, bobj, POS_INFTY, bobj, NULL);
+}
+
+static intervlist intersect_ray_cone(struct object * bobj,
+                                     struct point * p,
+                                     struct vector * v)
+{
+  return intersect_intervals(intersect_ray_infinite_cone(bobj, p, v),
+                             intersect_ray_slice(bobj, p->y, v->dy));
+}
+
+/* Approximate test based on bounding sphere */
+
+static inline int inside_bs(struct point * p,
+                            struct vector * v,
+                            struct object * obj)
+{
+  flt x, y, z, a, b, c;
+
+  assert(obj->radius >= 0.0);
+
+  if (obj->radius >= POS_INFTY) return 1;
+  /* Shift origin to obj.center */
+  x = p->x - obj->center.x;
+  y = p->y - obj->center.y;
+  z = p->z - obj->center.z;
+  /* Check whether quadratic has positive discriminant */
+  a = v->dx * v->dx + v->dy * v->dy + v->dz * v->dz;
+  b = x * v->dx + y * v->dy + z * v->dz;
+  c = x * x + y * y + z * z - obj->radius * obj->radius;
+  return (b * b > a * c);
+}
+
+/* Interval list representing the intersection of the given ray
+   and the given composite object */
+
+static intervlist intersect_ray_object(struct point * p,
+                                       struct vector * v,
+                                       struct object * obj)
+{
+#define CONVERT_V_P \
+    apply_to_point(obj->world2obj, p, &p2); \
+    apply_to_vect(obj->world2obj, v, &v2)
+  struct point p2;
+  struct vector v2;
+  intervlist res;
+
+  /* Fast, approximate test based on bounding sphere */
+  if (! inside_bs(p, v, obj)) return NULL;
+  /* Slow, exact test */
+  switch (obj->kind) {
+  case Cone:
+    CONVERT_V_P;
+    res = intersect_ray_cone(obj, &p2, &v2);
+    break;
+  case Cube:
+    CONVERT_V_P;
+    res = intersect_ray_cube(obj, &p2, &v2);
+    break;
+  case Cylinder:
+    CONVERT_V_P;
+    res = intersect_ray_cylinder(obj, &p2, &v2);
+    break;
+  case Plane:
+    CONVERT_V_P;
+    res = intersect_ray_plane(obj, &p2, &v2);
+    break;
+  case Sphere:
+    CONVERT_V_P;
+    res = intersect_ray_sphere(obj, &p2, &v2);
+    break;
+  case Union:
+    res = union_intervals(intersect_ray_object(p, v, obj->o1),
+                           intersect_ray_object(p, v, obj->o2));
+    break;
+  case Intersection:
+    res = intersect_intervals(intersect_ray_object(p, v, obj->o1),
+                               intersect_ray_object(p, v, obj->o2));
+    break;
+  case Difference:
+    res = difference_intervals(intersect_ray_object(p, v, obj->o1),
+                                intersect_ray_object(p, v, obj->o2));
+    break;
+  default:
+    assert(0);
+  }
+#undef CONVERT_V_P
+  return res;
+}
+
+/* Return the closest base object intersecting the given ray, and
+   the curvilinear abscissa t of the intersection point.
+   Return NULL if no intersection.
+   Return t = 0.0 if the viewpoint (origin of the ray) is inside an
+   object. */
+
+struct object * intersect_ray(struct point * p,
+                              struct vector * v,
+                              struct object * obj,
+                              int initial,
+                              /*out*/ flt * t)
+{
+  intervlist l = intersect_ray_object(p, v, obj);
+  if (l == NULL) return NULL;
+  assert(l->beg >= 0.0);
+  if (l->beg <= 0.0 && !initial) {
+    /* skip first intersection */
+    l = l->next;
+    if (l == NULL) return NULL;
+  }
+  *t = l->beg;
+  return l->obeg;
+}
diff --git a/test/raytracer/intersect.h b/test/raytracer/intersect.h
new file mode 100644
index 000000000..35241a604
--- /dev/null
+++ b/test/raytracer/intersect.h
@@ -0,0 +1,11 @@
+/* Return the closest base object intersecting the given ray, and
+   the curvilinear abscissa t of the intersection point.
+   Return NULL if no intersection.
+   Return t = 0.0 if the viewpoint (origin of the ray) is inside an
+   object. */
+
+struct object * intersect_ray(struct point * p,
+                              struct vector * v,
+                              struct object * obj,
+                              int initial,
+                              /*out*/ flt * t);
diff --git a/test/raytracer/kal.gml b/test/raytracer/kal.gml
new file mode 100644
index 000000000..e9d581525
--- /dev/null
+++ b/test/raytracer/kal.gml
@@ -0,0 +1,78 @@
+% kal.gml
+
+1.0 0.7 0.7 point /red
+0.7 1.0 0.7 point /green
+0.7 0.7 1.0 point /blue
+
+{ /color
+  { /v /u /face
+    color 0.1 0.99 1.0
+  }
+} /mirror
+
+                                % a cube
+{ /v /u /face                     % bind arguments
+  0.9 1.0 0.9 point               % surface color
+  0.1 0.9 1.0                     % kd ks n
+} cube
+
+0.0 -0.5 35.0 translate
+                                % a sphere
+{ /v /u /face                     % bind arguments
+  0.9 0.9 1.0 point               % surface color
+  0.1 0.9 1.0                     % kd ks n
+} sphere
+
+-1.0 0.0 20.0 translate
+
+union
+                                % a sphere
+{ /v /u /face                     % bind arguments
+  1.0 1.0 0.9 point               % surface color
+  0.1 0.9 1.0                     % kd ks n
+} cylinder
+20.0 rotatex
+0.5 0.5 15.0 translate
+union 
+
+blue mirror apply plane
+0.0 -2.0 0.0 translate
+
+union
+
+green mirror apply plane
+0.0 -2.0 0.0 translate
+120.0 rotatez
+
+union
+
+red mirror apply plane
+0.0 -2.0 0.0 translate
+240.0 rotatez
+
+union
+
+{ /v /u /face
+  0.4 0.4 0.4 point
+  0.0 0.0 1.0
+} plane
+
+90.0 rotatex
+0.0 0.0 3.0 translate
+difference
+
+0.0 -1.0 0.0 translate
+
+/scene
+
+0.0 0.0 -2.0 point 0.9 0.9 0.9 point pointlight /light1
+
+0.4 0.4	0.4 point	% ambient
+[light1]		% lights
+scene			% object
+1000			% depth
+90.0			% fov
+400 400 		% wid ht
+"kal.ppm"		% output file
+render
+
diff --git a/test/raytracer/light.c b/test/raytracer/light.c
new file mode 100644
index 000000000..1c90104ce
--- /dev/null
+++ b/test/raytracer/light.c
@@ -0,0 +1,126 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+#include "eval.h"
+#include "object.h"
+#include "intersect.h"
+#include "light.h"
+
+struct light * dirlight(struct point * dir, struct point * c)
+{
+  struct light * l = arena_alloc(sizeof(struct light));
+  l->kind = Directional;
+  ASSIGN(l->u.directional.dir, *dir);
+  ASSIGN(l->u.directional.col, *c);
+  return l;
+}
+
+struct light * pointlight(struct point * orig, struct point * c)
+{
+  struct light * l = arena_alloc(sizeof(struct light));
+  l->kind = Pointlight;
+  ASSIGN(l->u.point.orig, *orig);
+  ASSIGN(l->u.point.col, *c);
+  return l;
+}
+
+struct light * spotlight(struct point * pos, struct point * at,
+                         struct point * c, flt cutoff, flt exp)
+{
+  struct light * l = arena_alloc(sizeof(struct light));
+  struct vector uv;
+  l->kind = Spot;
+  ASSIGN(l->u.spot.orig, *pos);
+  ASSIGN(l->u.spot.at, *at);
+  ASSIGN(l->u.spot.col, *c);
+  l->u.spot.cutoff = cutoff;
+  l->u.spot.exponent = exp;
+  between(at, pos, &uv);
+  vnormalize(&uv, &(l->u.spot.unit));
+  return l;
+}
+
+static void color_from_light(struct object * scene,
+                             struct point * pt,
+                             struct vector * v,
+                             struct vector * n,
+                             flt kd, flt ks, flt phong,
+                             struct light * light,
+                             /*out*/ struct point * il)
+{
+  struct vector dir;
+  struct object * obj;
+  struct vector L, H, vnorm;
+  struct point i;
+  flt t, att, dp, m;
+
+  /* Compute direction towards light source */
+  switch (light->kind) {
+  case Directional:
+    dir.dx = - light->u.directional.dir.x;
+    dir.dy = - light->u.directional.dir.y;
+    dir.dz = - light->u.directional.dir.z;
+    break;
+  case Pointlight:
+    between(pt, &light->u.point.orig, &dir);
+    break;
+  case Spot:
+    between(pt, &light->u.spot.orig, &dir);
+    break;
+  }
+  /* Check that light source is not behind us */
+  if (dotproduct(&dir, n) <= 0.0) return; /*no contribution*/
+  /* Check that no object blocks the light */
+  obj = intersect_ray(pt, &dir, scene, 0 /*??*/, &t);
+  if (obj != NULL && (light->kind == Directional || t <= 1.0)) return;
+  /* Compute the L and H vectors */
+  vnormalize(&dir, &L);
+  vnormalize(v, &vnorm);
+  vsub(&L, &vnorm, &H);
+  vnormalize(&H, &H);
+  /* Intensity of light source at object */
+  switch (light->kind) {
+  case Directional:
+    i.x = light->u.directional.col.x;
+    i.y = light->u.directional.col.y;
+    i.z = light->u.directional.col.z;
+    break;
+  case Pointlight:
+    att = 100.0 / (99.0 + dist2(pt, &light->u.point.orig));
+    i.x = light->u.point.col.x * att;
+    i.y = light->u.point.col.y * att;
+    i.z = light->u.point.col.z * att;
+    break;
+  case Spot:
+    dp = dotproduct(&L, &light->u.spot.unit);
+    if (acos(dp) >= light->u.spot.cutoff) return;
+    att =
+      pow(dp, light->u.spot.exponent)
+      * (100.0 / (99.0 + dist2(pt, &light->u.spot.orig)));
+    i.x = light->u.spot.col.x * att;
+    i.y = light->u.spot.col.y * att;
+    i.z = light->u.spot.col.z * att;
+    break;
+  }
+  /* Add final contribution */
+  m = kd * dotproduct(n, &L) + ks * pow(dotproduct(n, &H), phong);
+  il->x += m * i.x;
+  il->y += m * i.y;
+  il->z += m * i.z;
+}
+
+void color_from_lights(struct object * scene,
+                       struct point * p,
+                       struct vector * v,
+                       struct vector * n,
+                       flt kd, flt ks, flt phong,
+                       struct light ** lights,
+                       /*out*/ struct point * il)
+{
+  int i;
+
+  il->x = il->y = il->z = 0.0;
+  for (i = 0; lights[i] != NULL; i++)
+    color_from_light(scene, p, v, n, kd, ks, phong, lights[i], il);
+}
+
diff --git a/test/raytracer/light.h b/test/raytracer/light.h
new file mode 100644
index 000000000..9a1d3e2aa
--- /dev/null
+++ b/test/raytracer/light.h
@@ -0,0 +1,33 @@
+struct light {
+  enum { Directional, Pointlight, Spot } kind;
+  union {
+    struct {
+      struct point dir;
+      struct point col;
+    } directional;
+    struct {
+      struct point orig;
+      struct point col;
+    } point;
+    struct {
+      struct point orig;
+      struct point at;
+      struct point col;
+      flt cutoff, exponent;
+      struct vector unit;
+    } spot;
+  } u;
+};
+
+struct light * dirlight(struct point * dir, struct point * c);
+struct light * pointlight(struct point * orig, struct point * c);
+struct light * spotlight(struct point * pos, struct point * at,
+                         struct point * c, flt a, flt cutoff);
+
+void color_from_lights(struct object * obj,
+                       struct point * p,
+                       struct vector * v,
+                       struct vector * n,
+                       flt kd, flt ks, flt phong,
+                       struct light ** lights,
+                       /*out*/ struct point * il);
diff --git a/test/raytracer/main.c b/test/raytracer/main.c
new file mode 100644
index 000000000..1efd396de
--- /dev/null
+++ b/test/raytracer/main.c
@@ -0,0 +1,13 @@
+#include "config.h"
+#include "arrays.h"
+#include "gmllexer.h"
+#include "gmlparser.h"
+#include "eval.h"
+
+int main(int argc, char ** argv)
+{
+  arena_init();
+  init_lexer();
+  execute_program(parse_program());
+  return 0;
+}
diff --git a/test/raytracer/matrix.c b/test/raytracer/matrix.c
new file mode 100644
index 000000000..881434689
--- /dev/null
+++ b/test/raytracer/matrix.c
@@ -0,0 +1,102 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+#include "matrix.h"
+
+struct matrix matrix_identity = {
+  1.0,  0.0,  0.0,  0.0,
+  0.0,  1.0,  0.0,  0.0,
+  0.0,  0.0,  1.0,  0.0,
+};
+
+void apply_to_point(struct matrix * m, struct point * p,
+                                  /*out*/ struct point * r)
+{
+  r->x = m->xx * p->x + m->xy * p->y + m->xz * p->z + m->xt;
+  r->y = m->yx * p->x + m->yy * p->y + m->yz * p->z + m->yt;
+  r->z = m->zx * p->x + m->zy * p->y + m->zz * p->z + m->zt;
+}
+
+void apply_to_vect(struct matrix * m, struct vector * v,
+                                  /*out*/ struct vector * r)
+{
+  r->dx = m->xx * v->dx + m->xy * v->dy + m->xz * v->dz;
+  r->dy = m->yx * v->dx + m->yy * v->dy + m->yz * v->dz;
+  r->dz = m->zx * v->dx + m->zy * v->dy + m->zz * v->dz;
+}
+
+
+
+struct matrix * mtranslate(flt sx, flt sy, flt sz)
+{
+  struct matrix * m = arena_alloc(sizeof(struct matrix));
+  m->xx = 1.0;  m->xy = 0.0;  m->xz = 0.0;  m->xt = sx;
+  m->yx = 0.0;  m->yy = 1.0;  m->yz = 0.0;  m->yt = sy;
+  m->zx = 0.0;  m->zy = 0.0;  m->zz = 1.0;  m->zt = sz;
+  return m;
+}
+
+struct matrix * mscale(flt sx, flt sy, flt sz)
+{
+  struct matrix * m = arena_alloc(sizeof(struct matrix));
+  m->xx = sx;   m->xy = 0.0;  m->xz = 0.0;  m->xt = 0.0;
+  m->yx = 0.0;  m->yy = sy ;  m->yz = 0.0;  m->yt = 0.0;
+  m->zx = 0.0;  m->zy = 0.0;  m->zz = sz ;  m->zt = 0.0;
+  return m;
+}
+
+struct matrix * mrotatex(flt a)
+{
+  struct matrix * m = arena_alloc(sizeof(struct matrix));
+  flt c = cos(a);
+  flt s = sin(a);
+  m->xx = 1.0;  m->xy = 0.0;  m->xz = 0.0;  m->xt = 0.0;
+  m->yx = 0.0;  m->yy = c;    m->yz = - s;  m->yt = 0.0;
+  m->zx = 0.0;  m->zy = s;    m->zz = c;    m->zt = 0.0;
+  return m;
+}
+
+struct matrix * mrotatey(flt a)
+{
+  struct matrix * m = arena_alloc(sizeof(struct matrix));
+  flt c = cos(a);
+  flt s = sin(a);
+  m->xx = c;    m->xy = 0.0;  m->xz = s;    m->xt = 0.0;
+  m->yx = 0.0;  m->yy = 1.0;  m->yz = 0.0;  m->yt = 0.0;
+  m->zx = - s;  m->zy = 0.0;  m->zz = c;    m->zt = 0.0;
+  return m;
+}
+
+struct matrix * mrotatez(flt a)
+{
+  struct matrix * m = arena_alloc(sizeof(struct matrix));
+  flt c = cos(a);
+  flt s = sin(a);
+  m->xx = c;    m->xy = - s;  m->xz = 0.0;  m->xt = 0.0;
+  m->yx = s;    m->yy = c;    m->yz = 0.0;  m->yt = 0.0;
+  m->zx = 0.0;  m->zy = 0.0;  m->zz = 1.0;  m->zt = 0.0;
+  return m;
+}
+
+struct matrix * mcompose(struct matrix * m, struct matrix * n)
+{
+  struct matrix * r = arena_alloc(sizeof(struct matrix));
+
+  r->xx = m->xx * n->xx + m->xy * n->yx + m->xz * n->zx;
+  r->xy = m->xx * n->xy + m->xy * n->yy + m->xz * n->zy;
+  r->xz = m->xx * n->xz + m->xy * n->yz + m->xz * n->zz;
+  r->xt = m->xx * n->xt + m->xy * n->yt + m->xz * n->zt + m->xt;
+
+  r->yx = m->yx * n->xx + m->yy * n->yx + m->yz * n->zx;
+  r->yy = m->yx * n->xy + m->yy * n->yy + m->yz * n->zy;
+  r->yz = m->yx * n->xz + m->yy * n->yz + m->yz * n->zz;
+  r->yt = m->yx * n->xt + m->yy * n->yt + m->yz * n->zt + m->yt;
+
+  r->zx = m->zx * n->xx + m->zy * n->yx + m->zz * n->zx;
+  r->zy = m->zx * n->xy + m->zy * n->yy + m->zz * n->zy;
+  r->zz = m->zx * n->xz + m->zy * n->yz + m->zz * n->zz;
+  r->zt = m->zx * n->xt + m->zy * n->yt + m->zz * n->zt + m->zt;
+
+  return r;
+}
+
diff --git a/test/raytracer/matrix.h b/test/raytracer/matrix.h
new file mode 100644
index 000000000..4b9f4346e
--- /dev/null
+++ b/test/raytracer/matrix.h
@@ -0,0 +1,18 @@
+struct matrix {
+  flt xx, xy, xz, xt;
+  flt yx, yy, yz, yt;
+  flt zx, zy, zz, zt;
+};
+
+void apply_to_point(struct matrix * m, struct point * p,
+                    /*out*/ struct point * r);
+void apply_to_vect(struct matrix * m, struct vector * v,
+                   /*out*/ struct vector * r);
+extern struct matrix matrix_identity;
+#define mid (&matrix_identity)
+struct matrix * mtranslate(flt sx, flt sy, flt sz);
+struct matrix * mscale(flt sx, flt sy, flt sz);
+struct matrix * mrotatex(flt a);
+struct matrix * mrotatey(flt a);
+struct matrix * mrotatez(flt a);
+struct matrix * mcompose(struct matrix * m, struct matrix * n);
diff --git a/test/raytracer/memory.c b/test/raytracer/memory.c
new file mode 100644
index 000000000..fe194cc91
--- /dev/null
+++ b/test/raytracer/memory.c
@@ -0,0 +1,60 @@
+/* Arena-based memory management */
+
+#include "config.h"
+
+#define SIZE_AREA 1024000
+
+struct area {
+  struct area * next;
+  char data[SIZE_AREA];
+};
+
+struct area * arena_head;
+struct area * arena_cur;
+int arena_cur_ofs;
+
+void arena_init(void)
+{
+  arena_head = malloc(sizeof(struct area));
+  if (arena_head == NULL) {
+    fprintf(stderr, "Out of memory\n");
+    exit(2);
+  }
+  arena_head->next = NULL;
+  arena_cur = arena_head;
+  arena_cur_ofs = 0;
+}
+
+void arena_clear(void)
+{
+  arena_cur = arena_head;
+  arena_cur_ofs = 0;
+}
+
+void * arena_alloc(int size)
+{
+  char * res;
+  struct area * n;
+
+  if (size >= SIZE_AREA) {
+    fprintf(stderr, "Memory request too large (%d)\n", size);
+    exit(2);
+  }
+  if (arena_cur_ofs + size <= SIZE_AREA) {
+    res = arena_cur->data + arena_cur_ofs;
+    arena_cur_ofs += size;
+    return res;
+  }
+  if (arena_cur->next == NULL) {
+    n = malloc(sizeof(struct area));
+    if (n == NULL) {
+      fprintf(stderr, "Out of memory\n");
+      exit(2);
+    }
+    n->next = NULL;
+    arena_cur->next = n;
+  }
+  arena_cur = arena_cur->next;
+  arena_cur_ofs = size;
+  return arena_cur->data;
+}
diff --git a/test/raytracer/object.c b/test/raytracer/object.c
new file mode 100644
index 000000000..3dd6e77e4
--- /dev/null
+++ b/test/raytracer/object.c
@@ -0,0 +1,214 @@
+#include "config.h"
+#include "arrays.h"
+#include "point.h"
+#include "vector.h"
+#include "eval.h"
+#include "object.h"
+#include "matrix.h"
+
+static struct object * new_object(int kind, struct closure * c)
+{
+  struct object * o = arena_alloc(sizeof(struct object));
+  o->kind = kind;
+  ASSIGN(o->surf, *c);
+  o->world2obj = o->obj2world = mid;
+  o->max_scale_applied = 1.0;
+  o->radius = BS_NOT_COMPUTED;
+  return o;
+}
+
+struct object * cone(struct closure * c)
+{ return new_object(Cone, c); }
+
+struct object * cube(struct closure * c)
+{ return new_object(Cube, c); }
+
+struct object * cylinder(struct closure * c)
+{ return new_object(Cylinder, c); }
+
+struct object * plane(struct closure * c)
+{ return new_object(Plane, c); }
+
+struct object * sphere(struct closure * c)
+{ return new_object(Sphere, c); }
+
+static struct object * transform(struct object * o,
+                                 struct matrix * t,
+                                 struct matrix * tinv,
+                                 flt scale)
+{
+  struct object * no = arena_alloc(sizeof(struct object));
+  no->kind = o->kind;
+  switch (o->kind) {
+  case Union:
+  case Intersection:
+  case Difference:
+    no->o1 = transform(o->o1, t, tinv, scale);
+    no->o2 = transform(o->o2, t, tinv, scale);
+    break;
+  default:
+    ASSIGN(no->surf, o->surf);
+    no->world2obj = mcompose(o->world2obj, tinv);
+    no->obj2world = mcompose(t, o->obj2world);
+    no->max_scale_applied = o->max_scale_applied * scale;
+  }
+  no->radius = BS_NOT_COMPUTED;
+  return no;
+}
+
+struct object * orotatex(struct object * o1, flt a)
+{
+  return transform(o1, mrotatex(a), mrotatex(-a), 1.0);
+}
+
+struct object * orotatey(struct object * o1, flt a)
+{
+  return transform(o1, mrotatey(a), mrotatey(-a), 1.0);
+}
+
+struct object * orotatez(struct object * o1, flt a)
+{
+  return transform(o1, mrotatez(a), mrotatez(-a), 1.0);
+}
+
+struct object * oscale(struct object * o1, flt sx, flt sy, flt sz)
+{
+  flt scale = sx;
+  if (sy > scale) scale = sy;
+  if (sz > scale) scale = sz;
+  return transform(o1, mscale(sx, sy, sz), mscale(1 / sx, 1 / sy, 1 / sz),
+                   scale);
+}
+
+struct object * otranslate(struct object * o1,
+                           flt tx, flt ty, flt tz)
+{
+  return transform(o1, mtranslate(tx, ty, tz), mtranslate(- tx, - ty, - tz), 
+                   1.0);
+}
+
+struct object * ouscale(struct object * o1, flt s)
+{
+  flt sinv = 1 / s;
+  return transform(o1, mscale(s, s, s), mscale(sinv, sinv, sinv), s);
+}
+
+struct object * odifference(struct object * o1, struct object * o2)
+{
+  struct object * o = arena_alloc(sizeof(struct object));
+  o->kind = Difference;
+  o->o1 = o1;
+  o->o2 = o2;
+  o->radius = BS_NOT_COMPUTED;
+  return o;
+}
+
+struct object * ointersect(struct object * o1, struct object * o2)
+{
+  struct object * o = arena_alloc(sizeof(struct object));
+  o->kind = Intersection;
+  o->o1 = o1;
+  o->o2 = o2;
+  o->radius = BS_NOT_COMPUTED;
+  return o;
+}
+
+struct object * ounion(struct object * o1, struct object * o2)
+{
+  struct object * o = arena_alloc(sizeof(struct object));
+  o->kind = Union;
+  o->o1 = o1;
+  o->o2 = o2;
+  o->radius = BS_NOT_COMPUTED;
+  return o;
+}
+
+static void normal_vector_object(struct object * obj, 
+                                 struct point * p,
+                                 int face,
+                                 /*out*/ struct vector * v)
+{
+  switch (obj->kind) {
+  case Cone:
+    if (face == 0) {
+      v->dx = p->x; v->dy = - p->y; v->dz = p->z;
+    } else {
+      v->dx = 0; v->dy = 1; v->dz = 0;
+    }
+    break;
+  case Cube:
+    switch (face) {
+    case 0:
+      v->dx = 0; v->dy = 0; v->dz = -1; break;
+    case 1:
+      v->dx = 0; v->dy = 0; v->dz = 1; break;
+    case 2:
+      v->dx = -1; v->dy = 0; v->dz = 0; break;
+    case 3:
+      v->dx = 1; v->dy = 0; v->dz = 0; break;
+    case 4:
+      v->dx = 0; v->dy = 1; v->dz = 0; break;
+    case 5:
+      v->dx = 0; v->dy = -1; v->dz = 0; break;
+    }
+    break;
+  case Cylinder:
+    switch (face) {
+    case 0:
+      v->dx = p->x; v->dy = 0; v->dz = p->z; break;
+    case 1:
+      v->dx = 0; v->dy = 1; v->dz = 0; break;
+    case 2:
+      v->dx = 0; v->dy = -1; v->dz = 0; break;
+    }
+    break;
+  case Plane:
+    v->dx = 0; v->dy = 1; v->dz = 0; break;
+  case Sphere:
+    v->dx = p->x; v->dy = p->y; v->dz = p->z; break;
+  default:
+    assert(0);
+  }
+}
+
+static void tangent_vectors(struct vector * n,
+                            /*out*/ struct vector * t1,
+                            /*out*/ struct vector * t2)
+{
+  if (n->dy > 0) {
+    t1->dx = n->dy; t1->dy = - n->dx; t1->dz = 0;
+    t2->dx = 0; t2->dy = n->dz; t2->dz = - n->dy;
+  }
+  /*   y   0      xy      x
+      -x ^ z   =  yy  = y y
+       0  -y      yz      z */
+  else if (n->dy == 0) {
+    t1->dx = n->dz; t1->dy = 0; t1->dz = - n->dx;
+    t2->dx = n->dz; t2->dy = 1; t2->dz = - n->dx;
+  }
+  /*   z   z       x      x
+       0 ^ 1   =   0  = 1 y
+      -x  -x       z      z */
+  else {
+    t1->dx = n->dy; t1->dy = - n->dx; t1->dz = 0;
+    t2->dx = 0; t2->dy = - n->dz; t2->dz = n->dy;
+  }
+  /*   y    0      -xy        x
+      -x ^ -z  =   -yy = (-y) y
+       0    y      -zy        z */
+}
+
+void normal_vector(struct object * obj, struct point * p, int face,
+                   /*out*/ struct vector * n)
+{
+  struct point pobj;
+  struct vector nobj, tang_obj1, tang_obj2, tang_world1, tang_world2;
+
+  apply_to_point(obj->world2obj, p, &pobj);
+  normal_vector_object(obj, &pobj, face, &nobj);
+  tangent_vectors(&nobj, &tang_obj1, &tang_obj2);
+  apply_to_vect(obj->obj2world, &tang_obj1, &tang_world1);
+  apply_to_vect(obj->obj2world, &tang_obj2, &tang_world2);
+  product(&tang_world1, &tang_world2, n);
+  vnormalize(n, n);
+}
diff --git a/test/raytracer/object.h b/test/raytracer/object.h
new file mode 100644
index 000000000..d890bda86
--- /dev/null
+++ b/test/raytracer/object.h
@@ -0,0 +1,36 @@
+struct object {
+  enum { Cone, Cube, Cylinder, Plane, Sphere,
+         Union, Intersection, Difference } kind;
+  /* For base objects: Cone, Cube, Cylinder, Plane, Sphere */
+  struct closure surf;          /* surface function */
+  struct matrix * world2obj, * obj2world;
+  flt max_scale_applied;
+  /* For compound objects: Union, Intersection, Difference */
+  struct object * o1, * o2;
+  /* For all objects: bounding sphere (center and radius) */
+  struct point center;
+  flt radius;
+};
+
+struct object * cone(struct closure * c);
+struct object * cube(struct closure * c);
+struct object * cylinder(struct closure * c);
+struct object * plane(struct closure * c);
+struct object * sphere(struct closure * c);
+
+struct object * orotatex(struct object * o1, flt a);
+struct object * orotatey(struct object * o1, flt a);
+struct object * orotatez(struct object * o1, flt a);
+struct object * oscale(struct object * o1, flt sx, flt sy, flt sz);
+struct object * otranslate(struct object * o1,
+                           flt tx, flt ty, flt tz);
+struct object * ouscale(struct object * o1, flt s);
+
+struct object * odifference(struct object * o1, struct object * o2);
+struct object * ointersect(struct object * o1, struct object * o2);
+struct object * ounion(struct object * o1, struct object * o2);
+
+void normal_vector(struct object * obj, struct point * p, int face,
+                   /*out*/ struct vector * n);
+
+#define BS_NOT_COMPUTED (-1.0)
diff --git a/test/raytracer/point.h b/test/raytracer/point.h
new file mode 100644
index 000000000..b8d6ff0fc
--- /dev/null
+++ b/test/raytracer/point.h
@@ -0,0 +1,18 @@
+struct point {
+  flt x, y, z;
+};
+
+#if 0
+static inline flt dist2(struct point * p1, struct point * p2)
+{
+  flt dx = p2->x - p1->x;
+  flt dy = p2->y - p1->y;
+  flt dz = p2->z - p1->z;
+  return dx * dx + dy * dy + dz * dz;
+}
+#else
+#define dist2(p1,p2) \
+  (((p2)->x - (p1)->x) * ((p2)->x - (p1)->x) + \
+   ((p2)->y - (p1)->y) * ((p2)->y - (p1)->y) + \
+   ((p2)->z - (p1)->z) * ((p2)->z - (p1)->z))
+#endif
diff --git a/test/raytracer/render.c b/test/raytracer/render.c
new file mode 100644
index 000000000..3678cb6fa
--- /dev/null
+++ b/test/raytracer/render.c
@@ -0,0 +1,128 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+#include "matrix.h"
+#include "eval.h"
+#include "object.h"
+#include "light.h"
+#include "intersect.h"
+#include "surface.h"
+#include "simplify.h"
+#include "render.h"
+
+static void render_ray(struct point * amb,
+                       int depth,
+                       int initial,
+                       struct object * obj,
+                       struct light ** lights,
+                       struct point * p,
+                       struct vector * v,
+                       flt multx,
+                       flt multy,
+                       flt multz,
+                       /*out*/ struct point * color)
+{
+  struct object * bobj;
+  flt t;
+  struct point inter_w, inter_o;
+  int face;
+  flt surf_u, surf_v;
+  struct surface_characteristics sc;
+  flt dotprod;
+  struct vector n, n2, s;
+  struct point il, is;
+
+  if (depth < 0) {
+    color->x = color->y = color->z = 0.0;
+    return;
+  }
+  bobj = intersect_ray(p, v, obj, initial, &t);
+  if (bobj == NULL || t == 0.0) {
+    color->x = color->y = color->z = 0.0;
+    return;
+  }
+  /* Compute surface characteristics */
+  point_along(p, v, t, &inter_w);
+  apply_to_point(bobj->world2obj, &inter_w, &inter_o);
+  surface_coords(bobj, &inter_o, &face, &surf_u, &surf_v);
+  surface_function(&bobj->surf, face, surf_u, surf_v, &sc);
+  /* Construct the vectors on figure 4 */
+  normal_vector(bobj, &inter_w, face, &n);
+  dotprod = dotproduct(v, &n);
+  vscale(&n, 2.0 * dotprod, &n2);
+  vsub(v, &n2, &s);
+  if (dotprod > 0.0) opposite(&n, &n);
+  /* Light sources */
+  color_from_lights(obj, &inter_w, v, &n,
+                    sc.kd, sc.ks, sc.phong,
+                    lights, &il);
+  /* Recursive call for ray along s */
+  multx = multx * sc.ks * sc.x;
+  multy = multy * sc.ks * sc.y;
+  multz = multz * sc.ks * sc.z;
+  if (multx < 0.1 && multy < 0.1 && multz < 0.1)
+    is.x = is.y = is.z = 0.0;
+  else
+    render_ray(amb, depth - 1, 0, obj, lights, &inter_w, &s,
+               multx, multy, multz, &is);
+  /* Compute final color */
+  color->x = (sc.kd * amb->x + il.x + sc.ks * is.x) * sc.x;
+  color->y = (sc.kd * amb->y + il.y + sc.ks * is.y) * sc.y;
+  color->z = (sc.kd * amb->z + il.z + sc.ks * is.z) * sc.z;
+}
+
+static int convert_color(flt c)
+{
+  int n = (int) (c * 255.0);
+  if (n < 0) n = 0;
+  if (n > 255) n = 255;
+  return n;
+}
+
+void render(struct point * amb,
+            int numlights,
+            struct light ** lights,
+            struct object * scene,
+            int depth,
+            flt fov,
+            int wid,
+            int ht,
+            char * filename)
+{
+  FILE * oc;
+  flt wid2, ht2, scale, x, y;
+  int i, j;
+  struct point p;
+  struct vector v;
+  struct point color;
+  char * command;
+
+  compute_bounding_spheres(scene);
+  wid2 = (flt) wid / 2.0;
+  ht2 = (flt) ht / 2.0;
+  scale = tan(fov / 2.0) / wid2;
+  oc = fopen(filename, "w");
+  fprintf(oc, "P6\n# Camls 'R Us\n%d %d\n255\n", wid, ht);
+  arena_init();
+  for (j = ht - 1; j >= 0; j--) {
+    y = ((flt) j - ht2 + 0.5) * scale;
+    for (i = 0; i < wid; i++) {
+      x = ((flt) i - wid2 + 0.5) * scale;
+      p.x = p.y = 0.0; p.z = -1.0;
+      v.dx = x; v.dy = y; v.dz = 1.0;
+      render_ray(amb, depth, 1, scene, lights, &p, &v, 255.0, 255.0, 255.0,
+                 &color);
+      fputc(convert_color(color.x), oc);
+      fputc(convert_color(color.y), oc);
+      fputc(convert_color(color.z), oc);
+      arena_clear();
+    }
+  }
+  fclose(oc);
+#ifdef XV
+  command = malloc(strlen(filename) + 20);
+  sprintf(command, "xv %s &", filename);
+  system(command);
+  free(command);
+#endif
+}
diff --git a/test/raytracer/render.h b/test/raytracer/render.h
new file mode 100644
index 000000000..558e2523a
--- /dev/null
+++ b/test/raytracer/render.h
@@ -0,0 +1,9 @@
+void render(struct point * amb,
+            int numlights,
+            struct light ** lights,
+            struct object * scene,
+            int depth,
+            flt fov,
+            int wid,
+            int ht,
+            char * filename);
diff --git a/test/raytracer/simplify.c b/test/raytracer/simplify.c
new file mode 100644
index 000000000..7a4a54571
--- /dev/null
+++ b/test/raytracer/simplify.c
@@ -0,0 +1,177 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+#include "arrays.h"
+#include "matrix.h"
+#include "eval.h"
+#include "object.h"
+#include "simplify.h"
+
+#define INFINITE_RADIUS 1e300
+
+static flt cone_radius = 1.0;
+static flt cube_radius = 0.86602540378443859659; /* sqrt(3)/2 */
+static flt cylinder_radius = 1.11803398874989490253; /* sqrt(5)/2 */
+static flt sphere_radius = 1.0;
+
+static struct point cone_center = { 0.0, 1.0, 0.0 };
+static struct point cube_center = { 0.5, 0.5, 0.5 };
+static struct point cylinder_center = { 0.0, 0.5, 0.0 };
+static struct point sphere_center = { 0, 0, 0 };
+static struct point plane_center = { 0, 0, 0 };
+
+static struct point origin = { 0, 0, 0 };
+
+static inline void set_infinite(struct object * t)
+{
+  t->radius = INFINITE_RADIUS;
+}
+
+static inline void union_bs(struct object * t1, struct object * t2,
+                            struct object * obj)
+{
+  struct vector c1c2;
+  flt dd2, rr, rr2, d, alpha;
+
+  if (t1->radius >= INFINITE_RADIUS || t2->radius >= INFINITE_RADIUS) {
+    obj->radius = INFINITE_RADIUS;
+    return;
+  }
+  between(&t1->center, &t2->center, &c1c2);
+  dd2 = vlength2(&c1c2);
+  rr = t2->radius - t1->radius;
+  rr2 = rr * rr;
+  if (dd2 <= rr2) {
+    /* take the biggest sphere */
+    if (t1->radius <= t2->radius) {
+      ASSIGN(obj->center, t2->center);
+      obj->radius = t2->radius;
+      set_infinite(t2);
+    } else {
+      ASSIGN(obj->center, t1->center);
+      obj->radius = t1->radius;
+      set_infinite(t1);
+    }
+    return;
+  }
+  d = sqrt(dd2);
+  alpha = rr / (2 * d) + 0.5;
+  point_along(&t1->center, &c1c2, alpha, &obj->center);
+  obj->radius = (d + t1->radius + t2->radius) / 2;
+}
+
+static inline void intersection_bs(struct object * t1, struct object * t2,
+                                   struct object * obj)
+{
+  struct vector c1c2;
+  flt dd2, rr, rr2, rpr, rpr2, diff, d, te1, te2, te3, te4, te, alpha;
+
+  if (t1->radius >= INFINITE_RADIUS) {
+    ASSIGN(obj->center, t2->center);
+    obj->radius = t2->radius;
+    return;
+  }
+  if (t2->radius >= INFINITE_RADIUS) {
+    ASSIGN(obj->center, t1->center);
+    obj->radius = t1->radius;
+    return;
+  }
+  between(&t1->center, &t2->center, &c1c2);
+  dd2 = vlength2(&c1c2);
+  rr = t1->radius - t2->radius;
+  rr2 = rr * rr;
+  if (dd2 <= rr2) {
+    /* take the smallest sphere */
+    if (t2->radius <= t1->radius) {
+      ASSIGN(obj->center, t2->center);
+      obj->radius = t2->radius;
+      set_infinite(t2);
+    } else {
+      ASSIGN(obj->center, t1->center);
+      obj->radius = t1->radius;
+      set_infinite(t1);
+    }
+    return;
+  }
+  rpr = t1->radius + t2->radius;
+  rpr2 = rpr * rpr;
+  if (dd2 > rpr2) {
+    ASSIGN(obj->center, origin);
+    obj->radius = 0.0;
+    return;
+  }
+  diff = t1->radius * t1->radius - t2->radius * t2->radius;
+  if (dd2 <= diff) {
+    ASSIGN(obj->center, t2->center);
+    obj->radius = t2->radius;
+    set_infinite(t2);
+    return;
+  }
+  if (dd2 <= -diff) {
+    ASSIGN(obj->center, t1->center);
+    obj->radius = t1->radius;
+    set_infinite(t1);
+    return;
+  }
+  d = sqrt(dd2);
+  te1 = t1->radius + d + t2->radius;
+  te2 = t1->radius + d - t2->radius;
+  te3 = t2->radius + d - t1->radius;
+  te4 = t1->radius + t2->radius - d;
+  te = (sqrt (te1 * te2 * te3 * te4)) / (2 * d);
+  alpha =
+    (t1->radius * t1->radius - t2->radius * t2->radius) / (2 * dd2) + 0.5;
+  point_along(&t1->center, &c1c2, alpha, &obj->center);
+  obj->radius = te;
+}
+
+static inline void difference_bs(struct object * t1, struct object * t2,
+                                   struct object * obj)
+{
+  ASSIGN(obj->center, t1->center);
+  obj->radius = t1->radius;
+  set_infinite(t1);
+}
+
+void compute_bounding_spheres(struct object * obj)
+{
+  if (obj->radius >= 0.0) return; /* already computed */
+  switch (obj->kind) {
+  case Cone:
+    apply_to_point(obj->obj2world, &cone_center, &obj->center);
+    obj->radius = obj->max_scale_applied * cone_radius;
+    break;
+  case Cube:
+    apply_to_point(obj->obj2world, &cube_center, &obj->center);
+    obj->radius = obj->max_scale_applied * cube_radius;
+    break;
+  case Cylinder:
+    apply_to_point(obj->obj2world, &cylinder_center, &obj->center);
+    obj->radius = obj->max_scale_applied * cylinder_radius;
+    break;
+  case Plane:
+    ASSIGN(obj->center, plane_center);
+    obj->radius = INFINITE_RADIUS;
+    break;
+  case Sphere:
+    apply_to_point(obj->obj2world, &sphere_center, &obj->center);
+    obj->radius = obj->max_scale_applied * sphere_radius;
+    break;
+  case Union:
+    compute_bounding_spheres(obj->o1);
+    compute_bounding_spheres(obj->o2);
+    union_bs(obj->o1, obj->o2, obj);
+    break;
+  case Intersection:
+    compute_bounding_spheres(obj->o1);
+    compute_bounding_spheres(obj->o2);
+    intersection_bs(obj->o1, obj->o2, obj);
+    break;
+  case Difference:
+    compute_bounding_spheres(obj->o1);
+    compute_bounding_spheres(obj->o2);
+    difference_bs(obj->o1, obj->o2, obj);
+    break;
+  }
+}
+
diff --git a/test/raytracer/simplify.h b/test/raytracer/simplify.h
new file mode 100644
index 000000000..b2d179574
--- /dev/null
+++ b/test/raytracer/simplify.h
@@ -0,0 +1,3 @@
+/* Add bounding sphere info to the given object (recursively) */
+
+void compute_bounding_spheres(struct object * obj);
diff --git a/test/raytracer/surface.c b/test/raytracer/surface.c
new file mode 100644
index 000000000..0059db101
--- /dev/null
+++ b/test/raytracer/surface.c
@@ -0,0 +1,148 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+#include "eval.h"
+#include "object.h"
+#include "surface.h"
+
+/* Sphere:
+      x = sqrt(1 - y^2) sin(360u)
+      y = 2v - 1
+      z = sqrt(1 - y^2) cos(360u)
+   hence v = (y+1)/2
+     and u = atan2_turns(x, z) (atan2 in "number of turns")
+*/
+
+#define INV2PI (1.0 / (2.0 * M_PI))
+
+static inline flt atan2_turns(flt x, flt z)
+{
+  flt r = atan2(x, z) * INV2PI;
+  if (r < 0.0) r += 1.0;
+  return r;
+}
+
+static inline void sphere_coords(flt x, flt y, flt z,
+                                 int * face, flt * u, flt * v)
+{
+  *face = 0;
+  *u = atan2_turns(x, z);
+  *v = (y + 1.0) * 0.5;
+}
+
+/* Cube:
+     (x, y, 0)  ->  (0, x, y)
+     (x, y, 1)  ->  (1, x, y)
+     (0, y, z)  ->  (2, z, y)
+     (1, y, z)  ->  (3, z, y)
+     (x, 1, z)  ->  (4, x, z)
+     (x, 0, z)  ->  (5, x, z)
+   Watch out for rounding errors when determining which coordinate is 0 or 1;
+   see which one is closest to 0 or 1 */
+
+static inline void cube_coords(flt x, flt y, flt z,
+                               int * face, flt * u, flt * v)
+{
+  flt dists[6] =
+  { fabs(z), fabs(1 - z), fabs(x), fabs(1 - x), fabs(y), fabs(1 - y) };
+  flt min;
+  int f, i;
+
+  f = 0; min = dists[0];
+  for (i = 1; i < 6; i++)
+    if (dists[i] < min) { min = dists[i]; f = i; }
+  *face = f;
+  switch (f) {
+  case 0: case 1:
+    *u = x; *v = y; break;
+  case 2: case 3:
+    *u = z; *v = y; break;
+  case 4: case 5:
+    *u = x; *v = z; break;
+  }
+}
+
+/* Cylinder:
+     (x, 0, z)  ->  (2, u, v) where x = 2u-1 and z = 2v-1 hence
+     (x, 1, z)  ->  (1, x, z)   u = (x+1)/2 and v = (z+1)/2
+     (x, y, z)  ->  (0, u, v) where x = sin(360u)  v = y  z = cos(360u)
+                              hence u = atan2_turns(x, z)  and v = y
+*/
+
+static inline void cylinder_coords(flt x, flt y, flt z,
+                                   int * face, flt * u, flt * v)
+{
+  flt min, d;
+  int f;
+
+  min = y * y; f = 0;
+  d = (1 - y) * (1 - y);
+  if (d < min) { min = d; f = 1; }
+  d = fabs(x * x + z * z - 1);
+  if (d < min) { min = d; f = 2; }
+  *face = 2 - f;
+  if (f < 2) {
+    *u = (x + 1) * 0.5;
+    *v = (z + 1) * 0.5;
+  } else {
+    *u = atan2_turns(x, z);
+    *v = y;
+  }
+}
+
+/* Cone:
+     (x, y, z)  ->  (0, u, v)  where x = v sin 360u
+                                     y = v
+                                     z = v cos 360u
+                               hence u = atan2_turns(x, z) and v = y
+     (x, 1, z)  ->  (1, u, v)  where x = 2u-1 and z = 2v-1
+                               hence u = (x+1)/2 and v = (z+1)/2
+*/
+
+static inline void cone_coords(flt x, flt y, flt z,
+                               int * face, flt * u, flt * v)
+{
+  if ((1 - y) * (1 - y) < fabs(x * x + z * z - y * y)) {
+    *face = 1;
+    *u = (x + 1) * 0.5;
+    *v = (z + 1) * 0.5;
+  } else {
+    *face = 0;
+    *u = atan2_turns(x, z);
+    *v = y;
+  }
+}
+
+/* Plane */
+
+static inline void plane_coords(flt x, flt y, flt z,
+                                int * face, flt * u, flt * v)
+{
+  *face = 0;
+  *u = x;
+  *v = z;
+}
+
+/* All together */
+
+void surface_coords(struct object * obj, struct point * p,
+                    /*out*/ int * face,
+                    /*out*/ flt * u,
+                    /*out*/ flt * v)
+{
+  switch (obj->kind) {
+  case Cone:
+    cone_coords(p->x, p->y, p->z, face, u, v); break;
+  case Cube:
+    cube_coords(p->x, p->y, p->z, face, u, v); break;
+  case Cylinder:
+    cylinder_coords(p->x, p->y, p->z, face, u, v); break;
+  case Plane:
+    plane_coords(p->x, p->y, p->z, face, u, v); break;
+  case Sphere:
+    sphere_coords(p->x, p->y, p->z, face, u, v); break;
+  default:
+    assert(0);
+  }
+}
+
diff --git a/test/raytracer/surface.h b/test/raytracer/surface.h
new file mode 100644
index 000000000..ac128b277
--- /dev/null
+++ b/test/raytracer/surface.h
@@ -0,0 +1,7 @@
+/* Return the texture coordinates for the given point on the given
+   object.  Point in object coords. */
+
+void surface_coords(struct object * obj, struct point * p,
+                    /*out*/ int * face,
+                    /*out*/ flt * u,
+                    /*out*/ flt * v);
diff --git a/test/raytracer/vector.c b/test/raytracer/vector.c
new file mode 100644
index 000000000..226b61a87
--- /dev/null
+++ b/test/raytracer/vector.c
@@ -0,0 +1,74 @@
+#include "config.h"
+#include "point.h"
+#include "vector.h"
+
+flt dotproduct(struct vector * a, struct vector * b)
+{
+  return a->dx * b->dx + a->dy * b->dy + a->dz * b->dz;
+}
+
+void between(struct point * p, struct point * q,
+                           /*out*/ struct vector * v)
+{
+  v->dx = q->x - p->x;
+  v->dy = q->y - p->y;
+  v->dz = q->z - p->z;
+}
+
+void opposite(struct vector * v,
+              /*out*/ struct vector * w)
+{
+  w->dx = - v->dx;
+  w->dy = - v->dy;
+  w->dz = - v->dz;
+}
+
+void point_along(struct point * p, struct vector * v,
+                 flt ac,
+                 /*out*/ struct point * q)
+{
+  q->x = p->x + v->dx * ac;
+  q->y = p->y + v->dy * ac;
+  q->z = p->z + v->dz * ac;
+}
+
+void product(struct vector * a, struct vector * b,
+                           /*out*/ struct vector * v)
+{
+  v->dx = a->dy * b->dz - a->dz * b->dy;
+  v->dy = a->dz * b->dx - a->dx * b->dz;
+  v->dz = a->dx * b->dy - a->dy * b->dx;
+}
+
+flt vlength2(struct vector * a)
+{
+  return a->dx * a->dx + a->dy * a->dy + a->dz * a->dz;
+}
+
+flt vlength(struct vector * a)
+{
+  return sqrt(vlength2(a));
+}
+
+void vscale(struct vector * a, flt s,
+                          /*out*/ struct vector * v)
+{
+  v->dx = a->dx * s;
+  v->dy = a->dy * s;
+  v->dz = a->dz * s;
+}
+
+void vnormalize(struct vector * a, /*out*/ struct vector * v)
+{
+  vscale(a, 1 / vlength(a), v);
+}
+
+void vsub(struct vector * a, struct vector * b,
+                        /*out*/ struct vector * v)
+{
+  v->dx = a->dx - b->dx;
+  v->dy = a->dy - b->dy;
+  v->dz = a->dz - b->dz;
+}
+
+  
diff --git a/test/raytracer/vector.h b/test/raytracer/vector.h
new file mode 100644
index 000000000..c5e400ada
--- /dev/null
+++ b/test/raytracer/vector.h
@@ -0,0 +1,23 @@
+struct vector {
+  flt dx, dy, dz;
+};
+
+flt dotproduct(struct vector * a, struct vector * b);
+void between(struct point * p, struct point * q,
+             /*out*/ struct vector * v);
+void opposite(struct vector * v,
+              /*out*/ struct vector * w);
+void point_along(struct point * p, struct vector * v,
+                 flt ac,
+                 /*out*/ struct point * q);
+void product(struct vector * a, struct vector * b,
+             /*out*/ struct vector * v);
+flt vlength2(struct vector * a);
+flt vlength(struct vector * a);
+
+void vscale(struct vector * a, flt s,
+            /*out*/ struct vector * v);
+void vnormalize(struct vector * a, /*out*/ struct vector * v);
+void vsub(struct vector * a, struct vector * b,
+          /*out*/ struct vector * v);
+  
-- 
GitLab