1  #import "@preview/ctheorems:1.1.3": *
  2  #import "@preview/cetz:0.4.0"
  3  #import "@preview/lovelace:0.3.0": *
  4  #set par(justify: true)
  5  #set math.equation(numbering: "(1)")
  6  
  7  #show: thmrules
  8  #let lemma = thmbox("lemma", "Lemma")
  9  #let def = thmbox("lemma", "Definition")
 10  #let proof = thmproof("proof", "Proof")
 11  #let remark = thmbox("lemma", "Remark")
 12  
 13  #let inexact(term) = {
 14    block(inset: 16pt,
 15      block(
 16        inset: 16pt,
 17        stroke: 1pt + gray,
 18        radius: 12pt,
 19        text(
 20          size: 10pt,
 21          [Note for inexact arithmetic: #term]
 22        )
 23      )
 24    )
 25  }
 26  
 27  // TODO: figure out how to get rid of the ":" in the caption
 28  #let invariant(term) = {
 29  figure(
 30    block(inset: 16pt,
 31        text(
 32          size: 10pt,
 33          [#term]
 34      )
 35    ),
 36    kind: "invariant",
 37    supplement: "invariant",
 38    caption: ""
 39  )
 40  }
 41  
 42  #let fl = text([fl])
 43  
 44  We'll be talking about sweep line algorithms, where the sweep line is horizontal and increasing in $y$.
 45  Therefore, every line segment "starts" at the coordinate with smaller $y$ and "ends" at the coordinate
 46  with larger $y$ (we'll assume for now that there are no horizontal segments). We'll parametrize each
 47  line segment as a function of $y$. So if $alpha: [y_0 (alpha), y_1 (alpha)] arrow bb(R)$ is a line segment
 48  then $alpha(y)$ is the $x$ coordinate at $y$-coordinate $y in [y_0 (alpha), y_1 (alpha)]$.
 49  We write $theta_alpha$ for the angle that $alpha$ makes with the positive horizontal axis.
 50  Let's have a picture. (In the discussion, it won't matter whether positive $y$ points up or down, but in the
 51  pictures we'll adopt the graphics convention of having positive $y$ point down.)
 52  
 53  #figure(
 54  cetz.canvas({
 55    import cetz.draw: *
 56  
 57    line((1, 3), (0, 0), name: "a")
 58    line((-4, 3), (4, 3), stroke: (dash: "dotted"))
 59    line((-4, 0), (4, 0), stroke: (dash: "dotted"))
 60    content((-4, 3), $y_0(alpha)$, anchor: "east")
 61    content((-4, 0), $y_1(alpha)$, anchor: "east")
 62    content((0.6, 1.5), $alpha$, anchor: "west")
 63  
 64    cetz.angle.angle("a.start", "a.end", (4, 3), label: $theta_alpha$, label-radius: 0.8)
 65  }),
 66  
 67  caption: "A segment and its angle"
 68  )
 69  
 70  We'll be dealing with inexact arithmetic, so let's define some "error bars" on our line segments.
 71  For an error parameter $epsilon > 0$, offsetting from $alpha$ by $plus.minus epsilon$ in the perpendicular-to-$alpha$ direction
 72  is the same as offsetting by $alpha plus.minus epsilon / (|sin theta_alpha|)$ in the horizontal direction.
 73  Roughly speaking, the "error bars" on $alpha$ amount to adding this horizontal error. But we'll be slightly
 74  more accurate around the corners, by truncating these error bars to the horizontal extents of $alpha$. Precisely, we define
 75  
 76  $
 77  alpha_(+,epsilon)(y) = min(alpha(y) + epsilon / (|sin theta_alpha|), max(alpha(y_0), alpha(y_1)) + epsilon) \
 78  alpha_(-,epsilon)(y) = max(alpha(y) - epsilon / (|sin theta_alpha|), min(alpha(y_0), alpha(y_1)) - epsilon) \
 79  $
 80  
 81  In pictures, the gray shaded region is the region between $alpha_(-,epsilon)$ and $alpha_(+,epsilon)$:
 82  The point of the scaling by $|sin theta_alpha|$ is to make this an approximation of an $epsilon$-neighborhood (in the
 83  perpendicular direction) of the line segment. The truncation near the corners ensures that if $x$ is between
 84  $alpha_(-,epsilon)(y)$ and $alpha_(+,epsilon)(y)$ then it is within $sqrt(2) epsilon$ of $alpha$.
 85  
 86  #figure(
 87  cetz.canvas({
 88    import cetz.draw: *
 89  
 90    line((0.5, 3), (-0.3, 0.6), (-0.3, 0), (0.5, 0), (1.3, 2.4), (1.3, 3), close: true, fill: gray, stroke: 0pt)
 91    line((0.5, 3), (-0.3, 0.6), (-0.3, 0), stroke: ( dash: "dashed" ))
 92    line((0.5, 0), (1.3, 2.4), (1.3, 3), stroke: ( dash: "dashed" ))
 93  
 94    line((1, 3), (0, 0), name: "a")
 95    line((-4, 3), (4, 3), stroke: (dash: "dotted"))
 96    line((-4, 0), (4, 0), stroke: (dash: "dotted"))
 97    content((-4, 3), $y_0(alpha)$, anchor: "east")
 98    content((-4, 0), $y_1(alpha)$, anchor: "east")
 99    content((0.6, 1.5), $alpha$, anchor: "west")
100  
101    content((0.8, 0.5), $alpha_(+,epsilon)$, anchor: "west")
102    content((0.2, 2.4), $alpha_(-,epsilon)$, anchor: "east")
103  }),
104  caption: "A segment and its error bars"
105  )<fig-error-bars>
106  
107  
108  Define a relation $prec_(y,epsilon)$ on line segments whose domain contains $y$, where
109  $alpha prec_(y,epsilon) beta$ if $alpha_(+,epsilon)(y) < beta_(-,epsilon)(y)$.
110  Intuitively, $alpha prec_(y,epsilon) beta$ if $alpha$ is definitely to the left of $beta$
111  at height $y$: $alpha$ is far enough to the left that their error bars don't overlap.
112  Clearly, for a given $y$ and $epsilon$ there are three mutually exclusive possibilities: either
113  $alpha prec_(y,epsilon) beta$ or $beta prec_(y,epsilon) alpha$ or neither of the two holds. We'll denote
114  this third possibility by $alpha approx_(y,epsilon) beta$, and we write
115  $alpha prec.tilde_(y,epsilon) beta$ for "$alpha prec_(y,epsilon) beta$ or $alpha approx_(y,epsilon) beta$."
116  
117  Here are a few basic properties of our definitions:
118  #lemma[
119  1. For any $y$ and any $epsilon > 0$, $prec_(y,epsilon)$ is transitive:
120    if $alpha prec_(y,epsilon) beta$
121    and $beta prec_(y,epsilon) gamma$ then $alpha prec_(y,epsilon) gamma$. (However, $prec.tilde_(y,epsilon)$ is not transitive.)
122  2. For any $y$ and any $epsilon > 0$,
123    if $alpha prec_(y,epsilon) beta$
124    and $beta prec.tilde_(y,epsilon) gamma$ then $alpha prec.tilde_(y,epsilon) gamma$.
125  3. For any $y$ and any $epsilon > 0$,
126    if $alpha prec.tilde_(y,epsilon) beta$
127    and $beta prec_(y,epsilon) gamma$ then $alpha prec.tilde_(y,epsilon) gamma$.
128  4. For any $y$, the relation $prec_(y,epsilon)$ is monotone in $epsilon$, in that if $alpha prec_(y,epsilon) beta$ then $alpha prec_(y,eta) beta$ for
129    any $eta in (0, epsilon)$.
130  ]<lem-basic-order-properties>
131  
132  Since $epsilon$ for us will usually be fixed, we will often drop it from the notation, and write $alpha_-$ and $alpha_+$ instead
133  of $alpha_(-,epsilon)$ and $alpha_(+,epsilon)$.
134  
135  #def[
136    Suppose $alpha$ and $beta$ are two segments whose domain contains $y$. We say that *$(alpha, beta)$
137    $epsilon$-cross by $y$* if $y$ belongs to both domains and $alpha succ_(y,epsilon) beta$.
138    We say that *$(alpha, beta)$ $epsilon$-cross* if they $epsilon$-cross by $min(y_1 (alpha), y_1 (beta))$.
139  
140    When $epsilon$ is zero, we leave it out: we say that $(alpha, beta)$ cross by $y$ if they $0$-cross by $y$.
141  ]
142  
143  Note that the definition of $epsilon$-crossing is not symmetric: $(alpha, beta)$ $epsilon$-crossing is
144  not the same as $(beta, alpha)$ $epsilon$-crossing. We will usually talk about $(alpha, beta)$ $epsilon$-crossing
145  in the context that $alpha$ starts off to the left of $beta$, and in this context "$(alpha, beta)$ $epsilon$-cross" means
146  that at some height before the end of $alpha$ and $beta$, $alpha$ has definitely crossed to the right of $beta$.
147  
148  == Partially ordered sweep-lines
149  
150  At our first pass, we won't try to detect intersections at all. Instead, we'll produce
151  a continuum of sweep-lines (constant except at a finite number of points) that *approximately*
152  track the horizontal order of the segments.
153  
154  #def[
155  The ordered collection $(alpha^1, ..., alpha^m)$ of line segments is #emph[$epsilon$-ordered at $y$]
156  if each $alpha^i$ has $y$ in its domain and $alpha^i prec.tilde_(y,epsilon) alpha^j$ for all $1 <= i < j <= m$.
157  ]
158  
159  Our algorithm will produce a family of sweep-lines that are $epsilon$-ordered at every $y$
160  (and also the sweep-lines will be #emph[complete] in the sense that the sweep-line at $y$ will contain all line segments whose
161  domain contains $y$). This seems weaker than finding all the intersections (for example, because if you
162  find all intersections you can use them to produce a completely ordered family of sweep-lines), but
163  in fact they're more-or-less equivalent: given weakly-ordered sweep-lines, you can perturb the lines so
164  that your weak order becomes the real order of the perturbed lines.
165  
166  #lemma[
167  If $(alpha^1, ..., alpha^m)$ is $epsilon$-ordered at $y$ then there exist $x^1 <= ... <= x^m$ such that
168  $alpha_(-,epsilon)^i (y) <= x^i <= alpha_(+,epsilon)^i (y)$ for all $i$.
169  ]<lem-horizontal-realization>
170  
171  #proof[
172  Define $x^i = max_(j <= i) alpha_(-,epsilon)^j (y)$.
173  Since we've included $j = i$, this clearly satisfies $x^i >= alpha_(-,epsilon)^j (y)$.
174  For the other inequality, the ordering condition implies that $alpha^j_(-,epsilon) (y) <= alpha^i_(+,epsilon) (y)$
175  for every $j <= i$. Therefore this inequality still holds for the maximum over these $j$. 
176  ]
177  
178  One consequence of our approximate approach is that we need to do a little extra bookkeeping to maintain
179  the continuity of the input paths: when one segment exits and its path-neighbor enters, we need to remember
180  that they are connected because the approximate sweep-line might not keep them together. We'll ignore this
181  bookkeeping for now;
182  the goal here is to get into detail about the sweep-line invariants, and prove that we can maintain them.
183  
184  == The sweep-line invariants
185  
186  We're going to have a sweep-line that depends on $y$. When we need to emphasize this, we'll use the
187  cumbersome but explicit notation
188  $(alpha_y^1, ..., alpha_y^(m_y))$.
189  In addition to the sweep-line, we'll maintain a queue of "events." The events are ordered by $y$-position
190  and they are similar to the classical Bentley-Ottmann events so we'll skimp on the details here. There
191  is an "enter" event, an "exit" event, and an "intersection" event.
192  
193  Our sweep-line will maintain two invariants:
194  + At every $y$, the sweep-line is $epsilon$-ordered at $y$. (We'll call this the "order" invariant.)
195  + For every $y$ and every $1 <= i < j <= m_y$, if $alpha_y^i$ and $alpha_y^j$ $epsilon$-cross
196    then the event queue contains an event between $i$ and $j$,
197    and at least one of these events occurs before the $epsilon$-crossing height, or at $y$.
198    (We'll call this the "crossing" invariant.)
199  
200  (When we say the event queue contains an event between $i$ and $j$, we mean that either there's
201  an exit event for some $alpha^k$ with $i < k < j$ or there's an intersection event for some $alpha^k$ and $alpha^ell$
202  with $i <= k < ell <= j$.)
203  
204  Hopefully the first invariant is already well-motivated, so let's discuss the second.
205  To naively ensure that we find
206  all intersections, we could adopt a stronger crossing invariant that requires all pairs of intersecting segments
207  in the sweep-line to immediately put an intersection event in the queue. This would be inefficient to maintain because
208  it would require testing all pairs of segments. The main Bentley-Ottmann observation is that it's enough to have intersection
209  events for all sweep-line-adjacent pairs, because any pair of lines will have to become sweep-line adjacent strictly
210  before they intersect. We can't rely only on sweep-line adjacency because of the $epsilon$ fudge, but our "crossing"
211  event essentially encodes the same property. Note that if $j = i+1$ and $y$ is just before the $epsilon$-crossing height and
212  there are no other segments nearby, then the mysterious $j'$ must be either $i$ or $j$ (because there is nothing in between)
213  and the mysterious queue event must be the crossing of $i$ and $j$ (it couldn't be an exit event, because we assumed they
214  cross and they must cross before they exit). In other cases, the crossing event ensures that even if we haven't
215  recorded the upcoming crossing of $alpha^i$ and $alpha^j$, something will happen in between them that will give us
216  a chance to test their crossing.
217  
218  == Sweeping the sweep line
219  
220  The first observation is that the sweep line invariants are maintained as we increase $y$ continuously up
221  to the next event:
222  + For the order invariant, first note that there is an event whenever $y$ leaves the domain of a segment, so $y$ remains in all domains until the next event.
223    Moreover, if at any $y' > y$ the ordering breaks, two line segments must have $epsilon$-crossed one another by $y'$.
224    The third invariant guarantees that there's an event before this happens, so by the contra-positive until an event happens the ordering
225    constraint is maintained.
226  + The crossing invariant is maintained because the set of things to check (i.e. the set of line segments that $epsilon$-cross
227    one another after $y$) only shrinks as $y$ increases.
228  
229  == Interaction with adjacent segments
230  
231  In vanilla Bentley-Ottmann, each segment gets compared to its two sweep-line neighbors; they can either intersect or not intersect.
232  When numerical errors are taken into account, we may need to compare to
233  more segments.
234  
235  #figure(
236  cetz.canvas({
237    import cetz.draw: *
238  
239    let calc_eps(x0, x1) = 0.5 * calc.sqrt(1 + calc.pow(calc.abs(x1 - x0) / 2, 2))
240  
241    let mkline(x0, x1) = {
242      let eps = calc_eps(x0, x1)
243      line((x0, 1), (x1, -1))
244      line((x0 + eps, 1), (x1 + eps, -1), stroke: ( thickness: 0.2pt, dash: "dashed" ))
245      line((x0 - eps, 1), (x1 - eps, -1), stroke: ( thickness: 0.2pt, dash: "dashed" ))
246    }
247  
248    let a0 = -8
249    let a1 = 8
250    let b0 = -2
251    let b1 = 2
252    let c0 = 0.5
253    let c1 = 0.5
254  
255    content((a0, 1), $alpha^1$, anchor: "south", padding: 5pt)
256    content((b0, 1), $alpha^2$, anchor: "south", padding: 5pt)
257    content((c0, 1), $alpha^3$, anchor: "south", padding: 5pt)
258  
259    mkline(a0, a1)
260    mkline(b0, b1)
261    mkline(c0, c1)
262  
263    let a_eps = calc_eps(a0, a1)
264    let b_eps = calc_eps(b0, b1)
265    let c_eps = calc_eps(c0, c1)
266    let cross12 = 1 - 2*(b0 + b_eps - a0 + a_eps) / (b0 - a0 - (b1 - a1))
267    let cross13 = 1 - 2*(c0 + c_eps - a0 + a_eps) / (c0 - a0 - (c1 - a1))
268  
269    line((-10, cross13), (10, cross13), stroke: ( thickness: 0.2pt ))
270    line((-10, cross12), (10, cross12), stroke: ( thickness: 0.2pt ))
271  }),
272  caption: [Three segments all mixed up. The upper horizontal line is the height at which $(alpha^1, alpha^3)$ $epsilon$-cross,
273  and the lower horizontal line is the height at which $(alpha^1, alpha^2)$ $epsilon$-cross.]
274  )<mixed-up-segments>
275  
276  For example, in @mixed-up-segments imagine that $alpha^1$ and $alpha^2$
277  are present in the sweep line, and we just added $alpha^3$. If we compare
278  $alpha^3$ to $alpha^2$, we'll see that the $epsilon$-cross and so we'll add an
279  intersection event between them sometime betfore they $epsilon$-cross. But if
280  that's all we do then we're in trouble, because by the time we come around to
281  looking at that intersection event $alpha^1$ may have already $epsilon$-crossed
282  $alpha^3$, breaking the ordering invariant. Note that $alpha^2$ doesn't
283  $epsilon$-cross $alpha^3$ at all in @mixed-up-segments, so there's no guarantee of
284  an intersection event between those two.
285  
286  Fix $y$ and $epsilon$, and
287  suppose we have a collection of lines $(alpha^i, ..., alpha^n)$ that satisfy the ordering invariant
288  at $y$, and suppose also that $(alpha^(i+1), ..., alpha^n)$
289  satisfy the crossing invariant at $y$.
290  To make the whole collection satisfy both invariants, we run the following algorithm.
291  We call this an *intersection scan to the right*.
292  
293  #figure(
294    pseudocode-list([
295      + *for* $j = i+1$ up to $n$
296        + #line-label(<w-def>) let $w^j$ be the smallest height of any event between $i$ and $j$
297  
298        + #line-label(<crossing-test>) *if* $(alpha^i, alpha^j_(+,epsilon))$ cross by $w^j$
299          + choose $z$ before the crossing, such that $alpha^i approx_(z,epsilon) alpha^j$
300          + insert an intersection event for $(alpha^i, alpha^j)$ at $z$
301  
302        + #line-label(<protect>) *if* $alpha^i (z) <= alpha^j_-(z)$ for all $z in [y, w^j]$
303          + *break*
304    ]),
305    caption: "Intersection scan to the right"
306  )
307  
308  #inexact[
309  The test at @protect can be seen as an early-stopping optimization, and is not necessary for correctness.
310  In particular, an approximation with no false positives
311  is also fine.
312  ]
313  
314  #lemma[
315  Suppose that $(alpha^i, ..., alpha^n)$ satisfy the ordering invariant
316  at $y$, and suppose also that $(alpha^(i+1), ..., alpha^n)$
317  satisfy the crossing invariant at $y$.
318  After running an intersection scan to the right,
319  $(alpha^i, ..., alpha^n)$ satisfy the crossing invariant at $y$.
320  
321  In fact, $(alpha^i, ..., alpha^n)$ satisfy a slightly stronger crossing invariant at $y$: for every $j > i$,
322  if $(alpha^i, alpha^j_(+,epsilon))$ cross then the event queue contains an event between $i$ and $j$, and before
323   the crossing height.
324  ]<lem-intersection-scan>
325  
326  (The special thing about the stronger crossing invariant is that it asks whether
327  $(alpha^i, alpha^j_(+,epsilon))$ cross, where the usual crossing invariant asks
328  whether 
329  $(alpha^i_(-,epsilon), alpha^j_(+,epsilon))$ cross.)
330  
331  #proof[
332    It suffices to check the stronger crossing invariant.
333    So take some $k > i$
334    that $(alpha^i, alpha^k_(+,epsilon))$ cross. We consider two cases: whether or not the loop terminated
335    before reaching $k$.
336  
337    - Suppose the loop terminated at some $j < k$. If $(alpha^i, alpha^k_(+,epsilon))$ cross
338      after $w^j$, then the definition of $w^j$ ensures that there is an event between $i$ and $j$ (and therefore between
339      $i$ and $k$) before the crossing. On the other hand, the termination condition ensures that
340      $alpha^i (z) <= alpha^j_-(z)$ until $w^j$, and so if $(alpha^i, alpha^k_(+,epsilon))$ cross before $w^j$ then
341      also $(alpha^j, alpha^k)$ cross before that. In this case, the crossing invariant for $(alpha^(i+1), ..., alpha^n)$
342      implies the existence of an event between $j$ and $k$ (and therefore between $i$ and $k$) before the crossing.
343    - If the loop included the case $j = k$, we break into two more cases:
344      - If $(alpha^i, alpha^k_(+,epsilon))$ cross by $w^j$, then the algorithm inserted a witnessing event between $i$ and $j$.
345      - Otherwise, the definition of $w^j$ ensures that there is an event between $i$ and $j$ before the crossing.
346  ]
347  
348  #remark[
349  We can tweak the algorithm a little to try and reduce the number of comparisons.
350  For a start, if we add an intersection event between $i$ and $j$ then we can
351  set $w^j$ to $z$, the height of that new event (because $z$ is the new smallest
352  height of any event between $i$ and $j$).
353  Therefore, if we broaden the test @crossing-test to check whether $(alpha^i,
354  alpha^j_-)$ cross, and we choose $z$ to be before the crossing of $(alpha^i,
355  alpha^j_-)$ then we can skip the test at @protect and terminate straight away.
356  Thus, there is no loop at all: we only need to consider $j = i + 1$.
357  
358  One potential issue with this is finding a height $z$ with $alpha^i approx_(z,epsilon) alpha^j$
359  that's before $(alpha^i, alpha^j_-)$ cross; in other words, we want a height after $(alpha^i_+, alpha^j_-)$
360  cross but before $(alpha^i, alpha^j_-)$ cross. If $alpha^j$ is almost horizontal, this window might be very
361  small (or maybe even not representable with our restricted precision).
362  It may still be worthwhile to have a fast path with a single comparison, and a slow path with a loop.
363  ]
364  
365  As you might have already guessed, we can also intersection scan to the left; it's pretty much a reflection
366  of the other one.
367  
368  #figure(
369    pseudocode-list([
370      + *for* $j = i$ down to $1$
371        + let $w^j$ be the smallest height of any event between $j$ and $i$
372  
373        + *if* $(alpha^j_(-,epsilon), alpha^i)$ cross by $w^j$
374          + choose $z$ before the crossing, such that $alpha^j approx_(z,epsilon) alpha^i$
375          + insert an intersection event for $(alpha^j, alpha^i)$ at $z$
376  
377        + *if* $alpha^j_+ (z) <= alpha^i (z)$ for all $z in [y, w^j]$
378          + *break*
379    ]),
380    caption: "Intersection scan to the left"
381  )
382  
383  We'll skip the proof of the following lemma, because it's basically the same as the last one.
384  
385  #lemma[
386  Suppose that $(alpha^1, ..., alpha^i)$ satisfy the ordering invariant
387  at $y$, and suppose also that $(alpha^1, ..., alpha^(i+1))$
388  satisfy the crossing invariant at $y$.
389  After running an intersection scan to the left,
390  $(alpha^1, ..., alpha^(i+1))$ satisfy a slightly stronger crossing invariant at $y$: for every $j <= i$,
391  if $(alpha^j_(-,epsilon), alpha^(i+1))$ cross then the event queue contains an event between $j$ and $i+1$, and before
392   the crossing height.
393  ]<lem-intersection-scan-left>
394  
395  The purpose of the stronger crossing invariant comes in the next lemma, which deals with scanning in both directions
396  and allows the insertion of a segment in the middle of a sweep-line.
397  
398  #lemma[
399  Suppose that $(alpha^1, ..., alpha^n)$ satisfy the ordering invariant at $y$, and suppose that
400  $(alpha^1, ..., alpha^i)$ and $(alpha^(i+1), ..., alpha^n)$ satisfy the crossing invariant at $y$.
401  Let $beta$ be another segment and assume that $(alpha^1, ..., alpha^i, beta, alpha^(i+1), ... alpha^n)$
402  satisfy the ordering invariant at $y$. After running an intersection scan to the left and an intersection
403  scan to the right from $beta$, 
404  $(alpha^1, ..., alpha^i, beta, alpha^(i+1), ... alpha^n)$ satisfies the crossing invariant at $y$.
405  ]<lem-intersection-scan-bidirectional>
406  
407  #proof[
408  @lem-intersection-scan implies that $(beta, alpha^(i+1), dots, alpha^n)$ satisfies the crossing invariant,
409  and
410  @lem-intersection-scan-left implies that $(alpha^1, ..., alpha^i, beta)$ does also. It only remains
411  to consider interactions between a segment before $beta$ and one after. So fix $j <= i < k$ and suppose
412  that $(alpha^j, alpha^k)$ $epsilon$-cross. If they $epsilon$-cross after $y_1(beta)$ then $beta$ exit
413  event witnesses the crossing invariant, so assume that $(alpha^j, alpha^k)$ $epsilon$-cross by $y_1(beta)$.
414  Then $(alpha^j_-, alpha^k_+)$ cross by $y_1(beta)$, and so one of them crosses $beta$ before $(alpha^j, alpha^k)$ $epsilon$-cross.
415  If $(alpha^j_-, beta)$ cross then @lem-intersection-scan-left implies that there is an event between $alpha^j$ and $beta$ (and
416  therefore between $alpha^j$) and $alpha^k$ before the crossing height; otherwise, $(beta, alpha^k_+)$ cross
417  and so @lem-intersection-scan provides the required event.
418  ]
419  
420  One last observation in this section, that follows trivially from the algorithm:
421  
422  #lemma[
423  If an intersection scan inserts an intersection event for $(alpha, beta)$
424  then the intersection event's height $z$ satisfies $alpha approx_(z,epsilon) beta$.
425  ]<lem-valid-intersection-events>
426  
427  == An "enter" event
428  
429  When inserting a new segment into the current sweep-line, we first choose its sweep-line position using
430  a binary search on its horizontal coordinate. Let's write $(alpha^1, dots, alpha^m)$ for the sweep-line
431  before inserting the new segment, and let's call the new segment $beta$. First, we observe that
432  there is a place to insert the new segment while preserving the ordering invariant.
433  
434  #lemma[
435  Suppose $(alpha^1, ..., alpha^m)$ is $epsilon$-ordered at $y$, and let $i$ be the largest $j$ for which
436  $alpha^j prec_(y,epsilon) beta$. Then
437  $(alpha^1, ..., alpha^i, beta, alpha^(i+1), ..., alpha^m)$ is $epsilon$-ordered at $y$.
438  (Here, we can allow the corner cases $i = 0$ and $i = m$ by declaring that
439  "$alpha^0$" is a vertical line at $x = -infinity$ and "$alpha^(m+1)$" is a vertical line at $x = infinity$).
440  ]<lem-insert-preserving-order>
441  
442  #proof[
443  Since $(alpha^1, ..., alpha^m)$ was $epsilon$-ordered at $y$, it suffices to compare beta with all $alpha^j$.
444  For $i + 1 <= j <= m$, our choice of $i$ immediately implies that $alpha^j succ.tilde_(y,epsilon) beta$.
445  So consider $1 <= j <= i$. Since $(alpha^1, ..., alpha^m)$ is $epsilon$-ordered, $alpha^j prec.tilde_(y,epsilon) alpha^i$.
446  Since $alpha^i prec_(y,epsilon) beta$, @lem-basic-order-properties implies that $alpha^j prec.tilde_(y,epsilon) beta$.
447  ]
448  
449  #inexact[
450  @lem-insert-preserving-order guarantees the existence of an insertion point, but it doesn't say how to
451  find it efficiently.
452  But consider a predicate $f(alpha^j)$ that returns true whenever $alpha^j prec_(y,epsilon) beta$, and
453  false whenever $alpha^j_+(y) > beta_+(y)$. Running a binary search with this predicate will find some $i$
454  for which $f(alpha^i)$ is true and $f(alpha^(i+1))$ is false. By scanning to the right from there, we can
455  find the largest such $i$.
456  
457  This $i$ is at least as large as the $i$ in @lem-insert-preserving-order, so
458  to check that it's a valid insertion point we only need to check that it isn't too large. So if $1 <= j <= i$ then
459  $alpha^j prec.tilde_(y,epsilon) alpha^i$ and so $alpha^j_-(y) <= alpha^i_+(y)$. On the other hand, $f(alpha^i)$
460  was true and so $alpha^i_+ <= beta_+(y)$. Putting these together shows that $alpha^j prec.tilde_(y,epsilon) beta$.
461  
462  This algorithm can be implemented with approximate arithmetic, and its running time is logarithmic in the
463  total length of the sweep line, plus linear in the number of elements that are very close to $beta$.
464  ]
465  
466  @lem-insert-preserving-order implies that we can insert a new segment while preserving the ordering invariant. By
467  @lem-intersection-scan-bidirectional, running an intersection scan restores the crossing invariant.
468  Thus, we can insert a new segment while preserving the sweep-line invariants.
469  
470  == An "exit" event
471  
472  When a segments exits the sweep-line, the ordering invariant clearly doesn't break.
473  Regarding the crossing invariant, it can only break because of $epsilon$-crossing pairs whose
474  crossing invariant was witnessed by the exit event that was just processed.
475  To restore the crossing invariant, we need to enqueue some new intersection events.
476  
477  Let $(alpha^1, ..., alpha^m)$ be the sweep-line after removing the just-exited segment
478  which, we assume, used to live between $alpha^i$ and $alpha^(i+1)$. Note that both
479  $(alpha^1, ..., alpha^i)$ and $(alpha^(i+1), ..., alpha^n)$ satisfy the crossing invariant.
480  By @lem-intersection-scan-bidirectional, running an intersection scan from $alpha^i$ in both
481  directions restores the crossing invariant. (Technically, this isn't covered by the statement
482  of @lem-intersection-scan-bidirectional -- which involves a new segment $beta$ -- but the proof
483  is basically the same.)
484  
485  == An "intersection" event
486  
487  Suppose our sweep-line is $(alpha^1, ..., alpha^m)$ and we've just encountered an intersection event for $(alpha^i, alpha^j)$
488  at height $y$.
489  If $i > j$ then they've already been swapped in our sweep-line, so we don't need to swap them again. If $i < j$, we need to
490  swap them. According to @lem-valid-intersection-events, $alpha^j prec.tilde_(y,epsilon) alpha^i$. It seems reasonable, therefore,
491  to reorder the sweep line by putting $alpha^i$ after $alpha^j$, like
492  
493  $
494  alpha^1, ..., alpha^(i-1), alpha^(i+1), ..., alpha^j, alpha^i, alpha^(j+1), ... alpha^n.
495  $
496  
497  The issue with this is that $prec.tilde_(y,epsilon)$ is that it's changed the order of pairs other than $alpha^i$ and $alpha^j$:
498  for every $i < k < j$, the ordering between $alpha^i$ and $alpha^k$ has been swapped. If $prec.tilde$ were transitive, this would
499  be fine. But it isn't.
500  
501  To fix this issue, we allow $alpha^i$ to "push" some of the intermediate segments along with it. It's a bit tedious to write (and read)
502  this precisely, so hopefully this description is enough: for each $k$ between $i$ and $j$, if $alpha^i prec_(y,epsilon) alpha^k$ then
503  we also move $alpha^k$ after $alpha^j$. Also, we preserve the relative order of all segments moved in this way.
504  To see that this "pushing" maintains the ordering invariants, note that by definition it preserves the ordering for all comparisons
505  with $alpha^i$: if some $alpha^ell$ was pushed along with $alpha^i$ then their relative orders haven't changed; and if $alpha^ell$
506  wasn't pushed then $alpha^ell prec.tilde_(y,epsilon) alpha^i$ and the new order is fine.
507  
508  What about other pairs? If $alpha^k$ and $alpha^ell$
509  changed orders then one of them ($alpha^ell$, say) was pushed and the other wasn't. Then $alpha^i prec_(y,epsilon) alpha^ell$
510  by our choice of the pushed segments, and $alpha^k prec.tilde_(y,epsilon) alpha^i$ by the previous paragraph. Putting these
511  together, @lem-basic-order-properties implies that $alpha^k prec.tilde_(y,epsilon) alpha^ell$ and so the new order is ok.
512  
513  #inexact[
514  We might not be able to tell exactly which segments $alpha^k$ between $alpha^i$ and $alpha^j$ satisfy
515  $alpha^i prec_(y,epsilon) alpha^k$. Fortunately, we can push a few too many segments while still maintaining
516  correctness: it suffices to
517  include all segments $alpha^k succ_(y,epsilon) alpha^i$, while also ensuring that we only include $alpha^k$ for which
518  $alpha^k_+(y) >= alpha^i_+(y)$.
519  ]
520  
521  Finally, we need to maintain the crossing invariant. We can probably be more efficient about this, but one
522  way to be correct is to treat the various swappings as a bunch of deletions from the sweep-line followed by
523  a bunch of insertions into the sweep-line. We have already shown that running an intersection scan
524  after each insertion and deletion correctly maintains the crossing invariant, so we can just do that.
525  
526  = Correctness of the output
527  
528  We've described how to turn a bunch of segments into a continuum of sweep-lines, but we probably actually wanted
529  to find the segments' intersection points. Let's define exactly what that means and how to get there.
530  We'll assume that we start with one or more polylines that may or may not be closed. Equivalently, each segment
531  comes with an orientation and (optionally) a "following" segment whose starting point is the ending point of the
532  current segment.
533  Our goal is to subdivide all segments at all intersection points, meaning that for each segment we'll produce
534  a polyline; we'll call the segments in each such polyline "output segments," to distinguish them from
535  the original segments. We require three properties:
536  
537  - "approximation": for each input segment, every vertex in its output polyline must be within $epsilon$ of the original segment.
538    Also the first and last vertices in the output polyline must be within $epsilon$ of the start and end of the segment, respectively.
539  - "continuity": if a segment has a following segment, then the last vertex in its output polyline must coincide
540    with the first vertex of the following segment's output polyline.
541  - "completeness": output segments intersect only at their endpoints, unless they are identical.
542  
543  In the solution we describe, we'll achieve both continuity and the second half
544  of approximation by insisting that output polylines have exactly the same start-
545  and endpoints as their input segment. This stronger condition makes the algorithm
546  simpler, but it also subjectively can hurt output "quality" by requiring more segments.
547  
548  == The "dense" version
549  
550  Suppose we've already found a weakly-ordered sweep line at every height. We say
551  a height is "important" if the sweep line changed at that height, either because
552  some segments entered or exited, or because the order changed. The naive version
553  of our subdivision algorithm will subdivide every segment at every important
554  height. Note that because the sweep line changes at a discrete set of heights,
555  at every important height $y$ we actually have two weakly-ordered sweep line:
556  the "old" one from infinitesmally before $y$ height and the "new" one from $y$
557  onwards. Because everything is continuous and the weak ordering conditions are
558  closed, both the old and new sweep lines satisfy the weak ordering condi0.2tions
559  at $y$.
560  
561  The algorithm goes like this: at every important height we use
562  @lem-horizontal-realization to assign horizontal positions to segments in the
563  old sweep line; we subdivide each segment at the resulting coordinate. Then we
564  use @lem-horizontal-realization again to assign horizontal positions to segments
565  in the new sweep line. If a segment gets assigned two different horizontal positions
566  by the old and new sweep lines, we add a horizontal segment between them. Also,
567  if a segment starts or ends at the current height and its assigned position is not the
568  same as its starting or ending position, we fix that up with another horizontal segment.
569  Finally, we subdivide all horizontal segments as necessary to ensure that they only
570  intersect at endpoints.
571  
572  I think it's pretty clear that this algorithm is correct. The approximation property
573  holds (with $sqrt 2 epsilon$ instead of $epsilon$, but that doesn't matter) because
574  @lem-horizontal-realization always puts a new point within its segment's upper and lower bounds,
575  and the bounds are chosen to be within $sqrt 2 epsilon$ of the segment; see @fig-error-bars.
576  As mentioned above, the "continuity" property and the other half of "approximation" follow
577  because we insisted on including the segment endpoints exactly.
578  
579  For the "completeness" property, let's focus on the bit between the important heights
580  (because everything gets subdivided at the important heights, and so maintaining the
581  completeness property there is just a matter of splitting into enough horizontal segments).
582  Between two important heights, the sweep line is constant, so the new sweep line
583  at the old height is the same as the old sweep line at the new height. Thus, all output
584  subsegments have the same order when leaving the old sweep line as they do when entering
585  the new sweep line. Therefore any two of them that intersect between the two important
586  heights must be identical.
587  
588  == The "sparse" version
589  
590  We'd like to avoid subdividing every segment at every important height. This basically
591  involves detecting which segments need to be divided and ignoring the rest. It's implemented
592  but not yet written up (TODO).
593  
594  == A more detailed insertion algorithm
595  
596  It took me a few tries to get the approximate insertion algorithm working, so
597  I figured it was worth writing up a few more details. First of all, a slight
598  tweak of @lem-insert-preserving-order shows that there is a position $i$ such
599  that $alpha^j_- (y) <= beta(y)$ for all $j <= i$ and $alpha^j_+ (y) >= beta(y)$
600  for all $j > i$. We will look for this $i$, but with some slight slack in our
601  arithmetic. The slack will be small enough to guarantee we find a position that
602  satisfies the weaker properties of @lem-insert-preserving-order.
603  
604  First, we consider a predicate $p(j)$ that returns `false` whenever $alpha^j_+
605  (y) > beta(y)$. We put no conditions on the values of $j$ for which $p(j)$
606  returns `true`; a tighter predicate will give a more efficient algorithm, but
607  doesn't affect correctness. Let $i_-$ be such that $p(i_-) = #true$ and $p(i_-
608  + 1) = #false$ (where we allow $i_- = -1$, and handle the boundary cases by
609  declaring that $p(-1) = #true$ and $p(m+1) = #false)$. Note that $i_-$ can be
610  found with a binary search.
611  
612  The definition of $p$ ensures that $alpha^(i_-)_+(y) <= beta(y)$. Under the
613  assumption that $(alpha^1, ..., alpha^m)$ is weakly ordered, it follows that
614  for every $j <= i_-$, $alpha^j_-(y) <= beta(y)$. That is, the $i$ that we are
615  looking for is larger than or equal to $i_-$.
616  
617  The next step is a linear scan up from $i_-$. Let $q(j)$ be a predicate
618  that returns `true` whenever $alpha^j_- (y) < beta^j (y)$ and `false` whenever
619  $alpha^j_- (y) > beta^j_+ (y)$. Take $i >= i_-$ to be the index just before $q$
620  first returns false; we will insert $beta$ after $alpha^i$.
621  
622  Let's check that this insertion position is valid: for every $j <= i$ we have
623  either $j <= i_-$ or $i_- < j <= i$. In the first case, we already observed
624  that $alpha^j_-(y) <= beta(y)$; in the second case, $q(j)$ returned `true` and
625  so $alpha^j_-(y) <= beta^j^+ (y)$. In either case, $alpha^j prec.tilde beta$.
626  On the other hand, $q(i+i)$ being `false` implies that $alpha^(i+1)_- (y) >=
627  beta(y)$. This implies that for all $j >= i+1$, $alpha^j_+ (y) >= beta(y)$ and
628  so $alpha^j succ.tilde_y beta$.
629  
630  = Accuracy analysis
631  
632  We're going to implement some algorithms in floating point. Suppose $p$ is the number of mantissa bits
633  of the floating-point type, and let $fl(t)$ denote the rounding function. Let $epsilon_0 = 2^(-1-p)$;
634  this is the relative error of the rounding function for all normal numbers: for a real number $x$,
635  as long as $fl(x)$ is finite and non-subnormal,
636  $
637  fl(x) = x plus.minus epsilon_0 |x|.
638  $
639  Now, if the result of floating point addition or subtraction is subnormal then it is exact; therefore,
640  for any floating-point $x$ and $y$, if $x - y$ and $x + y$ don't overflow then after rounding they
641  have relative error at most $epsilon_0$.
642  
643  #lemma[
644    If $y_0 <= y <= y_1$ are floating-point numbers and $epsilon_0 <= 1/2$ then computing $(y - y_0) / (y_1 - y_0)$
645    the naive way gives an absolute error of
646    at most $3 epsilon_0 + 4 epsilon_0^2$.
647  ]<lem-ratio-error>
648  
649  #proof[
650  The assumption $y_0 <= y <= y_1$ ensures that $y - y_0$ and $y_1 - y_0$ do not overflow. Therefore
651  
652  $
653  fl(y - y_0) &= (y - y_0) (1 plus.minus epsilon_0) \
654  fl(y + y_0) &= (y + y_0) (1 plus.minus epsilon_0)
655  $
656  and so
657  $
658  fl(y - y_0) /
659  fl(y_1 - y_0) = (y - y_0) / (y_1 - y_0) times
660   (1 plus.minus epsilon_0) / (1 plus.minus epsilon_0).
661  $
662  For the upper bound,
663  $
664   (1 + epsilon_0) / (1 - epsilon_0) = (1 + epsilon_0) (1 + epsilon_0 + epsilon_0^2 / (1 - epsilon_0))
665  = 1 + 2 epsilon_0 + 2 / (1 - epsilon_0) epsilon_0^2 <= 1 + 2 epsilon_0 + 4 epsilon_0^2.
666  $
667  For the lower bound,
668  $
669   (1 - epsilon_0) / (1 + epsilon_0) = 1 - epsilon_0 - (epsilon_0 (1-epsilon_0)) / (1+epsilon_0) >= 1 - 2 epsilon_0.
670  $
671  Putting these together,
672  $
673  fl(y - y_0) /
674  fl(y_1 - y_0) = (y - y_0) / (y_1 - y_0) (1 plus.minus (2 epsilon_0 + 4 epsilon_0^2)).
675  $
676  Since $0 <= (y - y_0) / (y_1 - y_0) <= 1$, this relative error is also an absolute error:
677  $
678  fl(y - y_0) /
679  fl(y_1 - y_0) = (y - y_0) / (y_1 - y_0) plus.minus (2 epsilon_0 + 4 epsilon_0^2).
680  $
681  The final rounding step adds another absolute error of at most $epsilon_0$.
682  ]
683  
684  #lemma[
685  Suppose that $x_0$, $x_1$, and $M$ are floating-point numbers with $|x_0|, |x_1| <= M/2$.
686  Let $0 <= t <= 1$ be a real number, and let $0 <= hat(t) <= 1$ be a floating point number
687  with $|t - hat(t)| <= epsilon$. Then computing $x_0 + t (x_1 - x_0)$ the obvious way gives an
688  absolute error of at most $(3 epsilon_0 + epsilon) M$.
689  ]<lem-convex-combination-error>
690  
691  #proof[
692  Let's first consider the error introduced by approximating $t (x_1 - x_0)$ with $hat(t) fl(x_1 - x_0)$:
693  $
694  hat(t) fl(x_1 - x_0) = hat(t) (x_1 - x_0)(1 plus.minus epsilon_0) = hat(t) (x_1 - x_0) plus.minus epsilon_0 M,
695  $
696  where the second approximation follows because $0 <= hat(t) <= 1$ and $|x_1 - x_0| <= M$.
697  Since $hat(t) = t plus.minus epsilon$, we continue with
698  $
699  hat(t) fl(x_1 - x_0)
700  = (t plus.minus epsilon) (x_1 - x_0) plus.minus epsilon_0 M,
701  = t (x_1 - x_0) plus.minus (epsilon_0 + epsilon) M.
702  $
703  Rounding the left hand side introduces a relative error of at most $epsilon_0$, which is an absolute error of at most $M epsilon_0$. Therefore,
704  $
705  fl(hat(t) fl(x_1 - x_0))
706  = t (x_1 - x_0) plus.minus (2 epsilon_0 + epsilon) M.
707  $
708  
709  Finally, we subtract both sides from $x_0$ and round one more time. This last rounding introduces a relative error of at most $epsilon_0$,
710  which again leads to an absolute error of at most $2 epsilon_0 M$ (and I think maybe the 2 is unnecessary, but let's not worry).
711  ]
712  
713  Next, let's analyze the computation of the $y$ intersection point.
714  
715  #lemma[
716  Suppose we have two segments $alpha$ and $beta$. Set $y_0 = y_0(alpha) or y_0(beta)$ and $y_1 = y_1(alpha) and y_1(beta)$
717  (these can be computed exactly). Suppose that we can compute $alpha(y)$ and $beta(y)$ to additive accuracy $epsilon$,
718  and assume that $beta(y_0) - alpha(y_0) >= 0$ and
719  $
720  alpha(y_1) - beta(y_1) > 16/3 (epsilon + 2 epsilon_0 M).
721  $
722  Define $t = (alpha(y_1) - beta(y_1)) / (alpha(y_1) - beta(y_1) + beta(y_0) - alpha(y_0)$ and set $y = y_0 + t(y_1 - y_0)$. The
723  natural way to compute $y$ produces an approximation $hat(y)$ satisfying
724  $
725  |alpha(hat(y)) - beta(hat(y))| <= 4(epsilon + 2 epsilon_0 M) (1 + max(1/m_alpha, 1/m_beta)),
726  $
727  where $m_alpha$ and $m_beta$ are the slopes of $alpha$ and $beta$.
728  ]<lem-intersection-height-error>
729  
730  Basically, this says we can find a good approximate crossing height as long as the two segments cross by at least $16/3 (epsilon + 2 epsilon_0 M)$.
731  
732  #proof[
733  First, note that at least one of $alpha(y_0)$ or $beta(y_0)$ is calculated exactly; similarly for $alpha(y_1)$ and $beta(y_1)$.
734  Therefore the absolute error in calculating $alpha(y_1) - beta(y_1)$ is at most $epsilon + epsilon_0 M$ (where the $epsilon_0 M$ comes
735  from rounding after the summation) and the absolute error in calculating
736  $(alpha(y_1) - beta(y_1)) + (beta(y_0) - alpha(y_0))$ is at most $2 epsilon + 3 epsilon_0 M$.
737  These errors aren't independent, though, since both terms involve a computation of $alpha(y_1) - beta(y_1)$.
738  Specifically, let $c = alpha(y_1) - beta(y_1)$ and let $d = (alpha(y_1) - beta(y_1)) + (beta(y_0) - alpha(y_0))$; suppose
739  our calculation of $c$ has error $epsilon_c$ and our calculation of $d$ has error $epsilon_d$. Then
740  
741  $
742  |epsilon_c - epsilon_d| <= epsilon + 2 epsilon_0 M,
743  $ <eq-error-cancellation>
744  
745  because this is the *additional* error introduced into $d$ by computing $beta(y_0) - alpha(y_0)$ and adding it to the other term.
746  
747  Now we consider the quotient
748  
749  $
750  (c + epsilon_c) / (d + epsilon_d)
751  &= c / (d (1 + epsilon_d/d)) + epsilon_c / (d (1 + epsilon_d/d)) \
752  &= c / d (1 - epsilon_d/d + (epsilon_d/d)^2 / (1 + epsilon_d/d)) + epsilon_c / d (1 - (epsilon_d/d)/(1 + epsilon_d / d)).
753  $
754  
755  We've assumed that $1 + epsilon_d/d >= 1/2$, and so the error
756  $
757  abs((c + epsilon_c) / (d + epsilon_d) - c/d)
758  $
759  is at most
760  $
761  & abs(c / d ( - epsilon_d/d + (epsilon_d/d)^2 / (1 + epsilon_d/d)) + epsilon_c / d (1 - (epsilon_d/d)/(1 + epsilon_d / d))) \
762  & <= abs(- c / d epsilon_d/d + epsilon_c / d ) + 2 c/d (epsilon_d/d)^2 + 2 abs(epsilon_c epsilon_d) / d^2 \
763  & <= c/d abs(epsilon_c/d - epsilon_d/d) + (1 - c/d) abs(epsilon_c/d) + 2 (epsilon_d/d)^2 + 2 abs(epsilon_c epsilon_d) / d^2,
764  $
765  where in the last line we recall that $0 <= c/d <= 1$. Now, recall from @eq-error-cancellation and the paragraph before it
766  that both $abs(epsilon_c)$ and $abs(epsilon_c - epsilon_d)$ are bounded by $epsilon + 2 epsilon_0 M$, and that
767  $abs(epsilon_d)$ is bounded by $2 epsilon + 3 epsilon_0 M$. Therefore, our computation of $c/d$ (which, recall, we called $t$ above) has absolute error of at most
768  
769  $
770  (epsilon + 2 epsilon_0 M)/d + 4 (2 epsilon + 3 epsilon_0 M)^2 / d^2.
771  $
772  
773  Let's call this expression $epsilon_t$.
774  
775  By @lem-convex-combination-error, the computation of $y = y_0 + t (y_1 - y_0)$ has absolute error at most $(2 epsilon_0 + epsilon_t) |y_1 - y_0| + epsilon_0 M$.
776  Let's see what this means for the horizontal accuracy of $alpha(y)$ and $beta(y)$. We can write
777  $
778  alpha(y) = alpha(y_0) + (y - y_0) / (y_1 - y_0) (alpha(y_1) - alpha(y_0)),
779  $
780  and similarly for $beta(y)$. Combining these and recalling our definitions of $c$ and $d$,
781  $
782  alpha(y) - beta(y) = -c + (y - y_0) / (y_1 - y_0) d.
783  $
784  For the true value of $y$, this is zero of course. With our error bound of $(2 epsilon_0 + epsilon_t) |y_1 - y_0| + epsilon_0 M$ on our
785  approximation of $y$ (let's call it $hat(y)$),
786  $
787  abs(alpha(hat(y)) - beta(hat(y))) <= (2 epsilon_0 + epsilon_t) d + (epsilon_0 M) / (y_1 - y_0) d
788  <= epsilon_t d + 3 epsilon_0 (M d) / (y_1 - y_0).
789  $
790  
791  If we assume that $16 (epsilon + 2 epsilon_0 M) / d <= 3$ then $epsilon_t d <= 4 epsilon + 8 epsilon_0 M$. 
792  Since $d/(y_1 - y_1)$ is the difference of the inverse-slopes of $alpha$ and $beta$, it's bounded by twice
793  the maximum of these inverse-slopes.
794  ]
795  
796  Putting all this together (and handling the "chamfers" yet):
797  #lemma[
798  Assume all line segments are inside the square $[-M/2, M/2]^2$,
799  and take $delta = 64 epsilon_0 M$. If $epsilon_0 <= 1/4$ then
800  
801  - We can compute $x$ coordinates of line segments to additive accuracy $delta/8$.
802  - If $alpha$ $(3 delta)/4$-crosses $beta$, we can compute a crossing height $hat(y)$ such
803    that
804    $
805    |alpha(hat(y)) - beta(hat(y))| <= (9 delta) /16 (1 + max(1/m_alpha, 1/m_beta)).
806    $
807  ]
808  
809  #proof[
810  We compute $x$ coordinates by first computing $t = (y - y_0) / (y_1 - y_0)$ (with error at most $3 epsilon_0 + 4 epsilon_0^2$ by
811  @lem-ratio-error). Then we apply @lem-convex-combination-error with $epsilon = 3 epsilon_0 + 4 epsilon_0^2$ to see that we
812  can compute the horizontal coordinate with error at most $(6 epsilon_0 + 4 epsilon_0^2)M$. If $epsilon_0 <= 1/4$ as we assumed,
813  this is at most $7 epsilon_0 M <= delta/8$. This proves the first claim.
814  
815  For the second claim, we apply @lem-intersection-height-error with $epsilon = (6 epsilon_0 + 4 epsilon_0^2) M$.
816  With this $epsilon$, $epsilon + 2 epsilon_0 M) <= 9 epsilon_0 M$ and so $16/3 (epsilon + 2 epsilon_0 M) <= 48 epsilon_0 M = (3 delta)/4$.
817  Thus, if $alpha$ $(3 delta)/4$-crosses $beta$ then the assumptions of @lem-intersection-height-error are satisfied
818  (and in the conclusion, $4 (epsilon + 2 epsilon_0) M <= 36 epsilon_0 M = (9 delta) /16 $.
819  ]
820  
821  = Curves etc.
822  
823  Let's move on from lines and talk about curves. Each curve is of the form $alpha: [0, 1] -> bb(R)^2$, which we'll sometimes
824  in components as $alpha_0: [0, 1] -> bb(R)$ and $alpha_1: [0, 1] -> bb(R)$. We'll assume that each $alpha_0$ is a strictly
825  increasing function. Write $S = [-1, 1] times [-1, 1]$ for the unit square in $bb(R)^2$. For a curve $alpha$ and a set $B subset bb(R)^2$, we write
826  $alpha + B$ for the set
827  
828  $
829  alpha + B = {(x, y) in bb(R)^2 : exists t in [0, 1] "with" (x, y) - alpha(t) in B}
830  $
831  
832  In other words, $alpha + B$ is the Minkowski sum between $B$ and the locus of $alpha$. For any set $B subset bb(R)^2$,
833  let
834  
835  $
836  B_y = {x in bb(R) : (x, y) in B}.
837  $
838  
839  For example, here is a picture of $(alpha + epsilon B)_y$: (TODO)
840  
841  == The curve orders
842  
843  When moving to curves, it seems useful to work with more abstract orders. For fixed $epsilon$ and $y$, we'll define
844  three of them: the "weak pre-order," the "strong order," and the "weak order."
845  
846  #let weakprelt = sym.lt.tri
847  #let weakpreleq = sym.lt.eq.tri
848  #let notweakprelt = sym.lt.tri.not
849  #let weakpreeq = sym.eq.dots
850  #let weaklt = sym.lt
851  #let stronglt = sym.lt.double
852  #let weaklt = sym.lt
853  #let weakeq = sym.tilde
854  
855  We begin with the weak pre-order.
856  
857  #def[
858    For each $epsilon >= 0$ and $y in bb(R)$, let $weakprelt_(y,epsilon)$ be a relation satisfying:
859  
860    - if $alpha weakprelt_(y,epsilon) beta$ then $(alpha + epsilon/2 S)_y <= (beta + epsilon/2 S)_y$,
861    - if $(alpha + epsilon S)_y <= (beta + epsilon S)_y$ then $alpha weakprelt_(y,epsilon) beta$,
862    - for any $k >= 1$, there are no $alpha_1, ..., alpha_k$ with $alpha_1 weakprelt_(y,epsilon) alpha_2 weakprelt_(y,epsilon) dots.h.c weakprelt_(y,epsilon) alpha_k weakprelt_(y,epsilon) alpha_1$.
863  
864    We call this last condition the "no cycles" condition.
865  ]
866  
867  Note that we *do not* assume transitivity: $alpha weakprelt_(y,epsilon) beta
868  weakprelt_(y,epsilon) gamma$ does not imply that $alpha weakprelt_(y,epsilon)
869  gamma$. We also do not assume anti-symmetry: it's possible to have neither
870  $alpha weakprelt_(y,epsilon) beta$ nor $beta weakprelt_(y,epsilon) alpha$. In
871  fact, if $(alpha + epsilon/2)_y$ and $(beta + epsilon/2 S)_y$ overlap, it is
872  guaranteed that neither of these orderings hold.
873  
874  #def[
875    If neither $alpha weakprelt_(y,epsilon) beta$ nor $beta weakprelt_(y,epsilon) alpha$, we say that
876  
877    $
878    alpha weakpreeq_(y,epsilon) beta.
879    $
880  
881    If either $alpha weakprelt_(y,epsilon) beta$ or 
882    $alpha weakpreeq_(y,epsilon) beta$, we write
883    $alpha weakpreleq_(y,epsilon) beta$. Equivalently, 
884    $alpha weakpreleq_(y,epsilon) beta$ is the negation of 
885    $beta weakprelt_(y,epsilon) alpha$.
886  ]
887  
888  TODO: draw a picture.
889  
890  The strong order is similar to the weak pre-order, but with bigger epsilons.
891  #def[
892    For each $epsilon >= 0$ and $y in bb(R)$, let $stronglt_(y,epsilon)$ be a relation satisfying:
893  
894    - if $alpha stronglt_(y,epsilon) beta$ then $(alpha + 2 epsilon S)_y <= (beta + 2 epsilon S)_y$,
895    - if $(alpha + 3 epsilon S)_y <= (beta + 3 epsilon S)_y$ then $alpha stronglt_(y,epsilon) beta$,
896    - for any $k >= 1$, there are no $alpha_1, ..., alpha_k$ with $alpha_1 stronglt_(y,epsilon) alpha_2 stronglt_(y,epsilon) dots.h.c stronglt_(y,epsilon) alpha_k stronglt_(y,epsilon) alpha_1$.
897  ]
898  
899  The definition of the weak order is a little strange at first. The initial idea
900  is that we want to ensure that our sweep-line is ordered by the weak pre-order:
901  we are going to try to ensure that if our sweep-line at $y$ is $alpha^1, ..., alpha^n$
902  then for all $i < j$, $alpha^i weakpreleq_y alpha^j$. Or in other words, we want to
903  ensure that whenever $alpha^j$ is definitely to the right of $alpha^i$, it comes after $alpha^i$
904  in the sweep-line.
905  
906  For sweep-line algorithms to be efficient, we need to avoid too many segment comparisons.
907  In the classical Bentley-Ottmann algorithm with exact arithmetic this is easy: whenever
908  we look for intersections against segment $alpha^i$, we compare it with its sweep-line
909  neighbors to the left and right (so, $alpha^(i-1)$ and $alpha^(i+1)$ in our notation).
910  Inexact arithmetic will make this a little trickier; we will typically need to look further
911  than just the immediate neighbors of $alpha^i$. However, we will need *some* mechanism
912  for early stopping, some way to say that after we have compared $alpha^i$ to, say, $alpha^(i-k)$
913  through $alpha^(i-1)$, we can be sure that no other $alpha^j$ (for $j < i - k$) needs
914  to be checked. The intuitive idea is that if $alpha^(i-k)$ is very far to the left of $alpha^i$
915  then we can stop, because then $alpha^(i-k)$ will "shield" $alpha^i$ from segments further
916  left, in much the same way that its immediate neighbors "shield" $alpha^i$ in the classical
917  sweep-line algorithm.
918  
919  We would like to use the strong order $stronglt_y$ to check whether $alpha^(i-k)$ is "very far
920  to the left" of $alpha^i$, but we run into trouble with examples like this:
921  
922  TODO: picture.
923  
924  Here, we clearly see that $alpha$ is to the left of $beta$; in situations like this, if we're
925  examining $beta$ for potential intersections and we see $alpha$ to its left, we'd really like
926  to stop looking further left. But $gamma$ is a problem: it's almost horizontal, and the
927  sweep-line at $y$ is at an ambiguous height, being more than $epsilon / 2$ but less than $epsilon$ below $gamma$.
928  According to the definition of our weak pre-order,
929  $alpha weakprelt_y gamma$ and $alpha weakpreeq_y gamma$ are
930  both allowed, as are both
931  $beta weakprelt_y gamma$ and $beta weakpreeq_y gamma$. So imagine that $alpha weakpreeq_y gamma$ and $beta weakprelt_y gamma$.
932  Before adding $beta$, the sweep-line could have been $(gamma, alpha)$. Then we could have tried to add $beta$ after $alpha$
933  (because clearly $beta$ is to the right of $alpha$), leading to the sweep-line $(gamma, alpha, beta)$.
934  Scanning for any intersections with the newly-added $beta$ won't find any because we'll look to the left, see
935  $alpha$, and stop. The result of all this is the illegal sweep-line $(gamma, alpha, beta)$, in which
936  $beta weakprelt_y gamma$ but $gamma$ comes first in the sweep-line.
937  
938  Basically, it seems hard to have both early stopping (which we want for efficiency) and the
939  ordering invariants that we need for correctness. Our "fix" is to change the definition of
940  our ordering (while still keeping it tight enough for the correctness guarantees that we need).
941  It looks almost like cheating, because we're going to redefine it precisely so that our early
942  stopping will work.
943  
944  #def[
945    Let $cal(S)_y$ be the set of all segments whose vertical range contains $y$. Say that $alpha weaklt_y beta$
946    if $alpha weakprelt_y beta$ and there do not exist $alpha', beta' in cal(S)_y$ with $beta weakpreleq_y beta' stronglt alpha' weakpreleq_y alpha$.
947  ]<def-weak-order>
948  
949  In other words, we define $weaklt$ to be the same as $weakprelt$, but if there's ever a situation where the (not-yet-specified)
950  early stopping algorithm could miss a wrongly ordered pair $alpha weakprelt beta$ by witnessing a gap between $beta'$ and $alpha'$,
951  then we weaken the $weakprelt$ order and declare instead that $alpha weakeq beta$ instead of $alpha weaklt beta$.
952  To see how this solves our tricky example above, note that when instantiating @def-weak-order with $alpha = beta$, $beta = gamma$,
953  $beta' = alpha$, and $alpha' = beta$, we see that there *do* exist $alpha', beta' in cal(S)_y$ with the required properties,
954  and so in fact $beta$ and $gamma$ are not strictly ordered in the weak order $weaklt$ even though they were strictly ordered
955  in the weak pre-order $weakprelt$. In particular, the sweep-line $(gamma, alpha, beta)$ is legal with respect to the weak order $weaklt$.
956  
957  You might object to @def-weak-order because of how hilariously inefficient it will be to implement. The
958  clever part, though, is that we won't need to implement it. Our algorithm will make its decisions based
959  only on the weak pre-order $weakprelt$ and the strong order $stronglt$. The inefficient weak order $weaklt$
960  will only be used in the correctness proof.