Ajout nouveaux tests

git-svn-id: https://yquem.inria.fr/compcert/svn/compcert/trunk@708 fca1b0fc-160b-0410-b1d3-a4f43f01ea2e

32 files changed, 3031 insertions, 0 deletions

diff --git a/test/raytracer/.depend b/test/raytracer/.depend
new file mode 100644
index 0000000..c31fbcd
--- /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 0000000..cd442c5
--- /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 0000000..43508ee
--- /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 0000000..47634f6
--- /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 0000000..ed5b847
--- /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 0000000..4d3f3dc
--- /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 0000000..610b5d9
--- /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 0000000..20ca3e9
--- /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 0000000..332126b
--- /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 0000000..9b8d630
--- /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 0000000..b0a6cbc
--- /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 0000000..6a31546
--- /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 0000000..c8a2126
--- /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 0000000..35241a6
--- /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 0000000..e9d5815
--- /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 0000000..1c90104
--- /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 0000000..9a1d3e2
--- /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 0000000..1efd396
--- /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 0000000..8814346
--- /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 0000000..4b9f434
--- /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 0000000..fe194cc
--- /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 0000000..3dd6e77
--- /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 0000000..d890bda
--- /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 0000000..b8d6ff0
--- /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 0000000..3678cb6
--- /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 0000000..558e252
--- /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 0000000..7a4a545
--- /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 0000000..b2d1795
--- /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 0000000..0059db1
--- /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 0000000..ac128b2
--- /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 0000000..226b61a
--- /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 0000000..c5e400a
--- /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);
+