[tessellator] Direct comparison of result in edges_compare_x_for_y.

We need to compare the x-coordinate of a line at a for a particular y,
without loss of precision.

The x-coordinate along an edge for a given y is:
  X = A_x + (Y - A_y) * A_dx / A_dy

So the inequality we wish to test is:
  A_x + (Y - A_y) * A_dx / A_dy -?- B_x + (Y - B_y) * B_dx / B_dy,
where -?- is our inequality operator.

By construction we know that A_dy and B_dy (and (Y - A_y), (Y - B_y)) are
all positive, so we can rearrange it thus without causing a sign
change:
  A_dy * B_dy * (A_x - B_x) -?- (Y - B_y) * B_dx * A_dy
                                - (Y - A_y) * A_dx * B_dy

Given the assumption that all the deltas fit within 32 bits, we can compute
this comparison directly using 128 bit arithmetic.

1 files changed, 55 insertions, 42 deletions

diff --git a/src/cairo-bentley-ottmann.c b/src/cairo-bentley-ottmann.c
index a78745d1b..ddcd9f0a1 100644
--- a/src/cairo-bentley-ottmann.c
+++ b/src/cairo-bentley-ottmann.c
@@ -219,42 +219,66 @@ _slope_compare (cairo_bo_edge_t *a,
     }
 }
 
-static cairo_quorem64_t
-edge_x_for_y (cairo_bo_edge_t *edge,
-	      int32_t y)
+/*
+ * We need to compare the x-coordinate of a line for a particular y,
+ * without loss of precision.
+ *
+ * The x-coordinate along an edge for a given y is:
+ *   X = A_x + (Y - A_y) * A_dx / A_dy
+ *
+ * So the inequality we wish to test is:
+ *   A_x + (Y - A_y) * A_dx / A_dy -?- B_x + (Y - B_y) * B_dx / B_dy,
+ * where -?- is our inequality operator.
+ *
+ * By construction, we know that A_dy and B_dy (and (Y - A_y), (Y - B_y)) are
+ * all positive, so we can rearrange it thus without causing a sign change:
+ *   A_dy * B_dy * (A_x - B_x) -?- (Y - B_y) * B_dx * A_dy
+ *                                 - (Y - A_y) * A_dx * B_dy
+ *
+ * Given the assumption that all the deltas fit within 32 bits, we can compute
+ * this comparison directly using 128 bit arithmetic.
+ *
+ * (And put the burden of the work on developing fast 128 bit ops, which are
+ * required throughout the tessellator.)
+ *
+ * See the similar discussion for _slope_compare().
+ */
+static int
+edges_compare_x_for_y (const cairo_bo_edge_t *a,
+		       const cairo_bo_edge_t *b,
+		       int32_t y)
 {
     /* XXX: We're assuming here that dx and dy will still fit in 32
      * bits. That's not true in general as there could be overflow. We
      * should prevent that before the tessellation algorithm
      * begins.
      */
-    int32_t dx = edge->bottom.x - edge->top.x;
-    int32_t dy = edge->bottom.y - edge->top.y;
-    int64_t numerator;
-    cairo_quorem64_t quorem;
-
-    if (edge->middle.y == y) {
-       quorem.quo = edge->middle.x;
-       quorem.rem = 0;
-       return quorem;
-    }
-    if (edge->bottom.y == y) {
-       quorem.quo = edge->bottom.x;
-       quorem.rem = 0;
-       return quorem;
-    }
-    if (dy == 0) {
-	quorem.quo = _cairo_int32_to_int64 (edge->top.x);
-	quorem.rem = 0;
-	return quorem;
-    }
+    int32_t adx, ady;
+    int32_t bdx, bdy;
+    cairo_int128_t L, R;
+
+    adx = a->bottom.x - a->top.x;
+    ady = a->bottom.y - a->top.y;
 
-    /* edge->top.x + (y - edge->top.y) * dx / dy */
-    numerator = _cairo_int32x32_64_mul ((y - edge->top.y), dx);
-    quorem = _cairo_int64_divrem (numerator, dy);
-    quorem.quo += edge->top.x;
+    bdx = b->bottom.x - b->top.x;
+    bdy = b->bottom.y - b->top.y;
 
-    return quorem;
+    L = _cairo_int64x32_128_mul (_cairo_int32x32_64_mul (ady, bdy),
+				 a->top.x - b->top.x);
+
+    R = _cairo_int128_sub (_cairo_int64x32_128_mul (_cairo_int32x32_64_mul (bdx,
+									    ady),
+						    y - b->top.y),
+			   _cairo_int64x32_128_mul (_cairo_int32x32_64_mul (adx,
+									    bdy),
+						    y - a->top.y));
+
+    /* return _cairo_int128_cmp (L, R); */
+    if (_cairo_int128_lt (L, R))
+	return -1;
+    if (_cairo_int128_gt (L, R))
+	return 1;
+    return 0;
 }
 
 static int
@@ -262,8 +286,6 @@ _cairo_bo_sweep_line_compare_edges (cairo_bo_sweep_line_t	*sweep_line,
 				    cairo_bo_edge_t		*a,
 				    cairo_bo_edge_t		*b)
 {
-    cairo_quorem64_t ax;
-    cairo_quorem64_t bx;
     int cmp;
 
     if (a == b)
@@ -292,18 +314,9 @@ _cairo_bo_sweep_line_compare_edges (cairo_bo_sweep_line_t	*sweep_line,
            if (amin > bmax) return +1;
     }
 
-    ax = edge_x_for_y (a, sweep_line->current_y);
-    bx = edge_x_for_y (b, sweep_line->current_y);
-    if (ax.quo > bx.quo)
-	return 1;
-    else if (ax.quo < bx.quo)
-	return -1;
-
-    /* Quotients are identical, test remainder. */
-    if (ax.rem > bx.rem)
-	return 1;
-    else if (ax.rem < bx.rem)
-	return -1;
+    cmp = edges_compare_x_for_y (a, b, sweep_line->current_y);
+    if (cmp)
+	return cmp;
 
     /* The two edges intersect exactly at y, so fall back on slope
      * comparison. We know that this compare_edges function will be