Fix perspective correctness, add -rgb, -grayscale, w parameter for vertex

git-svn-id: http://svn.code.sf.net/p/netpbm/code/trunk@3376 9d0c8265-081b-0410-96cb-a4ca84ce46f8

13 files changed, 431 insertions, 571 deletions

diff --git a/generator/pamtris/boundaries.c b/generator/pamtris/boundaries.c
index 8ea28682..7045cbc7 100644
--- a/generator/pamtris/boundaries.c
+++ b/generator/pamtris/boundaries.c
@@ -3,11 +3,11 @@
 ===============================================================================
    Boundary buffer functions
 
-   New triangles are drawn one scanline at a time, and for every such scanline
-   we have left and right boundary columns within the frame buffer such that
-   the fraction of the triangle's area within that scanline is enclosed
-   between those two points (inclusive). Those coordinates may correspond to
-   columns outside the frame buffer's actual limits, in which case proper
+   New triangles are drawn one row at a time, and for every such row we have
+   left and right boundary columns within the frame buffer such that the
+   fraction of the triangle's area within that scanline is enclosed between
+   those two points (inclusive). Those coordinates may correspond to columns
+   outside the frame buffer's actual limits, in which case proper
    post-processing should be made wherever such coordinates are used to
    actually plot anything into the frame buffer.
 =============================================================================*/
@@ -17,38 +17,24 @@
 #include <netpbm/mallocvar.h>
 #include <netpbm/pm.h>
 
+#include "varying.h"
 #include "utils.h"
-#include "fract.h"
 
 
 #include "boundaries.h"
 
 
 
-static fract
-make_pos_fract(int32_t const quotient,
-               int32_t const remainder) {
-
-    fract retval;
-
-    retval.q = quotient;
-    retval.r = remainder;
-    retval.negative_flag = 0;
-
-    return retval;
-}
-
-
-
 void
 init_boundary_buffer(boundary_info * const bi,
                      int16_t         const height) {
 
     MALLOCARRAY(bi->buffer, height * 2);
 
-    if (!bi->buffer)
-        pm_error("Unable to get memory for %u-row high boundary buffer",
+    if (!bi->buffer) {
+        pm_error("unable to get memory for %u-row high boundary buffer.",
                  height);
+    }
 }
 
 
@@ -68,39 +54,30 @@ gen_triangle_boundaries(Xy              const xy,
 /*----------------------------------------------------------------------------
   Generate an entry in the boundary buffer for the boundaries of every
   VISIBLE row of a particular triangle. In case there is no such row,
-  start_row is accordingly set to -1. The argument is a 3-element array
-  of pairs of int16_t's representing the coordinates of the vertices of
+  start_scanline is accordingly set to -1. "xy" is a 3-element array
+  of pairs of integers representing the coordinates of the vertices of
   a triangle. Those vertices MUST be already sorted in order from the
   uppermost to the lowermost vertex (which is what draw_triangle, the
   only function which uses this one, does with the help of sort3).
 
-  The return value indicates whether the middle vertex is to the left of the
-  line connecting the top vertex to the bottom vertex or not.
+  The return value indicates whether the middle vertex is to the left of
+  the line connecting the top vertex to the bottom vertex or not.
 -----------------------------------------------------------------------------*/
     int16_t leftmost_x;
     int16_t rightmost_x;
     int mid_is_to_the_left;
-    fract left_x;
-    fract right_x;
-    bool no_upper_part;
-    int32_t top2mid_delta;
-    int32_t top2bot_delta;
-    int32_t mid2bot_delta;
-    fract top2mid_step;
-    fract top2bot_step;
-    fract mid2bot_step;
-    fract* upper_left_step;
-    fract* lower_left_step;
-    fract* upper_right_step;
-    fract* lower_right_step;
-    int32_t upper_left_delta;
-    int32_t lower_left_delta;
-    int32_t upper_right_delta;
-    int32_t lower_right_delta;
-    fract* left_step[2];
-    fract* right_step[2];
-    int32_t left_delta[2];
-    int32_t right_delta[2];
+    varying top_x;
+    varying mid_x;
+    varying bot_x;
+    varying top2mid;
+    varying top2bot;
+    varying mid2bot;
+    varying* upper_left;
+    varying* lower_left;
+    varying* upper_right;
+    varying* lower_right;
+    varying* left[2];
+    varying* right[2];
     int16_t* num_rows_ptr[2];
     int32_t y;
     int32_t i;
@@ -114,8 +91,8 @@ gen_triangle_boundaries(Xy              const xy,
     bi->num_lower_rows = 0;
 
     if (xy._[2][1] < 0 || xy._[0][1] >= height) {
-        /* Triangle is either completely above the topmost scanline or
-           completely below the bottom scanline.
+        /* Triangle is either completely above the uppermost scanline or
+           completely below the lowermost scanline.
         */
 
         return false; /* Actual value doesn't matter. */
@@ -154,106 +131,57 @@ gen_triangle_boundaries(Xy              const xy,
 
     mid_is_to_the_left = 2;
 
-    left_x  = make_pos_fract(xy._[0][0], 0);
-    right_x = make_pos_fract(xy._[0][0], 0);
+    int32_to_varying_array(&xy._[0][0], &top_x, 1);
+    int32_to_varying_array(&xy._[1][0], &mid_x, 1);
+    int32_to_varying_array(&xy._[2][0], &bot_x, 1);
 
     if (xy._[0][1] == xy._[1][1]) {
         /* Triangle has only a lower part. */
+        k = 1;
 
         mid_is_to_the_left = 0;
+    } else {
+        k = 0;
 
-        right_x.q = xy._[1][0];
-    } else if (xy._[1][1] == xy._[2][1]) {
-        /* Triangle has only an upper part (plus the row of the middle
-           vertex).
-        */
-
-        mid_is_to_the_left = 1;
+        if (xy._[1][1] == xy._[2][1]) {
+            /* Triangle has only an upper part (plus the row of the middle
+               vertex).
+            */
+            mid_is_to_the_left = 1;
+        }
     }
 
-    no_upper_part = (xy._[1][1] == xy._[0][1]);
-
-    top2mid_delta = xy._[1][1] - xy._[0][1] + !no_upper_part;
-    top2bot_delta = xy._[2][1] - xy._[0][1] + 1;
-    mid2bot_delta = xy._[2][1] - xy._[1][1] + no_upper_part;
-
-    gen_steps(&xy._[0][0], &xy._[1][0], &top2mid_step, 1, top2mid_delta);
-    gen_steps(&xy._[0][0], &xy._[2][0], &top2bot_step, 1, top2bot_delta);
-    gen_steps(&xy._[1][0], &xy._[2][0], &mid2bot_step, 1, mid2bot_delta);
+    prepare_for_interpolation(&top_x, &mid_x, &top2mid, xy._[1][1] - xy._[0][1], 1);
+    prepare_for_interpolation(&top_x, &bot_x, &top2bot, xy._[2][1] - xy._[0][1], 1);
+    prepare_for_interpolation(&mid_x, &bot_x, &mid2bot, xy._[2][1] - xy._[1][1], 1);
 
     if (mid_is_to_the_left == 2) {
-        if (top2bot_step.negative_flag) {
-            if (top2mid_step.negative_flag) {
-                if (top2mid_step.q == top2bot_step.q) {
-                    mid_is_to_the_left =
-                        top2mid_step.r * top2bot_delta >
-                        top2bot_step.r * top2mid_delta;
-                } else {
-                    mid_is_to_the_left = top2mid_step.q < top2bot_step.q;
-                }
-            } else {
-                mid_is_to_the_left = 0;
-            }
-        } else {
-            if (!top2mid_step.negative_flag) {
-                if (top2mid_step.q == top2bot_step.q) {
-                    mid_is_to_the_left =
-                        top2mid_step.r * top2bot_delta <
-                        top2bot_step.r * top2mid_delta;
-                } else {
-                    mid_is_to_the_left = top2mid_step.q < top2bot_step.q;
-                }
-            } else {
-                mid_is_to_the_left = 1;
-            }
-        }
+        mid_is_to_the_left = top2mid.s < top2bot.s;
     }
+
     if (mid_is_to_the_left) {
-        upper_left_step     = &top2mid_step;
-        lower_left_step     = &mid2bot_step;
-        upper_right_step    = &top2bot_step;
-        lower_right_step    = upper_right_step;
-
-        upper_left_delta    = top2mid_delta;
-        lower_left_delta    = mid2bot_delta;
-        upper_right_delta   = top2bot_delta;
-        lower_right_delta   = upper_right_delta;
+        upper_left     = &top2mid;
+        lower_left     = &mid2bot;
+        upper_right    = &top2bot;
+        lower_right    = upper_right;
     } else {
-        upper_right_step    = &top2mid_step;
-        lower_right_step    = &mid2bot_step;
-        upper_left_step     = &top2bot_step;
-        lower_left_step     = upper_left_step;
-
-        upper_right_delta   = top2mid_delta;
-        lower_right_delta   = mid2bot_delta;
-        upper_left_delta    = top2bot_delta;
-        lower_left_delta    = upper_left_delta;
+        upper_right    = &top2mid;
+        lower_right    = &mid2bot;
+        upper_left     = &top2bot;
+        lower_left     = upper_left;
     }
 
-    left_step[0] = upper_left_step;
-    left_step[1] = lower_left_step;
-    right_step[0] = upper_right_step;
-    right_step[1] = lower_right_step;
-    left_delta[0] = upper_left_delta;
-    left_delta[1] = lower_left_delta;
-    right_delta[0] = upper_right_delta;
-    right_delta[1] = lower_right_delta;
+    left[0] = upper_left;
+    left[1] = lower_left;
+    right[0] = upper_right;
+    right[1] = lower_right;
+
     num_rows_ptr[0] = &bi->num_upper_rows;
     num_rows_ptr[1] = &bi->num_lower_rows;
 
     y = xy._[0][1];
 
     i = 0;
-    k = 0;
-
-    if (no_upper_part) {
-        k = 1;
-
-        right_x.q = xy._[1][0];
-    }
-
-    step_up(&left_x, left_step[k], 1, left_delta[k]);
-    step_up(&right_x, right_step[k], 1, right_delta[k]);
 
     while (k < 2) {
         int32_t end;
@@ -271,8 +199,8 @@ gen_triangle_boundaries(Xy              const xy,
 
             y += delta;
 
-            multi_step_up(&left_x, left_step[k], 1, delta, left_delta[k]);
-            multi_step_up(&right_x, right_step[k], 1, delta, right_delta[k]);
+            multi_step_up(left[k], delta, 1);
+            multi_step_up(right[k], delta, 1);
 
             if (y < 0) {
                 k++;
@@ -287,7 +215,7 @@ gen_triangle_boundaries(Xy              const xy,
         }
 
         while (y < end) {
-            if (left_x.q >= width || right_x.q < 0) {
+            if (round_varying(*left[k]) >= width || round_varying(*right[k]) < 0) {
                 if (bi->start_scanline > -1) {
                     return mid_is_to_the_left;
                 }
@@ -296,14 +224,14 @@ gen_triangle_boundaries(Xy              const xy,
                     bi->start_scanline = y;
                 }
 
-                bi->buffer[i++] = left_x.q;
-                bi->buffer[i++] = right_x.q;
+                bi->buffer[i++] = round_varying(*left[k]);
+                bi->buffer[i++] = round_varying(*right[k]);
 
                 (*(num_rows_ptr[k]))++;
             }
 
-            step_up(&left_x, left_step[k], 1, left_delta[k]);
-            step_up(&right_x, right_step[k], 1, right_delta[k]);
+            step_up(left[k], 1);
+            step_up(right[k], 1);
 
             y++;
         }
diff --git a/generator/pamtris/boundaries.h b/generator/pamtris/boundaries.h
index d41719cb..70f7f90d 100644
--- a/generator/pamtris/boundaries.h
+++ b/generator/pamtris/boundaries.h
@@ -1,6 +1,7 @@
 #ifndef BOUNDARIES_H_INCLUDED
 #define BOUNDARIES_H_INCLUDED
 
+#include <stdbool.h>
 #include <stdint.h>
 
 #include "triangle.h"
diff --git a/generator/pamtris/fract.h b/generator/pamtris/fract.h
deleted file mode 100644
index ff3c4402..00000000
--- a/generator/pamtris/fract.h
+++ /dev/null
@@ -1,54 +0,0 @@
-#ifndef FRACT_H_INCLUDED
-#define FRACT_H_INCLUDED
-
-#include <stdbool.h>
-#include <stdint.h>
-
-
-typedef struct {
-/*----------------------------------------------------------------------------
-    This struct and the functions that manipulate variables of this type act
-    as a substitute for floating point computations. Here, whenever we need a
-    value with a fractional component, we represent it using two parts: 1. An
-    integer part, called the "quotient", and 2. A fractional part, which is
-    itself composed of a "remainder" (or "numerator") and a "divisor" (or
-    "denominator"). The fract struct provides storage for the quotient and the
-    remainder, but the divisor must be given separately (because it often
-    happens in this program that whenever we are dealing with one variable of
-    type fract, we are dealing with more of them at the same time, and they
-    all have the same divisor).
-
-    To be more precise, the way we actually use variables of this type works
-    like this: We read integer values through standard input; When drawing
-    triangles, we need need to calculate differences between some pairs of
-    these input values and divide such differences by some other integer,
-    which is the above mentioned divisor. That result is then used to compute
-    successive interpolations between the two values for which we had
-    originally calculated the difference, and is therefore called the
-    "interpolation step". The values between which we wish to take successive
-    interpolations are called the "initial value" and the "final value". The
-    interpolation procedure works like this: First, we transform the initial
-    value into a fract variable by equating the quotient of that variable to
-    the initial value and assigning 0 to its remainder. Then, we successivelly
-    apply the interpolation step to that variable through successive calls to
-    step_up and/or multi_step_up until the quotient of the variable equals the
-    final value. Each application of step_up or multi_step_up yields a
-    particular linear interpolation between the initial and final values.
-
-    If and only if a particular fract variable represents an interpolation
-    step, the "negative_flag" field indicates whether the step is negative
-    (i. e. negative_flag == true) or not (negative_flag == false). This is
-    necessary in order to make sure that variables are "stepped up" in the
-    appropriate direction, so to speak, as the field which stores the
-    remainder in any fract variable, "r", is always equal to or above 0, and
-    the quotient of a step may be 0, so the actual sign of the step value is
-    not always discoverable through a simple examination of the sign of the
-    quotient. On the other hand, if the variable does not represent an
-    interpolation step, the negative_flag is meaningless.
------------------------------------------------------------------------------*/
-    int32_t q;     /* Quotient */
-    int32_t r: 31; /* Remainder */
-    bool    negative_flag: 1;
-} fract;
-
-#endif
diff --git a/generator/pamtris/framebuffer.c b/generator/pamtris/framebuffer.c
index 03cd720c..93263c91 100644
--- a/generator/pamtris/framebuffer.c
+++ b/generator/pamtris/framebuffer.c
@@ -37,9 +37,10 @@
 #include <stdint.h>
 #include <stdlib.h>
 #include <string.h>
+#include <math.h>
 
 #include "utils.h"
-#include "fract.h"
+#include "varying.h"
 #include "limits_pamtris.h"
 
 #include "framebuffer.h"
@@ -60,8 +61,6 @@ set_tupletype(const char * const str,
 -----------------------------------------------------------------------------*/
     if (str == NULL) {
         memset(tupletype, 0, 256);
-
-        return 1;
     } else {
         size_t len;
 
@@ -80,9 +79,9 @@ set_tupletype(const char * const str,
         while(len > 0 && isspace(tupletype[len])) {
             tupletype[len--] = '\0';
         }
-
-        return 1;
     }
+
+    return 1;
 }
 
 
@@ -279,9 +278,7 @@ clear_framebuffer(bool               const clear_image_buffer,
 void
 draw_span(uint32_t           const base,
           uint16_t           const length,
-          fract *            const attribs_start,
-          const fract *      const attribs_steps,
-          int32_t            const div,
+          varying *          const attribs,
           framebuffer_info * const fbi) {
 /*----------------------------------------------------------------------------
   Draw a horizontal span of "length" pixels into the frame buffer, performing
@@ -291,42 +288,50 @@ draw_span(uint32_t           const base,
 
   This function does not perform any kind of bounds checking.
 -----------------------------------------------------------------------------*/
-    uint8_t const num_planes = fbi->num_attribs + 1;
+    static double const depth_range = MAX_Z;
+
+    uint16_t const maxval = fbi->maxval;
+    uint8_t  const z      = fbi->num_attribs;
+    uint8_t  const w      = z + 1;
+    uint8_t  const n      = w + 1;
+
+    uint8_t  const num_planes = w;
 
     unsigned int i;
 
     /* Process each pixel in the span: */
 
     for (i = 0; i < length; i++) {
-        int32_t  const z        = MAX_Z - attribs_start[fbi->num_attribs].q;
-        uint32_t const z_mask   = -(~(z - fbi->z.buffer[base + i]) >> 31);
-        uint32_t const n_z_mask = ~z_mask;
+        int32_t  const d      = round(depth_range * attribs[z].v);
+        uint32_t const d_mask = geq_mask64(d, fbi->z.buffer[base + i]);
 
         uint32_t const j = base + i;
         uint32_t const k = j * num_planes;
 
+        varying const inverse_w = inverse_varying(attribs[w]);
+
         unsigned int l;
 
         /* The following statements will only have any effect if the depth
            test, performed above, has suceeded. I. e. if the depth test fails,
-           no changes will be made on the framebuffer; otherwise, the
-           framebuffer will be updated with the new values.
+           no changes will be made on the frame buffer; otherwise, the
+           frame buffer will be updated with the new values.
         */
-        fbi->z.buffer[j] = (fbi->z.buffer[j] & n_z_mask) | (z & z_mask);
+        fbi->z.buffer[j] = (fbi->z.buffer[j] & ~d_mask) | (d & d_mask);
+
+        for (l = 0; l < z; l++) {
+	    varying const newval = multiply_varyings(attribs[l], inverse_w);
 
-        for (l = 0; l < fbi->num_attribs; l++) {
             fbi->img.buffer[k + l] =
-                (fbi->img.buffer[k + l] & n_z_mask) |
-                (attribs_start[l].q & z_mask);
+                (fbi->img.buffer[k + l] & ~d_mask) |
+                (round_varying(newval) &  d_mask);
         }
 
-        fbi->img.buffer[k + fbi->num_attribs] =
-            (fbi->img.buffer[k + fbi->num_attribs] & n_z_mask) |
-            (fbi->maxval & z_mask);
+        fbi->img.buffer[k + z] |= (maxval & d_mask);
 
         /* Compute the attribute values for the next pixel: */
 
-        step_up(attribs_start, attribs_steps, num_planes, div);
+        step_up(attribs, n);
     }
 }
 
diff --git a/generator/pamtris/framebuffer.h b/generator/pamtris/framebuffer.h
index 73ec96be..b3d4f7a3 100644
--- a/generator/pamtris/framebuffer.h
+++ b/generator/pamtris/framebuffer.h
@@ -3,7 +3,8 @@
 
 #include <stdint.h>
 #include <stdbool.h>
-#include "fract.h"
+
+#include "varying.h"
 #include "netpbm/pam.h"
 
 typedef struct framebuffer_info {
@@ -12,7 +13,8 @@ typedef struct framebuffer_info {
 -----------------------------------------------------------------------------*/
     /* These fields are initialized once by reading the command line
        arguments. "maxval" and "num_attribs" may be modified later
-       through "realloc_image_buffer".
+       through "realloc_image_buffer"; "correct" may also be modified
+       if the eponymous command is given.
     */
     int32_t width;
     int32_t height;
@@ -67,9 +69,7 @@ realloc_image_buffer(int32_t            const new_maxval,
 void
 draw_span(uint32_t           const base,
           uint16_t           const length,
-          fract *            const attribs_start,
-          const fract *      const attribs_steps,
-          int32_t            const divisor,
+          varying *          const attribs,
           framebuffer_info * const fbi);
 
 #endif
diff --git a/generator/pamtris/input.c b/generator/pamtris/input.c
index 166d6db5..c35359b4 100644
--- a/generator/pamtris/input.c
+++ b/generator/pamtris/input.c
@@ -17,37 +17,32 @@
 
 #include "input.h"
 
-#define MAX_COORD       32767
-#define MIN_COORD       -MAX_COORD
-
 #define DRAW_MODE_TRIANGLES 1
 #define DRAW_MODE_STRIP     2
 #define DRAW_MODE_FAN       3
 
 #define CMD_SET_MODE        "mode"
 #define CMD_SET_ATTRIBS     "attribs"
-#define CMD_VERTEX      "vertex"
-#define CMD_PRINT       "print"
-#define CMD_CLEAR       "clear"
-#define CMD_RESET       "reset"
-#define CMD_QUIT        "quit"
+#define CMD_VERTEX          "vertex"
+#define CMD_PRINT           "print"
+#define CMD_CLEAR           "clear"
+#define CMD_RESET           "reset"
+#define CMD_QUIT            "quit"
 
 #define ARG_TRIANGLES       "triangles"
-#define ARG_STRIP       "strip"
-#define ARG_FAN         "fan"
-#define ARG_IMAGE       "image"
-#define ARG_DEPTH       "depth"
+#define ARG_STRIP           "strip"
+#define ARG_FAN             "fan"
+#define ARG_IMAGE           "image"
+#define ARG_DEPTH           "depth"
 
 #define WARNING_EXCESS_ARGS "warning: ignoring excess arguments: line %lu."
 #define SYNTAX_ERROR        "syntax error: line %lu."
-
+    
 typedef struct {
     Xy v_xy;
         /* X- and Y-coordinates of the vertices for the current triangle.
-           int32_t v_attribs[3][MAX_NUM_ATTRIBS + 1]; // Vertex attributes for
-           the current triangle. Includes the Z-coordinates.
         */
-	Attribs v_attribs;
+    Attribs v_attribs;
         /* Vertex attributes for the current triangle. Includes the
            Z-coordinates.
         */
@@ -371,20 +366,32 @@ process_next_command(input_info           * const line,
             break;
         }
 
-        for (i = 0; i < 3; i++) {
+        for (i = 0; i < 4; i++) {
             nt = next_token(nt.end);
 
             i_args[i] = strtol(nt.begin, &strtol_end, 10);
 
-            if (*nt.begin == '\0' || strtol_end != nt.end) {
-                pm_errormsg(SYNTAX_ERROR, line->number);
+            if (*nt.begin == '\0') {
+		if(i != 3) {
+                    pm_errormsg(SYNTAX_ERROR, line->number);
 
-                must_break_out = true;
+                    must_break_out = true;
 
-                break;
+                    break;
+                } else {
+                    i_args[i] = 1;
+                }
+            } else {
+                if (strtol_end != nt.end) {
+                    pm_errormsg(SYNTAX_ERROR, line->number);
+
+                    must_break_out = true;
+
+                    break;
+                }
             }
 
-            if (i < 2) {
+            if (i < 3) {
                 if (i_args[i] < MIN_COORD || i_args[i] > MAX_COORD) {
                     pm_errormsg(
                         "error: coordinates out of bounds: line %lu.",
@@ -395,9 +402,9 @@ process_next_command(input_info           * const line,
                     break;
                 }
             } else {
-                if (i_args[i] < 0 || i_args[i] > MAX_Z) {
+                if (i_args[i] < MIN_INPUT_W || i_args[i] > MAX_INPUT_W) {
                     pm_errormsg(
-                        "error: Z component out of bounds: line %lu.",
+                        "error: perspective correction factor (w) out of bounds: line %lu.",
                         line->number);
 
                     must_break_out = true;
@@ -416,7 +423,8 @@ process_next_command(input_info           * const line,
             state.v_attribs._[state.next][i] = state.curr_attribs[i];
         }
 
-        state.v_attribs._[state.next][fbi->num_attribs] = i_args[2];
+        state.v_attribs._[state.next][fbi->num_attribs + 0] = i_args[2];
+        state.v_attribs._[state.next][fbi->num_attribs + 1] = i_args[3];
 
         state.v_xy._[state.next][0] = i_args[0];
         state.v_xy._[state.next][1] = i_args[1];
diff --git a/generator/pamtris/limits_pamtris.h b/generator/pamtris/limits_pamtris.h
index dcf1f1e6..a7ed503f 100644
--- a/generator/pamtris/limits_pamtris.h
+++ b/generator/pamtris/limits_pamtris.h
@@ -1,7 +1,11 @@
 #ifndef LIMITS_H_INCLUDED
 #define LIMITS_H_INCLUDED
 
-#define MAX_NUM_ATTRIBS     20
-#define MAX_Z           ((1 << 30) - 1)
+#define MAX_NUM_ATTRIBS 20
+#define MAX_COORD       32767
+#define MIN_COORD       (-MAX_COORD)
+#define MAX_INPUT_W     1048575
+#define MIN_INPUT_W     1
+#define MAX_Z           0x3FFFFFFF
 
 #endif
diff --git a/generator/pamtris/pamtris.c b/generator/pamtris/pamtris.c
index 74663531..d31e07a6 100644
--- a/generator/pamtris/pamtris.c
+++ b/generator/pamtris/pamtris.c
@@ -29,19 +29,21 @@ parse_command_line(int *         const argc_ptr,
         /* Instructions to pm_optParseOptions3 on how to parse our options */
     unsigned int option_def_index;
 
-    char * tupletype_ptr;
+    char * tupletype_tmp;
 
-    unsigned int width_spec, height_spec, maxval_spec, attribs_spec;
-    unsigned int tupletype_spec;
+    unsigned int width_spec, height_spec, attribs_spec, tupletype_spec;
+    unsigned int rgb_spec, grayscale_spec, maxval_spec;
 
     MALLOCARRAY_NOFAIL(option_def, 100);
 
     option_def_index = 0;  /* incremented by OPTENT3 */
     OPTENT3(0, "width",       OPT_INT,    width,          &width_spec,      0);
     OPTENT3(0, "height",      OPT_INT,    height,         &height_spec,     0);
-    OPTENT3(0, "maxval",      OPT_INT,    maxval,         &maxval_spec,     0);
     OPTENT3(0, "num_attribs", OPT_INT,    num_attribs,    &attribs_spec,    0);
-    OPTENT3(0, "tupletype",   OPT_STRING, &tupletype_ptr, &tupletype_spec,  0);
+    OPTENT3(0, "tupletype",   OPT_STRING, &tupletype_tmp, &tupletype_spec,  0);
+    OPTENT3(0, "rgb",         OPT_FLAG,   NULL,           &rgb_spec,        0);
+    OPTENT3(0, "grayscale",   OPT_FLAG,   NULL,           &grayscale_spec,  0);
+    OPTENT3(0, "maxval",      OPT_INT,    maxval,         &maxval_spec,     0);
 
     opt.opt_table     = option_def;
     opt.short_allowed = false;
@@ -49,9 +51,17 @@ parse_command_line(int *         const argc_ptr,
 
     pm_optParseOptions3(argc_ptr, (char **)argv, opt, sizeof(opt), 0);
 
-    if (!width_spec || !height_spec || !attribs_spec) {
+    if (!width_spec || !height_spec || (!attribs_spec && !(rgb_spec || grayscale_spec))) {
         pm_errormsg(
-            "you must at least specify -width, -height and -num_attribs.");
+            "you must at least specify -width, -height and "
+	    "either -num_attribs, -rgb or -grayscale.");
+
+        return 0;
+    }
+
+    if (rgb_spec + grayscale_spec + attribs_spec != 1) {
+        pm_errormsg("you must provide either only -num_attribs, "
+                    "-rgb or -grayscale; not a combination of those.");
 
         return 0;
     }
@@ -70,8 +80,6 @@ parse_command_line(int *         const argc_ptr,
 
     if (maxval_spec) {
         if (*maxval < 1 || *maxval > PAM_OVERALL_MAXVAL) {
-
-
             pm_errormsg("invalid maxval.");
 
             return 0;
@@ -80,6 +88,16 @@ parse_command_line(int *         const argc_ptr,
         *maxval = 255;
     }
 
+    if (rgb_spec) {
+        *num_attribs = 3;
+        set_tupletype("RGB_ALPHA", tupletype);
+    }
+
+    if (grayscale_spec) {
+        *num_attribs = 1;
+        set_tupletype("GRAYSCALE_ALPHA", tupletype);
+    }
+
     if (*num_attribs < 1 || *num_attribs > MAX_NUM_ATTRIBS) {
         pm_errormsg("invalid number of generic attributes per vertex.");
 
@@ -87,8 +105,15 @@ parse_command_line(int *         const argc_ptr,
     }
 
     if (tupletype_spec) {
-        if (!set_tupletype(tupletype_ptr, tupletype)) {
-            pm_errormsg("warning: invalid tuple type; using the null string.");
+        if(rgb_spec || grayscale_spec) {
+            pm_errormsg("you may not provide -tupletype together with "
+                        "-rgb or -grayscale.");
+
+            return 0;
+        }
+
+        if (!set_tupletype(tupletype_tmp, tupletype)) {
+            pm_errormsg("warning: invalid tuple type; using empty string.");
 
             set_tupletype(NULL, tupletype);
         }
@@ -112,9 +137,13 @@ main(int argc, const char ** argv) {
 
     set_tupletype(NULL, fbi.outpam.tuple_type);
 
-    if (!parse_command_line(&argc, argv,
-                            &fbi.width, &fbi.height, &fbi.maxval,
-                            &fbi.num_attribs, fbi.outpam.tuple_type)) {
+    if (!parse_command_line(&argc,
+                            argv,
+                            &fbi.width,
+                            &fbi.height,
+                            &fbi.maxval,
+                            &fbi.num_attribs,
+                            fbi.outpam.tuple_type)) {
         return 1;
     }
 
diff --git a/generator/pamtris/triangle.c b/generator/pamtris/triangle.c
index 09d821e0..5143f9ee 100644
--- a/generator/pamtris/triangle.c
+++ b/generator/pamtris/triangle.c
@@ -9,7 +9,7 @@
 #include "netpbm/mallocvar.h"
 
 #include "utils.h"
-#include "fract.h"
+#include "varying.h"
 #include "boundaries.h"
 #include "framebuffer.h"
 
@@ -17,29 +17,30 @@
 
 static void
 draw_partial_triangle(
-    const fract *         const left_attribs_input,
-    const fract *         const left_attribs_steps,
-    const fract *         const rght_attribs_input,
-    const fract *         const rght_attribs_steps,
-    int32_t               const left_div,
-    int32_t               const rght_div,
+    const varying *       const left_attribs_input,
+    const varying *       const rght_attribs_input,
     bool                  const upper_part,
     const boundary_info * const bi,
     framebuffer_info *    const fbi) {
 
-    uint8_t const num_planes = fbi->num_attribs + 1;
+    uint8_t const z = fbi->num_attribs;
+    uint8_t const w = z + 1;
+    uint8_t const n = w + 1;
 
-    fract * left_attribs;
-    fract * rght_attribs;
+    varying * left_attribs;
+    varying * rght_attribs;
+
+    varying * attribs;
 
     int32_t first_row_index;
     int32_t last_row_index;
 
-    MALLOCARRAY_NOFAIL(left_attribs, num_planes);
-    MALLOCARRAY_NOFAIL(rght_attribs, num_planes);
+    MALLOCARRAY_NOFAIL(left_attribs, n);
+    MALLOCARRAY_NOFAIL(rght_attribs, n);
+    MALLOCARRAY_NOFAIL(attribs, n);
 
-    memcpy(left_attribs, left_attribs_input, num_planes * sizeof(fract));
-    memcpy(rght_attribs, rght_attribs_input, num_planes * sizeof(fract));
+    memcpy(left_attribs, left_attribs_input, n * sizeof(varying));
+    memcpy(rght_attribs, rght_attribs_input, n * sizeof(varying));
 
     if (upper_part) {
         first_row_index = 0;
@@ -57,41 +58,26 @@ draw_partial_triangle(
         int32_t left_boundary;
         int32_t rght_boundary;
 
-        for (row = first_row_index; row <= last_row_index; ) {
+        for (row = first_row_index; row <= last_row_index; row++) {
             get_triangle_boundaries(row, &left_boundary, &rght_boundary, bi);
             {
                 int32_t const column_delta = rght_boundary - left_boundary;
                 int32_t start_column;
                 int32_t span_length;
 
-                fract   * attribs_start;
-                int32_t * attribs_begin;
-                int32_t * attribs_end;
-                fract   * attribs_steps;
-
-                MALLOCARRAY_NOFAIL(attribs_start, num_planes);
-                MALLOCARRAY_NOFAIL(attribs_begin, num_planes);
-                MALLOCARRAY_NOFAIL(attribs_end,   num_planes);
-                MALLOCARRAY_NOFAIL(attribs_steps, num_planes);
-
                 start_column = left_boundary;  /* initial value */
                 span_length = column_delta;    /* initial value */
 
-                fract_to_int32_array(left_attribs, attribs_begin, num_planes);
-                fract_to_int32_array(rght_attribs, attribs_end, num_planes);
-
-                int32_to_fract_array(attribs_begin, attribs_start, num_planes);
-
-                gen_steps(attribs_begin, attribs_end, attribs_steps,
-                          num_planes, column_delta);
+                prepare_for_interpolation(left_attribs, rght_attribs,
+                                          attribs, column_delta,
+                                          n);
 
                 if (left_boundary < 0) {
                     start_column = 0;
 
                     span_length += left_boundary;
 
-                    multi_step_up(attribs_start, attribs_steps, num_planes,
-                                  -left_boundary, column_delta);
+                    multi_step_up(attribs, -left_boundary, n);
                 }
 
                 if (rght_boundary >= fbi->width) {
@@ -101,24 +87,17 @@ draw_partial_triangle(
                 }
 
                 draw_span(
-                    ((bi->start_scanline + row) * fbi->width) + start_column,
-                    span_length, attribs_start, attribs_steps, column_delta,
-                    fbi);
+                    (bi->start_scanline + row) * fbi->width + start_column,
+                    span_length, attribs, fbi);
 
                 if (row_delta > 0) {
-                    step_up(left_attribs, left_attribs_steps, num_planes,
-                            left_div);
-                    step_up(rght_attribs, rght_attribs_steps, num_planes,
-                            rght_div);
+                    step_up(left_attribs, n);
+                    step_up(rght_attribs, n);
                 }
-                row++;
-                free(attribs_steps);
-                free(attribs_end);
-                free(attribs_begin);
-                free(attribs_start);
             }
         }
     }
+    free(attribs);
     free(rght_attribs);
     free(left_attribs);
 }
@@ -127,31 +106,19 @@ draw_partial_triangle(
 
 static void
 draw_degenerate_horizontal(Xy                 const xy,
-                           fract *            const attribs_left,
-                           fract *            const attribs_mid,
-                           const fract *      const top2mid_steps,
-                           const fract *      const top2bot_steps,
-                           const fract *      const mid2bot_steps,
-                           int32_t            const top2mid_delta,
-                           int32_t            const top2bot_delta,
-                           int32_t            const mid2bot_delta,
+                           varying *          const top2mid,
+                           varying *          const top2bot,
+                           varying *          const mid2bot,
                            framebuffer_info * const fbi) {
 
-    uint8_t const num_planes = fbi->num_attribs + 1;
-
-    fract * attribs_left_bkup;
-
-    MALLOCARRAY_NOFAIL(attribs_left_bkup, num_planes);
-
-    memcpy(attribs_left_bkup, attribs_left, num_planes * sizeof(fract));
+    uint8_t const n = fbi->num_attribs + 2;
 
     {
         int16_t const y = xy._[0][1];
 
         int16_t x[3];
         int16_t x_start[3];
-        fract * attribs[3];
-        const fract * steps[3];
+        varying * attribs[3];
         int32_t span_length[3];
         unsigned int i;
 
@@ -163,28 +130,21 @@ draw_degenerate_horizontal(Xy                 const xy,
         x_start[1] = x[0];
         x_start[2] = x[1];
 
-        attribs[0] = attribs_left;
-        attribs[1] = attribs_left_bkup;
-        attribs[2] = attribs_mid;
-
-        steps[0] = top2bot_steps;
-        steps[1] = top2mid_steps;
-        steps[2] = mid2bot_steps;
+        attribs[0] = top2bot;
+        attribs[1] = top2mid;
+        attribs[2] = mid2bot;
 
         span_length[0] = x[2] - x[0];
         span_length[1] = x[1] - x[0];
         span_length[2] = x[2] - x[1];
 
         for (i = 0; i < 3; i++) {
-            int32_t const column_delta = span_length[i];
-
             if (x_start[i] >= fbi->width || x_start[i] + span_length[i] < 0) {
                 continue;
             }
 
             if (x_start[i] < 0) {
-                multi_step_up(attribs[i], steps[i], num_planes, -x_start[i],
-                              column_delta);
+                multi_step_up(attribs[i], -x_start[i], n);
 
                 span_length[i] += x_start[i];
 
@@ -197,11 +157,10 @@ draw_degenerate_horizontal(Xy                 const xy,
                 span_length[i]++;
             }
 
-            draw_span((y * fbi->width) + x_start[i], span_length[i],
-                      attribs[i], steps[i], column_delta, fbi);
+            draw_span(y * fbi->width + x_start[i], span_length[i],
+                      attribs[i], fbi);
         }
     }
-    free(attribs_left_bkup);
 }
 
 
@@ -212,23 +171,28 @@ draw_triangle(Xy                 const xy_input,
               boundary_info *    const bi,
               framebuffer_info * const fbi) {
 
-    uint8_t const num_planes = fbi->num_attribs + 1;
+    uint8_t const z = fbi->num_attribs;
+    uint8_t const w = z + 1;
+    uint8_t const n = w + 1;
 
     Xy xy;
-    int32_t * attribs[3];
+    varying * attribs[3];
     unsigned int i;
     uint8_t index_array[3];
     int32_t y_array[3];
     int32_t x_array[3];
 
-    MALLOCARRAY_NOFAIL(attribs[0], num_planes);
-    MALLOCARRAY_NOFAIL(attribs[1], num_planes);
-    MALLOCARRAY_NOFAIL(attribs[2], num_planes);
+    MALLOCARRAY_NOFAIL(attribs[0], n);
+    MALLOCARRAY_NOFAIL(attribs[1], n);
+    MALLOCARRAY_NOFAIL(attribs[2], n);
 
     xy = xy_input;
 
     for (i = 0; i < 3; i++) {
-        memcpy(attribs[i], attribs_input._[i], num_planes * sizeof(int32_t));
+        int32_to_varying_array(attribs_input._[i], attribs[i], n);
+	attribs[i][z] = compute_varying_z(attribs_input._[i][z]);
+	attribs[i][w] = inverse_varying(attribs[i][w]);
+        multiply_varying_array_by_varying(attribs[i], attribs[i][w], z);
     }
 
     /* Argument preparations for sort3: */
@@ -259,122 +223,65 @@ draw_triangle(Xy                 const xy_input,
             gen_triangle_boundaries(xy_sorted, bi, fbi->width, fbi->height);
 
         if (bi->start_scanline == -1) {
-            /* Triangle is completely out of the bounds of the framebuffer. */
+            /* Triangle is completely out of the bounds of the frame buffer. */
         } else {
             bool const no_upper_part =
                 (xy_sorted._[1][1] == xy_sorted._[0][1]);
 
             bool const horizontal =
                 (xy._[0][1] == xy._[1][1] && xy._[1][1] == xy._[2][1]);
-                /* We are dealing with a degenerate horizontal triangle */
+                /* Tells whether we are dealing with a degenerate
+                 * horizontal triangle */
 
-            uint8_t t = ~horizontal & 1;
+            uint8_t const t = horizontal ^ 1;
 
             int32_t top2mid_delta = xy._[mid][t] - xy._[top][t];
             int32_t top2bot_delta = xy._[bot][t] - xy._[top][t];
             int32_t mid2bot_delta = xy._[bot][t] - xy._[mid][t];
 
-            fract * top2mid_steps;
-            fract * top2bot_steps;
-            fract * mid2bot_steps;
+            varying * top2mid;
+            varying * top2bot;
+            varying * mid2bot;
 
-            fract * upper_left_attribs_steps;
-            fract * lower_left_attribs_steps;
-            fract * upper_rght_attribs_steps;
-            fract * lower_rght_attribs_steps;
+            varying * upper_left_attribs;
+            varying * lower_left_attribs;
+            varying * upper_rght_attribs;
+            varying * lower_rght_attribs;
 
-            int32_t upper_left_delta;
-            int32_t lower_left_delta;
-            int32_t upper_rght_delta;
-            int32_t lower_rght_delta;
+            MALLOCARRAY_NOFAIL(top2mid, n);
+            MALLOCARRAY_NOFAIL(top2bot, n);
+            MALLOCARRAY_NOFAIL(mid2bot, n);
 
-            fract * left_attribs;
-            fract * rght_attribs;
-
-            bool degenerate_horizontal;
-
-            MALLOCARRAY_NOFAIL(top2mid_steps, num_planes);
-            MALLOCARRAY_NOFAIL(top2bot_steps, num_planes);
-            MALLOCARRAY_NOFAIL(mid2bot_steps, num_planes);
-            MALLOCARRAY_NOFAIL(left_attribs, num_planes);
-            MALLOCARRAY_NOFAIL(rght_attribs, num_planes);
-
-            if (!horizontal) {
-                top2mid_delta += !no_upper_part;
-                top2bot_delta += 1;
-                mid2bot_delta += no_upper_part;
-            }
-
-            gen_steps(attribs[top], attribs[mid], top2mid_steps, num_planes,
-                      top2mid_delta);
-            gen_steps(attribs[top], attribs[bot], top2bot_steps, num_planes,
-                      top2bot_delta);
-            gen_steps(attribs[mid], attribs[bot], mid2bot_steps, num_planes,
-                      mid2bot_delta);
-
-            int32_to_fract_array(attribs[top], left_attribs, num_planes);
-            int32_to_fract_array(attribs[top], rght_attribs, num_planes);
+            prepare_for_interpolation(attribs[top], attribs[mid], top2mid, top2mid_delta, n);
+            prepare_for_interpolation(attribs[top], attribs[bot], top2bot, top2bot_delta, n);
+            prepare_for_interpolation(attribs[mid], attribs[bot], mid2bot, mid2bot_delta, n);
 
             if (mid_is_to_the_left) {
-                upper_left_attribs_steps = top2mid_steps;
-                lower_left_attribs_steps = mid2bot_steps;
-                upper_rght_attribs_steps = top2bot_steps;
-                lower_rght_attribs_steps = upper_rght_attribs_steps;
-
-                upper_left_delta = top2mid_delta;
-                lower_left_delta = mid2bot_delta;
-                upper_rght_delta = top2bot_delta;
-                lower_rght_delta = upper_rght_delta;
+                upper_left_attribs = top2mid;
+                lower_left_attribs = mid2bot;
+                upper_rght_attribs = top2bot;
+                lower_rght_attribs = upper_rght_attribs;
             } else {
-                upper_rght_attribs_steps = top2mid_steps;
-                lower_rght_attribs_steps = mid2bot_steps;
-                upper_left_attribs_steps = top2bot_steps;
-                lower_left_attribs_steps = upper_left_attribs_steps;
-
-                upper_rght_delta = top2mid_delta;
-                lower_rght_delta = mid2bot_delta;
-                upper_left_delta = top2bot_delta;
-                lower_left_delta = upper_left_delta;
+                upper_rght_attribs = top2mid;
+                lower_rght_attribs = mid2bot;
+                upper_left_attribs = top2bot;
+                lower_left_attribs = upper_left_attribs;
             }
 
-            if (no_upper_part) {
-                int32_to_fract_array(attribs[mid], rght_attribs, num_planes);
-
-                if (horizontal) {
-                    degenerate_horizontal = true;
-                } else {
-                    degenerate_horizontal = false;
-
-                    step_up(left_attribs, lower_left_attribs_steps, num_planes,
-                            lower_left_delta);
-                    step_up(rght_attribs, lower_rght_attribs_steps, num_planes,
-                            lower_rght_delta);
-                }
-            } else {
+            if (!(horizontal || no_upper_part)) {
                 int32_t delta;
 
-                degenerate_horizontal = false;
-
-                step_up(left_attribs, upper_left_attribs_steps, num_planes,
-                        upper_left_delta);
-                step_up(rght_attribs, upper_rght_attribs_steps, num_planes,
-                        upper_rght_delta);
-
                 if (bi->num_upper_rows > 0) {
-
                     if (bi->start_scanline > xy._[top][1]) {
                         delta = bi->start_scanline - xy._[top][1];
 
-                        multi_step_up(left_attribs, upper_left_attribs_steps,
-                                      num_planes, delta, upper_left_delta);
-                        multi_step_up(rght_attribs, upper_rght_attribs_steps,
-                                      num_planes, delta, upper_rght_delta);
+                        multi_step_up(upper_left_attribs, delta, n);
+                        multi_step_up(upper_rght_attribs, delta, n);
                     }
 
                     draw_partial_triangle(
-                        left_attribs, upper_left_attribs_steps,
-                        rght_attribs, upper_rght_attribs_steps,
-                        upper_left_delta, upper_rght_delta,
+                        upper_left_attribs,
+                        upper_rght_attribs,
                         true,
                         bi,
                         fbi
@@ -385,40 +292,33 @@ draw_triangle(Xy                 const xy_input,
                     delta = top2mid_delta;
                 }
 
-                multi_step_up(left_attribs, upper_left_attribs_steps,
-                              num_planes, delta, upper_left_delta);
-                multi_step_up(rght_attribs, upper_rght_attribs_steps,
-                              num_planes, delta, upper_rght_delta);
+                multi_step_up(upper_left_attribs, delta, n);
+                multi_step_up(upper_rght_attribs, delta, n);
             }
-            if (degenerate_horizontal) {
+
+            if (horizontal) {
                 draw_degenerate_horizontal(
                     xy_sorted,
-                    left_attribs, rght_attribs,
-                    top2mid_steps, top2bot_steps, mid2bot_steps,
-                    top2mid_delta, top2bot_delta, mid2bot_delta,
+                    top2mid, top2bot, mid2bot,
                     fbi
                     );
             } else {
                 if (bi->start_scanline > xy._[mid][1]) {
                     int32_t const delta = bi->start_scanline - xy._[mid][1];
 
-                    multi_step_up(left_attribs, lower_left_attribs_steps,
-                                  num_planes, delta, lower_left_delta);
-                    multi_step_up(rght_attribs, lower_rght_attribs_steps,
-                                  num_planes, delta, lower_rght_delta);
+                    multi_step_up(lower_left_attribs, delta, n);
+                    multi_step_up(lower_rght_attribs, delta, n);
                 }
 
                 draw_partial_triangle(
-                    left_attribs, lower_left_attribs_steps,
-                    rght_attribs, lower_rght_attribs_steps,
-                    lower_left_delta, lower_rght_delta,
+                    lower_left_attribs,
+                    lower_rght_attribs,
                     false,
                     bi,
                     fbi
                     );
             }
-            free(rght_attribs); free(left_attribs);
-            free(mid2bot_steps); free(top2bot_steps); free(top2mid_steps);
+            free(mid2bot); free(top2bot); free(top2mid);
         }
     }
     free(attribs[2]); free(attribs[1]); free(attribs[0]);
diff --git a/generator/pamtris/triangle.h b/generator/pamtris/triangle.h
index 79178ad0..e043e95c 100644
--- a/generator/pamtris/triangle.h
+++ b/generator/pamtris/triangle.h
@@ -13,7 +13,7 @@ typedef struct {
 } Xy;
 
 typedef struct {
-    int32_t _[3][MAX_NUM_ATTRIBS + 1];
+    int32_t _[3][MAX_NUM_ATTRIBS + 2];
 } Attribs;
 
 void
diff --git a/generator/pamtris/utils.c b/generator/pamtris/utils.c
index 09c9b4d0..a6b6e7d4 100644
--- a/generator/pamtris/utils.c
+++ b/generator/pamtris/utils.c
@@ -6,171 +6,181 @@
 
 #include <stdlib.h>
 #include <stdint.h>
+#include <math.h>
 
-#include "fract.h"
+#include "limits_pamtris.h"
+#include "varying.h"
 
 #include "utils.h"
 
 
 
 void
-step_up(fract *       const vars,
-        const fract * const steps,
-        uint8_t       const element_ct,
-        int32_t       const divisor) {
-/*----------------------------------------------------------------------------
-  Apply interpolation steps steps[] to a collection of fract variables vars[]
-  once.  I.e. add each steps[i] to vars[i].
+prepare_for_interpolation(const varying * const begin,
+                          const varying * const end,
+                          varying *       const out,
+                          int32_t               num_steps,
+                          uint8_t         const elements) {
 
-  'element_ct' is the number of elements in 'vars' and 'steps'.
+    double inverse_num_steps;
+    unsigned int i;
 
-  'divisor' is the divisor used to interpret the fractions.
+    if (num_steps < 1) {
+        num_steps = 1;
+    }
 
-  It *is* safe to pass a 0 divisor to this function.
------------------------------------------------------------------------------*/
-    unsigned int i;
+    inverse_num_steps = 1.0 / num_steps;
 
-    for (i = 0; i < element_ct; ++i) {
-        /* To add the fraction steps[i] to the fraction vars[i]: add the
-           quotient of step steps[i] to the quotient of variable vars[i] and
-           the remainder of the step to the remainder of the variable. If this
-           makes the agumented remainder equal to or larger than the divisor,
-           increment the quotient of the variable if the step is positive or
-           decrement it if the step is negative, and subtract the divisor from
-           the remainder of the variable (in either case).
-        */
+    for (i = 0; i < elements; i++) {
+        out[i].v = begin[i].v;
+        out[i].s = (end[i].v - begin[i].v) * inverse_num_steps;
+    }
+}
 
-        vars[i].q += steps[i].q;
-        vars[i].r += steps[i].r;
 
-        {
-            uint32_t const negative_mask = -steps[i].negative_flag;
-                /* (-1 if the step is negative; 1 otherwise) */
 
-            uint32_t const overdiv_mask =
-                -(((uint32_t)~(vars[i].r - divisor)) >> 31);
-                /*  = ~0 if var->r >= div; 0 otherwise. */
+varying
+compute_varying_z(int32_t const input_z) {
 
-            vars[i].q += (negative_mask | 1) & overdiv_mask;
-            vars[i].r -= divisor & overdiv_mask;
-        }
-    }
+    varying retval;
+
+    retval.v = 1.0 / (1 + input_z - MIN_COORD);
+    retval.s = 0.0;
+
+    return retval;
 }
 
 
 
 void
-multi_step_up(fract *       const vars,
-              const fract * const steps,
-              uint8_t       const elements,
-              int32_t       const times,
-              int32_t       const div) {
-/*----------------------------------------------------------------------------
-  Similar to step_up, but apply the interpolation step an arbitrary number
-  of times, instead of just once.
+multiply_varying_array_by_varying(varying * const vars,
+                                  varying   const multiplier,
+                                  uint8_t   const elements) {
 
-  It *is* also safe to pass a 0 divisor to this function.
------------------------------------------------------------------------------*/
     unsigned int i;
 
     for (i = 0; i < elements; i++) {
-        uint32_t const negative_mask = -steps[i].negative_flag;
+        vars[i].v *= multiplier.v;
+	vars[i].s  = 0.0;
+    }
+}
 
-        vars[i].q += times * steps[i].q;
-        vars[i].r += times * steps[i].r;
 
-        if(vars[i].r >= div && div != 0) {
-            int32_t const r_q = vars[i].r / div;
-            int32_t const r_r = vars[i].r % div;
+void
+divide_varying_array_by_varying(varying * const vars,
+                                varying   const divisor,
+                                uint8_t   const elements) {
 
-            vars[i].q += (-r_q & negative_mask) | (r_q & ~negative_mask);
-                /* = -r_q if the step is negative; r_q, otherwise. */
-            vars[i].r = r_r;
-        }
+    double const inverse_divisor = 1.0 / divisor.v;
+
+    unsigned int i;
+
+    for (i = 0; i < elements; i++) {
+        vars[i].v *= inverse_divisor;
+	vars[i].s  = 0.0;
     }
 }
 
 
 
-void
-gen_steps(const int32_t * const begin,
-          const int32_t * const end,
-          fract         * const out,
-          uint8_t         const elements,
-          int32_t         const div) {
-/*----------------------------------------------------------------------------
-  Generate the interpolation steps for a collection of initial and final
-  values. "begin" points to an array of initial values, "end" points to the
-  array of corresponding final values; each interpolation step is stored in
-  the appropriate position in the array pointed by "out"; "elements" indicates
-  the number of elements in each of the previously mentioned arrays and
-  "divisor" is the common value by which we want to divide the difference
-  between each element in the array pointed to by "end" and the corresponding
-  element in the array pointed to by "begin".  After an execution of this
-  function, for each out[i], with 0 <= i < elements, the following will hold:
-
-    1. If divisor > 1:
-      out[i].q = (end[i] - begin[i]) / divisor
-      out[i].r = abs((end[i] - begin[i]) % divisor)
-
-    2. If divisor == 1 || divisor == 0:
-      out[i].q = end[i] - begin[i]
-      out[i].r = 0
------------------------------------------------------------------------------*/
-    if (div > 1) {
-        unsigned int i;
+varying
+inverse_varying(varying const var) {
 
-        for (i = 0; i < elements; i++) {
-            int32_t const delta = end[i] - begin[i];
+    varying retval;
 
-            out[i].q = delta / div;
-            out[i].r = abs(delta % div);
-            out[i].negative_flag = ((uint32_t)delta) >> 31;
-        }
-    } else {
-        unsigned int i;
+    retval.v = 1.0 / var.v;
+    retval.s = 0.0;
 
-        for (i = 0; i < elements; i++) {
-            int32_t const delta = end[i] - begin[i];
+    return retval;
+}
 
-            out[i].q = delta;
-            out[i].r = 0;
-            out[i].negative_flag = ((uint32_t)delta) >> 31;
-        }
+
+
+varying
+multiply_varyings(varying const a,
+                  varying const b) {
+
+    varying retval;
+
+    retval.v = a.v * b.v;
+    retval.s = 0.0;
+
+    return retval;
+}
+
+
+
+void
+step_up(varying * const vars,
+        uint8_t   const elements) {
+
+    unsigned int i;
+
+    for (i = 0; i < elements; i++) {
+        vars[i].v += vars[i].s;
     }
 }
 
 
 
 void
-fract_to_int32_array(const fract * const in,
-                     int32_t *     const out,
-                     uint8_t       const elements) {
+multi_step_up(varying * const vars,
+              int32_t   const times,
+              uint8_t   const elements) {
 
     unsigned int i;
 
     for (i = 0; i < elements; i++) {
-        out[i] = in[i].q;
+        vars[i].v += times * vars[i].s;
     }
 }
 
 
 
 void
-int32_to_fract_array(const int32_t * const in,
-                     fract *         const out,
-                     uint8_t         const elements) {
+int32_to_varying_array(const int32_t * const in,
+                       varying *       const out,
+                       uint8_t         const elements) {
 
     unsigned int i;
 
     for (i = 0; i < elements; i++) {
-        out[i].q = in[i];
-        out[i].r = 0;
+        out[i].v = in[i];
+        out[i].s = 0.0;
     }
 }
 
 
 
+/* static int64_t
+abs64(int64_t x)
+{
+
+    int64_t const nm = ~geq_mask64(x, 0);
+
+    return (-x & nm) | (x & ~nm);
+} */
+
+
+
+int32_t
+round_varying(varying const var) {
+
+    return round(var.v);
+}
+
+
+
+int64_t
+geq_mask64(int64_t a, int64_t b) {
+
+    uint64_t const diff = a - b;
+
+    return -((~diff) >> 63);
+}
+
+
+
 static void
 swap(uint8_t * const a,
      uint8_t * const b) {
diff --git a/generator/pamtris/utils.h b/generator/pamtris/utils.h
index 3b7cfbe4..bd9dcdbe 100644
--- a/generator/pamtris/utils.h
+++ b/generator/pamtris/utils.h
@@ -1,37 +1,54 @@
 #ifndef UTIL_H_INCLUDED
 #define UTIL_H_INCLUDED
 
-#include "fract.h"
+#include "varying.h"
 
 void
-gen_steps(const int32_t * const begin,
-          const int32_t * const end,
-          fract *         const out,
-          uint8_t         const elements,
-          int32_t         const divisor);
+prepare_for_interpolation(const varying * const begin,
+                          const varying * const end,
+                          varying *       const out,
+                          int32_t               num_steps,
+                          uint8_t         const elements);
+
+varying
+compute_varying_z(int32_t const input_z);
 
 void
-step_up(fract *       const vars,
-        const fract * const steps,
-        uint8_t       const elements,
-        int32_t       const divisor);
+multiply_varying_array_by_varying(varying * const vars,
+                                  varying   const divisor,
+                                  uint8_t   const elements);
 
 void
-multi_step_up(fract *       const vars,
-              const fract * const steps,
-              uint8_t       const elements,
-              int32_t       const times,
-              int32_t       const divisor);
+divide_varying_array_by_varying(varying * const vars,
+                                varying   const divisor,
+                                uint8_t   const elements);
+
+varying
+inverse_varying(varying const var);
+
+varying
+multiply_varyings(varying const a,
+                  varying const b);
 
 void
-fract_to_int32_array(const fract * const in,
-                     int32_t     * const out,
-                     uint8_t       const elements);
+step_up(varying * const vars,
+       uint8_t    const elements);
 
 void
-int32_to_fract_array(const int32_t * const in,
-                     fract *         const out,
-                     uint8_t         const elements);
+multi_step_up(varying * const vars,
+             int32_t    const times,
+             uint8_t    const elements);
+
+void
+int32_to_varying_array(const int32_t * const in,
+                       varying *       const out,
+                       uint8_t         const elements);
+
+int32_t
+round_varying(varying const var);
+
+int64_t
+geq_mask64(int64_t a, int64_t b);
 
 void
 sort3(uint8_t *       const index_array,
diff --git a/generator/pamtris/varying.h b/generator/pamtris/varying.h
new file mode 100644
index 00000000..6605f02d
--- /dev/null
+++ b/generator/pamtris/varying.h
@@ -0,0 +1,12 @@
+#ifndef VARYING_H_INCLUDED
+#define VARYING_H_INCLUDED
+
+#include <stdint.h>
+
+
+typedef struct {
+    double v; /* Value */
+    double s; /* Step */
+} varying;
+
+#endif