From a7eba53efc1b48f357119e9ca7c00b9bdce124a3 Mon Sep 17 00:00:00 2001
From: Neil Toronto <neil.toronto@gmail.com>
Date: Fri, 4 Apr 2014 11:23:20 -0600
Subject: [PATCH] Removed all flonum speed hacks; plots are usually no slower
 and are sometimes faster now

In an ultimately vain attempt to make plotting fast but still allow
plots with bounds that have very few flonums in them, Plot used to try
to detect when each axis's bounds had enough flonums in it and switch to
flonum-only math. There are two problems with this:

 * It was trying to determine whether *future computations* would have
   enough precision, which is very difficult to know for sure, so it
   only *usually* worked. (Example fail: contour-intervals3d with 0 and
   +max.0 endpoints. Half the isobands weren't drawn.)

 * It caused extra conversions between flonums and exact rationals,
   which are very slow.

Here's the new approach:

 1. Always send exact rationals to user functions (e.g. the lambda
    arguments to `function' and `contour-intervals3d').

 2. Immediately convert user-given plot coordinates to exact
    rationals (if rational; e.g. not +nan.0)

 3. Always represent normalized coordinates (e.g. in [0,1]^2 for 2D
    plots and [-1/2,1/2]^3 for 3D plots) using flonums.

 4. Reduce the number of exact operations as much as possible.

IOW, plot coordinates, which can be incredibly huge or tiny and don't
always fit in flonums, are always exact internally. Normalized, view,
and device coordinates, which are used only internally, always have
bounds with plenty of flonums in them, so Plot uses flonums.

Doing #4 accounts for most of the changes; e.g. to the marching
squares and marching cubes implementations, and axial plane clipping.

Most plots' speeds seem to be unaffected by the changes. Surfaces and
contour intervals are sometimes faster. Isosurfaces are a tad slower
in some cases and faster in others. Points are about 50% slower,
partly undoing the speedup from a recent commit. Plots with axis
transforms can be much slower (4x), when long, straight lines are
subdivided many times. Plots with bounds outside flonum range seem
to be about 3x faster.
---
 .../plot/private/common/marching-cubes.rkt    | 167 +++--
 .../plot/private/common/marching-squares.rkt  | 584 +++++++++---------
 .../plot/private/common/marching-utils.rkt    |   8 +-
 .../plot-lib/plot/private/common/math.rkt     |  23 -
 .../plot-lib/plot/private/common/sample.rkt   |  78 +--
 .../plot-lib/plot/private/contracted/math.rkt |   2 +-
 .../plot/private/contracted/sample.rkt        |   6 -
 .../plot-lib/plot/private/plot2d/clip.rkt     | 295 +++++----
 .../plot-lib/plot/private/plot2d/contour.rkt  |  23 +-
 .../plot/private/plot2d/decoration.rkt        |   7 +-
 .../plot-lib/plot/private/plot2d/interval.rkt |   2 +-
 .../plot/private/plot2d/plot-area.rkt         | 155 ++---
 .../plot-lib/plot/private/plot2d/point.rkt    |   2 +-
 .../plot-lib/plot/private/plot2d/vector.rkt   |  46 ++
 .../plot-lib/plot/private/plot3d/bsp.rkt      |   2 -
 .../plot-lib/plot/private/plot3d/clip.rkt     | 354 ++++++-----
 .../plot-lib/plot/private/plot3d/contour.rkt  |  25 +-
 .../plot/private/plot3d/isosurface.rkt        |  55 +-
 .../plot-lib/plot/private/plot3d/matrix.rkt   |  61 --
 .../plot/private/plot3d/plot-area.rkt         | 229 ++++---
 .../plot-lib/plot/private/plot3d/split.rkt    |   3 +-
 .../plot-lib/plot/private/plot3d/surface.rkt  |  12 +-
 .../plot-lib/plot/private/plot3d/vector.rkt   | 144 +++++
 .../plot/tests/extreme-bounds-tests.rkt       |   6 +-
 .../plot-test/plot/tests/plot2d-tests.rkt     |   1 +
 .../plot-test/plot/tests/plot3d-bsp-tests.rkt |  19 +-
 .../typed/parameter-regression-tests.rkt      |  12 +-
 27 files changed, 1184 insertions(+), 1137 deletions(-)
 create mode 100644 pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt
 delete mode 100644 pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt
 create mode 100644 pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt

diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt
index 15ab082e62..11b1d12457 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-cubes.rkt
@@ -17,20 +17,20 @@
 
 ;; edge vertexes 
 
-(define-syntax-rule (edge-1-2 d d1 d2) (vector (unsafe-solve-t d d1 d2) 0.0 0.0))
-(define-syntax-rule (edge-2-3 d d2 d3) (vector 1.0 (unsafe-solve-t d d2 d3) 0.0))
-(define-syntax-rule (edge-3-4 d d3 d4) (vector (unsafe-solve-t d d4 d3) 1.0 0.0))
-(define-syntax-rule (edge-1-4 d d1 d4) (vector 0.0 (unsafe-solve-t d d1 d4) 0.0))
+(define-syntax-rule (edge-1-2 d d1 d2) (vector (solve-t d d1 d2) 0 0))
+(define-syntax-rule (edge-2-3 d d2 d3) (vector 1 (solve-t d d2 d3) 0))
+(define-syntax-rule (edge-3-4 d d3 d4) (vector (solve-t d d4 d3) 1 0))
+(define-syntax-rule (edge-1-4 d d1 d4) (vector 0 (solve-t d d1 d4) 0))
 
-(define-syntax-rule (edge-5-6 d d5 d6) (vector (unsafe-solve-t d d5 d6) 0.0 1.0))
-(define-syntax-rule (edge-6-7 d d6 d7) (vector 1.0 (unsafe-solve-t d d6 d7) 1.0))
-(define-syntax-rule (edge-7-8 d d7 d8) (vector (unsafe-solve-t d d7 d8) 1.0 1.0))
-(define-syntax-rule (edge-5-8 d d5 d8) (vector 0.0 (unsafe-solve-t d d5 d8) 1.0))
+(define-syntax-rule (edge-5-6 d d5 d6) (vector (solve-t d d5 d6) 0 1))
+(define-syntax-rule (edge-6-7 d d6 d7) (vector 1 (solve-t d d6 d7) 1))
+(define-syntax-rule (edge-7-8 d d7 d8) (vector (solve-t d d7 d8) 1 1))
+(define-syntax-rule (edge-5-8 d d5 d8) (vector 0 (solve-t d d5 d8) 1))
 
-(define-syntax-rule (edge-1-5 d d1 d5) (vector 0.0 0.0 (unsafe-solve-t d d1 d5)))
-(define-syntax-rule (edge-2-6 d d2 d6) (vector 1.0 0.0 (unsafe-solve-t d d2 d6)))
-(define-syntax-rule (edge-3-7 d d3 d7) (vector 1.0 1.0 (unsafe-solve-t d d3 d7)))
-(define-syntax-rule (edge-4-8 d d4 d8) (vector 0.0 1.0 (unsafe-solve-t d d4 d8)))
+(define-syntax-rule (edge-1-5 d d1 d5) (vector 0 0 (solve-t d d1 d5)))
+(define-syntax-rule (edge-2-6 d d2 d6) (vector 1 0 (solve-t d d2 d6)))
+(define-syntax-rule (edge-3-7 d d3 d7) (vector 1 1 (solve-t d d3 d7)))
+(define-syntax-rule (edge-4-8 d d4 d8) (vector 0 1 (solve-t d d4 d8)))
 
 #|
 Cube vertex numbers:
@@ -72,7 +72,7 @@ Cube vertex numbers:
 (define known-cube-0000-1111 known-cube-1111-0000)
 
 (define ((make-known-cube-1010-0000 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d1 d2 d3 d4))
+  (define da (/ (+ d1 d2 d3 d4) 4))
   (cond
     [(test? da d)
      (list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3)
@@ -83,8 +83,8 @@ Cube vertex numbers:
      (list (list (edge-1-2 d d1 d2) (edge-1-5 d d1 d5) (edge-1-4 d d1 d4))
            (list (edge-2-3 d d2 d3) (edge-3-4 d d3 d4) (edge-3-7 d d3 d7)))]))
 
-(define known-cube-1010-0000 (make-known-cube-1010-0000 unsafe-fl>=))
-(define known-cube-0101-1111 (make-known-cube-1010-0000 unsafe-fl<))
+(define known-cube-1010-0000 (make-known-cube-1010-0000 >=))
+(define known-cube-0101-1111 (make-known-cube-1010-0000 <))
 
 (define (known-cube-1000-0010 d d1 d2 d3 d4 d5 d6 d7 d8)
   (list (list (edge-1-2 d d1 d2) (edge-1-5 d d1 d5) (edge-1-4 d d1 d4))
@@ -93,7 +93,7 @@ Cube vertex numbers:
 (define known-cube-0111-1101 known-cube-1000-0010)
 
 (define ((make-known-cube-1100-0010 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d2 d6 d7 d3))
+  (define da (/ (+ d2 d6 d7 d3) 4))
   (cond
     [(test? da d)
      (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)
@@ -106,12 +106,12 @@ Cube vertex numbers:
                  (edge-2-3 d d2 d3) (edge-1-4 d d1 d4))
            (list (edge-6-7 d d6 d7) (edge-3-7 d d3 d7) (edge-7-8 d d7 d8)))]))
 
-(define known-cube-1100-0010 (make-known-cube-1100-0010 unsafe-fl>=))
-(define known-cube-0011-1101 (make-known-cube-1100-0010 unsafe-fl<))
+(define known-cube-1100-0010 (make-known-cube-1100-0010 >=))
+(define known-cube-0011-1101 (make-known-cube-1100-0010 <))
 
 (define ((make-known-cube-1100-0011 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d1 d5 d8 d4))
-  (define db (unsafe-flavg4 d2 d6 d7 d3))
+  (define da (/ (+ d1 d5 d8 d4) 4))
+  (define db (/ (+ d2 d6 d7 d3) 4))
   (cond
     [(and (test? da d) (test? db d))
      (list (list (edge-1-5 d d1 d5) (edge-5-8 d d5 d8)
@@ -123,7 +123,7 @@ Cube vertex numbers:
                             (v+ (edge-4-8 d d4 d8) (edge-5-8 d d5 d8)))
                         (v+ (v+ (edge-6-7 d d6 d7) (edge-2-6 d d2 d6))
                             (v+ (edge-2-3 d d2 d3) (edge-3-7 d d3 d7))))
-                    0.125))
+                    1/8))
      (list (list ec (edge-5-8 d d5 d8) (edge-6-7 d d6 d7))
            (list ec (edge-1-5 d d1 d5) (edge-2-6 d d2 d6))
            (list ec (edge-4-8 d d4 d8) (edge-3-7 d d3 d7))
@@ -137,7 +137,7 @@ Cube vertex numbers:
                             (v+ (edge-4-8 d d4 d8) (edge-5-8 d d5 d8)))
                         (v+ (v+ (edge-6-7 d d6 d7) (edge-2-6 d d2 d6))
                             (v+ (edge-2-3 d d2 d3) (edge-3-7 d d3 d7))))
-                    0.125))
+                    1/8))
      (list (list ec (edge-5-8 d d5 d8) (edge-6-7 d d6 d7))
            (list ec (edge-1-5 d d1 d5) (edge-2-6 d d2 d6))
            (list ec (edge-4-8 d d4 d8) (edge-3-7 d d3 d7))
@@ -152,8 +152,8 @@ Cube vertex numbers:
            (list (edge-5-8 d d5 d8) (edge-6-7 d d6 d7)
                  (edge-3-7 d d3 d7) (edge-4-8 d d4 d8)))]))
 
-(define known-cube-1100-0011 (make-known-cube-1100-0011 unsafe-fl>=))
-(define known-cube-0011-1100 (make-known-cube-1100-0011 unsafe-fl<))
+(define known-cube-1100-0011 (make-known-cube-1100-0011 >=))
+(define known-cube-0011-1100 (make-known-cube-1100-0011 <))
 
 #|
 Cube vertex numbers:
@@ -169,12 +169,12 @@ Cube vertex numbers:
 
 
 (define ((make-known-cube-1010-0101 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d1 d2 d3 d4))
-  (define db (unsafe-flavg4 d1 d5 d8 d4))
-  (define dc (unsafe-flavg4 d3 d4 d8 d7))
-  (define dd (unsafe-flavg4 d1 d2 d6 d5))
-  (define de (unsafe-flavg4 d2 d3 d7 d6))
-  (define df (unsafe-flavg4 d5 d6 d7 d8))
+  (define da (/ (+ d1 d2 d3 d4) 4))
+  (define db (/ (+ d1 d5 d8 d4) 4))
+  (define dc (/ (+ d3 d4 d8 d7) 4))
+  (define dd (/ (+ d1 d2 d6 d5) 4))
+  (define de (/ (+ d2 d3 d7 d6) 4))
+  (define df (/ (+ d5 d6 d7 d8) 4))
   (append
    (list (list (edge-1-5 d d1 d5) (edge-1-2 d d1 d2) (edge-1-4 d d1 d4))
          (list (edge-2-3 d d2 d3) (edge-3-7 d d3 d7) (edge-3-4 d d3 d4))
@@ -205,26 +205,26 @@ Cube vertex numbers:
                    (edge-7-8 d d7 d8) (edge-5-8 d d5 d8)))
        empty)))
 
-(define known-cube-1010-0101 (make-known-cube-1010-0101 unsafe-fl>=))
-(define known-cube-0101-1010 (make-known-cube-1010-0101 unsafe-fl<))
+(define known-cube-1010-0101 (make-known-cube-1010-0101 >=))
+(define known-cube-0101-1010 (make-known-cube-1010-0101 <))
 
-(define ((make-known-cube-1110-0001 unsafe-fl>=) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d1 d5 d8 d4))
-  (define db (unsafe-flavg4 d7 d8 d4 d3))
+(define ((make-known-cube-1110-0001 >=) d d1 d2 d3 d4 d5 d6 d7 d8)
+  (define da (/ (+ d1 d5 d8 d4) 4))
+  (define db (/ (+ d7 d8 d4 d3) 4))
   (cond
-    [(and (da . unsafe-fl>= . d) (db . unsafe-fl>= . d))
+    [(and (da . >= . d) (db . >= . d))
      (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7))
            (list (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)
                  (edge-7-8 d d7 d8) (edge-5-8 d d5 d8))
            (list (edge-1-4 d d1 d4) (edge-3-4 d d3 d4) (edge-4-8 d d4 d8)))]
-    [(da . unsafe-fl>= . d)
-     (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5))
+    [(da . >= . d)
+     (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 1/2))
      (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7))
            (list ec (edge-3-7 d d3 d7) (edge-3-4 d d3 d4) (edge-1-4 d d1 d4))
            (list ec (edge-1-5 d d1 d5) (edge-5-8 d d5 d8) (edge-7-8 d d7 d8))
            (list ec (edge-1-4 d d1 d4) (edge-4-8 d d4 d8) (edge-7-8 d d7 d8)))]
-    [(db . unsafe-fl>= . d)
-     (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 0.5))
+    [(db . >= . d)
+     (define ec (v* (v+ (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)) 1/2))
      (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6) (edge-3-7 d d3 d7))
            (list ec (edge-1-5 d d1 d5) (edge-1-4 d d1 d4) (edge-3-4 d d3 d4))
            (list ec (edge-3-7 d d3 d7) (edge-7-8 d d7 d8) (edge-5-8 d d5 d8))
@@ -235,8 +235,8 @@ Cube vertex numbers:
                  (edge-3-4 d d3 d4) (edge-1-4 d d1 d4))
            (list (edge-7-8 d d7 d8) (edge-5-8 d d5 d8) (edge-4-8 d d4 d8)))]))
 
-(define known-cube-1110-0001 (make-known-cube-1110-0001 unsafe-fl>=))
-(define known-cube-0001-1110 (make-known-cube-1110-0001 unsafe-fl<))
+(define known-cube-1110-0001 (make-known-cube-1110-0001 >=))
+(define known-cube-0001-1110 (make-known-cube-1110-0001 <))
 
 (define (known-cube-1110-0100 d d1 d2 d3 d4 d5 d6 d7 d8)
   (list (list (edge-1-5 d d1 d5) (edge-5-6 d d5 d6) (edge-6-7 d d6 d7)
@@ -251,9 +251,9 @@ Cube vertex numbers:
 (define known-cube-0001-1101 known-cube-1110-0010)
 
 (define ((make-known-cube-1010-0001 test?) d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define da (unsafe-flavg4 d1 d2 d3 d4))
-  (define db (unsafe-flavg4 d1 d5 d8 d4))
-  (define dc (unsafe-flavg4 d3 d4 d8 d7))
+  (define da (/ (+ d1 d2 d3 d4) 4))
+  (define db (/ (+ d1 d5 d8 d4) 4))
+  (define dc (/ (+ d3 d4 d8 d7) 4))
   (append
    (list (list (edge-1-5 d d1 d5) (edge-1-2 d d1 d2) (edge-1-4 d d1 d4))
          (list (edge-2-3 d d2 d3) (edge-3-7 d d3 d7) (edge-3-4 d d3 d4))
@@ -271,8 +271,8 @@ Cube vertex numbers:
                    (edge-7-8 d d7 d8) (edge-4-8 d d4 d8)))
        empty)))
 
-(define known-cube-1010-0001 (make-known-cube-1010-0001 unsafe-fl>=))
-(define known-cube-0101-1110 (make-known-cube-1010-0001 unsafe-fl<))
+(define known-cube-1010-0001 (make-known-cube-1010-0001 >=))
+(define known-cube-0101-1110 (make-known-cube-1010-0001 <))
 
 (define-for-syntax known-cubes
   '((0 0 0 0 0 0 0 0)
@@ -308,7 +308,7 @@ Cube vertex numbers:
 
 (define (mirror-vec-d v)
   (match-define (vector x y d) v)
-  (vector x y (unsafe-fl- 1.0 d)))
+  (vector x y (- 1 d)))
 
 (define ((mirror-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map mirror-vec-d poly))
@@ -316,7 +316,7 @@ Cube vertex numbers:
 
 (define (mirror-vec-y v)
   (match-define (vector x y d) v)
-  (vector x (unsafe-fl- 1.0 y) d))
+  (vector x (- 1 y) d))
 
 (define ((mirror-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map mirror-vec-y poly))
@@ -324,7 +324,7 @@ Cube vertex numbers:
 
 (define (mirror-vec-x v)
   (match-define (vector x y d) v)
-  (vector (unsafe-fl- 1.0 x) y d))
+  (vector (- 1 x) y d))
 
 (define ((mirror-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map mirror-vec-x poly))
@@ -334,7 +334,7 @@ Cube vertex numbers:
 
 (define (unrotate-vec-d v)
   (match-define (vector x y d) v)
-  (vector (unsafe-fl- 1.0 y) x d))
+  (vector (- 1 y) x d))
 
 (define ((rotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map unrotate-vec-d poly))
@@ -342,7 +342,7 @@ Cube vertex numbers:
 
 (define (unrotate-vec-y v)
   (match-define (vector x y d) v)
-  (vector d y (unsafe-fl- 1.0 x)))
+  (vector d y (- 1 x)))
 
 (define ((rotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map unrotate-vec-y poly))
@@ -350,7 +350,7 @@ Cube vertex numbers:
 
 (define (unrotate-vec-x v)
   (match-define (vector x y d) v)
-  (vector x (unsafe-fl- 1.0 d) y))
+  (vector x (- 1 d) y))
 
 (define ((rotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map unrotate-vec-x poly))
@@ -360,7 +360,7 @@ Cube vertex numbers:
 
 (define (rotate-vec-d v)
   (match-define (vector x y d) v)
-  (vector y (unsafe-fl- 1.0 x) d))
+  (vector y (- 1 x) d))
 
 (define ((unrotate-cube-d f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map rotate-vec-d poly))
@@ -368,7 +368,7 @@ Cube vertex numbers:
 
 (define (rotate-vec-y v)
   (match-define (vector x y d) v)
-  (vector (unsafe-fl- 1.0 d) y x))
+  (vector (- 1 d) y x))
 
 (define ((unrotate-cube-y f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map rotate-vec-y poly))
@@ -376,7 +376,7 @@ Cube vertex numbers:
 
 (define (rotate-vec-x v)
   (match-define (vector x y d) v)
-  (vector x d (unsafe-fl- 1.0 y)))
+  (vector x d (- 1 y)))
 
 (define ((unrotate-cube-x f) d d1 d2 d3 d4 d5 d6 d7 d8)
   (map (λ (poly) (map rotate-vec-x poly))
@@ -511,15 +511,15 @@ Cube vertex numbers:
 (define-syntax-rule (add-digit j idx)
   (unsafe-fx+ (unsafe-fx* idx 2) j))
 
-(define (unsafe-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8)
-  (define j1 (if (d1 . unsafe-fl< . d) 0 1))
-  (define j2 (if (d2 . unsafe-fl< . d) 0 1))
-  (define j3 (if (d3 . unsafe-fl< . d) 0 1))
-  (define j4 (if (d4 . unsafe-fl< . d) 0 1))
-  (define j5 (if (d5 . unsafe-fl< . d) 0 1))
-  (define j6 (if (d6 . unsafe-fl< . d) 0 1))
-  (define j7 (if (d7 . unsafe-fl< . d) 0 1))
-  (define j8 (if (d8 . unsafe-fl< . d) 0 1))
+(define (do-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8)
+  (define j1 (if (d1 . < . d) 0 1))
+  (define j2 (if (d2 . < . d) 0 1))
+  (define j3 (if (d3 . < . d) 0 1))
+  (define j4 (if (d4 . < . d) 0 1))
+  (define j5 (if (d5 . < . d) 0 1))
+  (define j6 (if (d6 . < . d) 0 1))
+  (define j7 (if (d7 . < . d) 0 1))
+  (define j8 (if (d8 . < . d) 0 1))
   (define facet-num
     (add-digit
      j1 (add-digit
@@ -536,37 +536,12 @@ Cube vertex numbers:
                               [d1 real?] [d2 real?] [d3 real?] [d4 real?]
                               [d5 real?] [d6 real?] [d7 real?] [d8 real?]
                               ) (listof (listof (vector/c real? real? real?)))
-  (cond [(all inexact-real? xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)
-         (define polys (unsafe-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8))
+  (cond [(= d d1 d2 d3 d4 d5 d6 d7 d8)  empty]
+        [else
+         (define polys (do-heights->cube-polys d d1 d2 d3 d4 d5 d6 d7 d8))
          (for/list ([poly  (in-list polys)])
            (for/list ([uvw  (in-list poly)])
              (match-define (vector u v w) uvw)
-             (vector (unsafe-unsolve-t xa xb u)
-                     (unsafe-unsolve-t ya yb v)
-                     (unsafe-unsolve-t za zb w))))]
-        [(find-failure-index real? xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)
-         => (λ (i) (raise-type-error 'heights->polys "real number"
-                                     i xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))]
-        [(= d d1 d2 d3 d4 d5 d6 d7 d8)  empty]
-        [else
-         (let-map
-          (d d1 d2 d3 d4 d5 d6 d7 d8) inexact->exact
-          (define d-min (min d d1 d2 d3 d4 d5 d6 d7 d8))
-          (define d-max (max d d1 d2 d3 d4 d5 d6 d7 d8))
-          (define d-scale (- d-max d-min))
-          (define polys
-            (unsafe-heights->cube-polys (exact->inexact (/ (- d d-min) d-scale))
-                                        (exact->inexact (/ (- d1 d-min) d-scale))
-                                        (exact->inexact (/ (- d2 d-min) d-scale))
-                                        (exact->inexact (/ (- d3 d-min) d-scale))
-                                        (exact->inexact (/ (- d4 d-min) d-scale))
-                                        (exact->inexact (/ (- d5 d-min) d-scale))
-                                        (exact->inexact (/ (- d6 d-min) d-scale))
-                                        (exact->inexact (/ (- d7 d-min) d-scale))
-                                        (exact->inexact (/ (- d8 d-min) d-scale))))
-          (for/list ([poly  (in-list polys)])
-            (for/list ([uvw  (in-list poly)])
-              (match-define (vector u v w) uvw)
-              (vector (unsolve-t xa xb u)
-                      (unsolve-t ya yb v)
-                      (unsolve-t za zb w)))))]))
+             (vector (unsolve-t xa xb u)
+                     (unsolve-t ya yb v)
+                     (unsolve-t za zb w))))]))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt
index 4bed76be26..09598c2826 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-squares.rkt
@@ -2,6 +2,7 @@
 
 (require racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops racket/contract
          unstable/latent-contract/defthing
+         (only-in math/flonum fl)
          (for-syntax racket/base)
          "math.rkt"
          "utils.rkt"
@@ -29,18 +30,6 @@ first argument is the contour value; the rest are z coordinates arranged as
 above.
 |#
 
-(define-syntax-rule (unrotate-vec v)
-  (match-let ([(vector x y z)  v])
-    (vector (unsafe-fl- 1.0 y) x z)))
-
-(define-syntax-rule (mirror-x-vec v)
-  (match-let ([(vector x y z)  v])
-    (vector (fl- 1.0 x) y z)))
-
-(define-syntax-rule (mirror-y-vec v)
-  (match-let ([(vector x y z)  v])
-    (vector x (fl- 1.0 y) z)))
-
 ;; =============================================================================
 ;; Contour lines
 
@@ -59,23 +48,23 @@ above.
 
 ;(: lines1000 (FType Lines))
 (define-syntax-rule (lines1000 z z1 z2 z3 z4)
-  (list (vector (unsafe-solve-t z z1 z2) 0.0
-                0.0 (unsafe-solve-t z z1 z4))))
+  (list (vector (solve-t z z1 z2) 0
+                0 (solve-t z z1 z4))))
 
 ;(: lines0100 (FType Lines))
 (define-syntax-rule (lines0100 z z1 z2 z3 z4)
-  (list (vector (unsafe-solve-t z z1 z2) 0.0
-                1.0 (unsafe-solve-t z z2 z3))))
+  (list (vector (solve-t z z1 z2) 0
+                1 (solve-t z z2 z3))))
 
 ;(: lines0010 (FType Lines))
 (define-syntax-rule (lines0010 z z1 z2 z3 z4)
-  (list (vector 1.0 (unsafe-solve-t z z2 z3)
-                (unsafe-solve-t z z4 z3) 1.0)))
+  (list (vector 1 (solve-t z z2 z3)
+                (solve-t z z4 z3) 1)))
 
 ;(: lines0001 (FType Lines))
 (define-syntax-rule (lines0001 z z1 z2 z3 z4)
-  (list (vector 0.0 (unsafe-solve-t z z1 z4)
-                (unsafe-solve-t z z4 z3) 1.0)))
+  (list (vector 0 (solve-t z z1 z4)
+                (solve-t z z4 z3) 1)))
 
 (define-syntax-rule (lines0111 z z1 z2 z3 z4) (lines1000 z z1 z2 z3 z4))
 (define-syntax-rule (lines1011 z z1 z2 z3 z4) (lines0100 z z1 z2 z3 z4))
@@ -87,13 +76,13 @@ above.
 
 ;(: lines1100 (FType Lines))
 (define-syntax-rule (lines1100 z z1 z2 z3 z4)
-  (list (vector 0.0 (unsafe-solve-t z z1 z4)
-                1.0 (unsafe-solve-t z z2 z3))))
+  (list (vector 0 (solve-t z z1 z4)
+                1 (solve-t z z2 z3))))
 
 ;(: lines0110 (FType Lines))
 (define-syntax-rule (lines0110 z z1 z2 z3 z4)
-  (list (vector (unsafe-solve-t z z1 z2) 0.0
-                (unsafe-solve-t z z4 z3) 1.0)))
+  (list (vector (solve-t z z1 z2) 0
+                (solve-t z z4 z3) 1)))
 
 (define-syntax-rule (lines0011 z z1 z2 z3 z4) (lines1100 z z1 z2 z3 z4))
 (define-syntax-rule (lines1001 z z1 z2 z3 z4) (lines0110 z z1 z2 z3 z4))
@@ -104,25 +93,27 @@ above.
 ;(: lines-opposite ((Float Float -> Boolean) -> (FType Lines)))
 (define-syntax-rule (lines-opposite test? z z1 z2 z3 z4)
   ; disambiguate using average of corners as guess for center value
-  (let ([z5  (unsafe-flavg4 z1 z2 z3 z4)])
+  (let ([z5  (/ (+ z1 z2 z3 z4) 4)])
     (if (test? z5 z)
-        (list (vector (unsafe-solve-t z z1 z2) 0.0
-                      1.0 (unsafe-solve-t z z2 z3))
-              (vector 0.0 (unsafe-solve-t z z1 z4)
-                      (unsafe-solve-t z z4 z3) 1.0))
-        (list (vector (unsafe-solve-t z z1 z2) 0.0
-                      0.0 (unsafe-solve-t z z1 z4))
-              (vector 1.0 (unsafe-solve-t z z2 z3)
-                      (unsafe-solve-t z z4 z3) 1.0)))))
+        (list (vector (solve-t z z1 z2) 0
+                      1 (solve-t z z2 z3))
+              (vector 0 (solve-t z z1 z4)
+                      (solve-t z z4 z3) 1))
+        (list (vector (solve-t z z1 z2) 0
+                      0 (solve-t z z1 z4))
+              (vector 1 (solve-t z z2 z3)
+                      (solve-t z z4 z3) 1)))))
 
-(define-syntax-rule (lines1010 z z1 z2 z3 z4) (lines-opposite unsafe-fl>= z z1 z2 z3 z4))
-(define-syntax-rule (lines0101 z z1 z2 z3 z4) (lines-opposite unsafe-fl< z z1 z2 z3 z4))
+(define-syntax-rule (lines1010 z z1 z2 z3 z4) (lines-opposite >= z z1 z2 z3 z4))
+(define-syntax-rule (lines0101 z z1 z2 z3 z4) (lines-opposite <  z z1 z2 z3 z4))
 
-(define (unsafe-heights->lines z z1 z2 z3 z4)
-  (define p1 (z1 . unsafe-fl>= . z))
-  (define p2 (z2 . unsafe-fl>= . z))
-  (define p3 (z3 . unsafe-fl>= . z))
-  (define p4 (z4 . unsafe-fl>= . z))
+(define (do-heights->lines z z1 z2 z3 z4)
+  ;(printf "~v ~v ~v ~v ~v~n" z z1 z2 z3 z4)
+  (define p1 (z1 . >= . z))
+  (define p2 (z2 . >= . z))
+  (define p3 (z3 . >= . z))
+  (define p4 (z4 . >= . z))
+  ;(printf "~v ~v ~v ~v~n~n" p1 p2 p3 p4)
   (if p1
       (if p2
           (if p3
@@ -158,46 +149,45 @@ above.
 (defproc (heights->lines [xa real?] [xb real?] [ya real?] [yb real?]
                          [z real?] [z1 real?] [z2 real?] [z3 real?] [z4 real?]
                          ) (listof (list/c (vector/c real? real? real?) (vector/c real? real? real?)))
-  (cond [(all inexact-real? xa xb ya yb z z1 z2 z3 z4)
-         (define lines (unsafe-heights->lines z z1 z2 z3 z4))
-         (for/list ([line  (in-list lines)])
-           (match-define (vector u1 v1 u2 v2) line)
-           (list (vector (unsafe-unsolve-t xa xb u1) (unsafe-unsolve-t ya yb v1) z)
-                 (vector (unsafe-unsolve-t xa xb u2) (unsafe-unsolve-t ya yb v2) z)))]
-        [(find-failure-index real? xa xb ya yb z z1 z2 z3 z4)
-         => (λ (i) (raise-type-error 'heights->liens "real number" i xa xb ya yb z z1 z2 z3 z4))]
+  (cond [(and (z . < . z1) (z . < . z2) (z . < . z3) (z . < . z4))  empty]
+        [(and (z . > . z1) (z . > . z2) (z . > . z3) (z . > . z4))  empty]
         [(= z z1 z2 z3 z4)  empty]
         [else
-         (let-map
-          (z z1 z2 z3 z4) inexact->exact
-          (define z-min (min z z1 z2 z3 z4))
-          (define z-max (max z z1 z2 z3 z4))
-          (define z-scale (- z-max z-min))
-          (define lines
-            (unsafe-heights->lines (exact->inexact (/ (- z z-min) z-scale))
-                                   (exact->inexact (/ (- z1 z-min) z-scale))
-                                   (exact->inexact (/ (- z2 z-min) z-scale))
-                                   (exact->inexact (/ (- z3 z-min) z-scale))
-                                   (exact->inexact (/ (- z4 z-min) z-scale))))
-          (for/list ([line  (in-list lines)])
-            (match-define (vector u1 v1 u2 v2) line)
-            (list (vector (unsolve-t xa xb u1) (unsolve-t ya yb v1) z)
-                  (vector (unsolve-t xa xb u2) (unsolve-t ya yb v2) z))))]))
+         (define lines (do-heights->lines z z1 z2 z3 z4))
+         (for/list ([line  (in-list lines)])
+           (match-define (vector u1 v1 u2 v2) line)
+           (list (vector (unsolve-t xa xb u1) (unsolve-t ya yb v1) z)
+                 (vector (unsolve-t xa xb u2) (unsolve-t ya yb v2) z)))]))
 
 ;; =============================================================================
 ;; Isoband marching squares: polygonizes contour between two isoline values
 
-(define ((rotate-facet f) za zb z1 z2 z3 z4)
-  (map (λ (poly) (map (λ (v) (unrotate-vec v)) poly))
-       (f za zb z2 z3 z4 z1)))
+(define-syntax-rule (unrotate-vec v)
+  (match-let ([(vector x y z)  v])
+    (vector (- 1 y) x z)))
 
-(define ((mirror-x-facet f) za zb z1 z2 z3 z4)
-  (map (λ (poly) (map (λ (v) (mirror-x-vec v)) poly))
-       (f za zb z2 z1 z4 z3)))
+(define-syntax-rule (unrotate-vec2 v)
+  (match-let ([(vector x y z)  v])
+    (vector (- 1 x) (- 1 y) z)))
 
-(define ((mirror-y-facet f) za zb z1 z2 z3 z4)
-  (map (λ (poly) (map (λ (v) (mirror-y-vec v)) poly))
-       (f za zb z4 z3 z2 z1)))
+(define-syntax-rule (unrotate-vec3 v)
+  (match-let ([(vector x y z)  v])
+    (vector y (- 1 x) z)))
+
+(define-syntax-rule (rotate-facet f)
+  (λ (za zb z1 z2 z3 z4)
+    (map (λ (vs) (map (λ (v) (unrotate-vec v)) vs))
+         (f za zb z2 z3 z4 z1))))
+
+(define-syntax-rule (rotate-facet2 f)
+  (λ (za zb z1 z2 z3 z4)
+    (map (λ (vs) (map (λ (v) (unrotate-vec2 v)) vs))
+         (f za zb z3 z4 z1 z2))))
+
+(define-syntax-rule (rotate-facet3 f)
+  (λ (za zb z1 z2 z3 z4)
+    (map (λ (vs) (map (λ (v) (unrotate-vec3 v)) vs))
+         (f za zb z4 z1 z2 z3))))
 
 ;; -----------------------------------------------------------------------------
 ;; all corners same
@@ -210,233 +200,233 @@ above.
 ;; single triangle
 
 (define (polys1000 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys0100 (rotate-facet polys1000))
-(define polys0010 (rotate-facet polys0100))
-(define polys0001 (rotate-facet polys0010))
+(define polys0010 (rotate-facet2 polys1000))
+(define polys0001 (rotate-facet3 polys1000))
 
 (define (polys1222 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys2122 (rotate-facet polys1222))
-(define polys2212 (rotate-facet polys2122))
-(define polys2221 (rotate-facet polys2212))
+(define polys2212 (rotate-facet2 polys1222))
+(define polys2221 (rotate-facet3 polys1222))
 
 ;; -----------------------------------------------------------------------------
 ;; single trapezoid
 
 (define (polys2000 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector (solve-t zb z1 z2) 0 zb)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 0 (solve-t za z1 z4) za)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys0200 (rotate-facet polys2000))
-(define polys0020 (rotate-facet polys0200))
-(define polys0002 (rotate-facet polys0020))
+(define polys0020 (rotate-facet2 polys2000))
+(define polys0002 (rotate-facet3 polys2000))
 
 (define (polys0222 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys2022 (rotate-facet polys0222))
-(define polys2202 (rotate-facet polys2022))
-(define polys2220 (rotate-facet polys2202))
+(define polys2202 (rotate-facet2 polys0222))
+(define polys2220 (rotate-facet3 polys0222))
 
 ;; -----------------------------------------------------------------------------
 ;; single rectangle
 
 (define (polys1100 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector 0 0 z1)
+              (vector 1 0 z2)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys0110 (rotate-facet polys1100))
-(define polys0011 (rotate-facet polys0110))
-(define polys1001 (rotate-facet polys0011))
+(define polys0011 (rotate-facet2 polys1100))
+(define polys1001 (rotate-facet3 polys1100))
 
 (define (polys1122 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 0 z1)
+              (vector 1 0 z2)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys2112 (rotate-facet polys1122))
-(define polys2211 (rotate-facet polys2112))
-(define polys1221 (rotate-facet polys2211))
+(define polys2211 (rotate-facet2 polys1122))
+(define polys1221 (rotate-facet3 polys1122))
 
 (define (polys0022 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 (solve-t za z1 z4) za)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys2002 (rotate-facet polys0022))
-(define polys2200 (rotate-facet polys2002))
-(define polys0220 (rotate-facet polys2200))
+(define polys2200 (rotate-facet2 polys0022))
+(define polys0220 (rotate-facet3 polys0022))
 
 ;; -----------------------------------------------------------------------------
 ;; single pentagon
 
 (define (polys0111 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 1.0 z3)
-              (vector 0.0 1.0 z4)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector 1 0 z2)
+              (vector 1 1 z3)
+              (vector 0 1 z4)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys1011 (rotate-facet polys0111))
-(define polys1101 (rotate-facet polys1011))
-(define polys1110 (rotate-facet polys1101))
+(define polys1101 (rotate-facet2 polys0111))
+(define polys1110 (rotate-facet3 polys0111))
 
 (define (polys2111 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 1.0 z3)
-              (vector 0.0 1.0 z4)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 0 z2)
+              (vector 1 1 z3)
+              (vector 0 1 z4)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys1211 (rotate-facet polys2111))
-(define polys1121 (rotate-facet polys1211))
-(define polys1112 (rotate-facet polys1121))
+(define polys1121 (rotate-facet2 polys2111))
+(define polys1112 (rotate-facet3 polys2111))
 
 (define (polys1002 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys2100 (rotate-facet polys1002))
-(define polys0210 (rotate-facet polys2100))
-(define polys0021 (rotate-facet polys0210))
+(define polys0210 (rotate-facet2 polys1002))
+(define polys0021 (rotate-facet3 polys1002))
 
 (define (polys1220 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys0122 (rotate-facet polys1220))
-(define polys2012 (rotate-facet polys0122))
-(define polys2201 (rotate-facet polys2012))
+(define polys2012 (rotate-facet2 polys1220))
+(define polys2201 (rotate-facet3 polys1220))
 
 (define (polys1200 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys0120 (rotate-facet polys1200))
-(define polys0012 (rotate-facet polys0120))
-(define polys2001 (rotate-facet polys0012))
+(define polys0012 (rotate-facet2 polys1200))
+(define polys2001 (rotate-facet3 polys1200))
 
 (define (polys1022 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys2102 (rotate-facet polys1022))
-(define polys2210 (rotate-facet polys2102))
-(define polys0221 (rotate-facet polys2210))
+(define polys2210 (rotate-facet2 polys1022))
+(define polys0221 (rotate-facet3 polys1022))
 
 ;; -----------------------------------------------------------------------------
 ;; single hexagon
 
 (define (polys0112 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector 1 0 z2)
+              (vector 1 1 z3)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys2011 (rotate-facet polys0112))
-(define polys1201 (rotate-facet polys2011))
-(define polys1120 (rotate-facet polys1201))
+(define polys1201 (rotate-facet2 polys0112))
+(define polys1120 (rotate-facet3 polys0112))
 
 (define (polys2110 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 0 z2)
+              (vector 1 1 z3)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector 0 (solve-t za z1 z4) za)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define polys0211 (rotate-facet polys2110))
-(define polys1021 (rotate-facet polys0211))
-(define polys1102 (rotate-facet polys1021))
+(define polys1021 (rotate-facet2 polys2110))
+(define polys1102 (rotate-facet3 polys2110))
 
 (define (polys0121 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 0.0 z2)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 1.0 z4)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector 1 0 z2)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 1 z4)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define polys1012 (rotate-facet polys0121))
-(define polys2101 (rotate-facet polys1012))
-(define polys1210 (rotate-facet polys2101))
+(define polys2101 (rotate-facet2 polys0121))
+(define polys1210 (rotate-facet3 polys0121))
 
 ;; -----------------------------------------------------------------------------
 ;; 6-sided saddle
 
 (define (polys10100 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))
-        (list (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 0 (solve-t za z1 z4) za))
+        (list (vector 1 1 z3)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector 1 (solve-t za z2 z3) za))))
 
 (define (polys10101 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 1 1 z3)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define (polys1010 za zb z1 z2 z3 z4)
-  (define z5 (unsafe-flavg4 z1 z2 z3 z4))
-  (cond [(z5 . unsafe-fl< . za)  (polys10100 za zb z1 z2 z3 z4)]
+  (define z5 (/ (+ z1 z2 z3 z4) 4))
+  (cond [(z5 . < . za)  (polys10100 za zb z1 z2 z3 z4)]
         ; (z5 . >= . zb) is impossible
         [else  (polys10101 za zb z1 z2 z3 z4)]))
 
 (define polys0101 (rotate-facet polys1010))
 
 (define (polys1212-2 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))
-        (list (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 0 (solve-t zb z1 z4) zb))
+        (list (vector 1 1 z3)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 1 (solve-t zb z2 z3) zb))))
 
 (define (polys1212-1 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 0.0 z1)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector 0 0 z1)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 1 z3)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define (polys1212 za zb z1 z2 z3 z4)
-  (define z5 (unsafe-flavg4 z1 z2 z3 z4))
-  (cond [(z5 . unsafe-fl>= . zb)  (polys1212-2 za zb z1 z2 z3 z4)]
+  (define z5 (/ (+ z1 z2 z3 z4) 4))
+  (cond [(z5 . >= . zb)  (polys1212-2 za zb z1 z2 z3 z4)]
         ; (z5 . < . za) is impossible
         [else  (polys1212-1 za zb z1 z2 z3 z4)]))
 
@@ -446,98 +436,98 @@ above.
 ;; 7-sided saddle
 
 (define (polys0212-1 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 1 z3)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define (polys0212-2 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))
-        (list (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))
+        (list (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 1 z3)
+              (vector (solve-t zb z4 z3) 1 zb))))
 
 (define (polys0212 za zb z1 z2 z3 z4)
-  (define z5 (unsafe-flavg4 z1 z2 z3 z4))
-  (cond [(z5 . unsafe-fl< . zb)  (polys0212-1 za zb z1 z2 z3 z4)]
+  (define z5 (/ (+ z1 z2 z3 z4) 4))
+  (cond [(z5 . < . zb)  (polys0212-1 za zb z1 z2 z3 z4)]
         ; handling (z5 . < . za) separately results in a non-convex polygon
         [else  (polys0212-2 za zb z1 z2 z3 z4)]))
 
 (define polys2021 (rotate-facet polys0212))
-(define polys1202 (rotate-facet polys2021))
-(define polys2120 (rotate-facet polys1202))
+(define polys1202 (rotate-facet2 polys0212))
+(define polys2120 (rotate-facet3 polys0212))
 
 (define (polys2010-1 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector (solve-t zb z1 z2) 0 zb)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector 1 1 z3)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector 0 (solve-t za z1 z4) za)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define (polys2010-0 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))
-        (list (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector 1.0 1.0 z3)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za))))
+  (list (list (vector (solve-t zb z1 z2) 0 zb)
+              (vector (solve-t za z1 z2) 0 za)
+              (vector 0 (solve-t za z1 z4) za)
+              (vector 0 (solve-t zb z1 z4) zb))
+        (list (vector 1 (solve-t za z2 z3) za)
+              (vector 1 1 z3)
+              (vector (solve-t za z4 z3) 1 za))))
 
 (define (polys2010 za zb z1 z2 z3 z4)
-  (define z5 (unsafe-flavg4  z1 z2 z3 z4))
-  (cond [(z5 . unsafe-fl>= . za)  (polys2010-1 za zb z1 z2 z3 z4)]
+  (define z5 (/ (+ z1 z2 z3 z4) 4))
+  (cond [(z5 . >= . za)  (polys2010-1 za zb z1 z2 z3 z4)]
         ; handling (z5 . >= . zb) separately results in a non-convex polygon
         [else  (polys2010-0 za zb z1 z2 z3 z4)]))
 
 (define polys0201 (rotate-facet polys2010))
-(define polys1020 (rotate-facet polys0201))
-(define polys0102 (rotate-facet polys1020))
+(define polys1020 (rotate-facet2 polys2010))
+(define polys0102 (rotate-facet3 polys2010))
 
 ;; -----------------------------------------------------------------------------
 ;; 8-sided saddle
 
 (define (polys0202-0 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za))
-        (list (vector 0.0 (unsafe-solve-t za z1 z4) za)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 (solve-t za z2 z3) za))
+        (list (vector 0 (solve-t za z1 z4) za)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb))))
 
 (define (polys0202-1 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector (solve-t zb z4 z3) 1 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))))
 
 (define (polys0202-2 za zb z1 z2 z3 z4)
-  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
-              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (unsafe-solve-t zb z1 z4) zb)
-              (vector 0.0 (unsafe-solve-t za z1 z4) za))
-        (list (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
-              (vector 1.0 (unsafe-solve-t za z2 z3) za)
-              (vector (unsafe-solve-t za z4 z3) 1.0 za)
-              (vector (unsafe-solve-t zb z4 z3) 1.0 zb))))
+  (list (list (vector (solve-t za z1 z2) 0 za)
+              (vector (solve-t zb z1 z2) 0 zb)
+              (vector 0 (solve-t zb z1 z4) zb)
+              (vector 0 (solve-t za z1 z4) za))
+        (list (vector 1 (solve-t zb z2 z3) zb)
+              (vector 1 (solve-t za z2 z3) za)
+              (vector (solve-t za z4 z3) 1 za)
+              (vector (solve-t zb z4 z3) 1 zb))))
 
 (define (polys0202 za zb z1 z2 z3 z4)
-  (define z5 (unsafe-flavg4 z1 z2 z3 z4))
-  (cond [(z5 . unsafe-fl< . za)  (polys0202-0 za zb z1 z2 z3 z4)]
-        [(z5 . unsafe-fl< . zb)  (polys0202-1 za zb z1 z2 z3 z4)]
+  (define z5 (/ (+ z1 z2 z3 z4) 4))
+  (cond [(z5 . < . za)  (polys0202-0 za zb z1 z2 z3 z4)]
+        [(z5 . < . zb)  (polys0202-1 za zb z1 z2 z3 z4)]
         [else  (polys0202-2 za zb z1 z2 z3 z4)]))
 
 (define polys2020 (rotate-facet polys0202))
@@ -583,11 +573,11 @@ above.
           polys2210 polys2211 polys2212 
           polys2220 polys2221 polys2222))
 
-(define (unsafe-heights->polys za zb z1 z2 z3 z4)
-  (define t1 (if (z1 . unsafe-fl< . za) 0 (if (z1 . unsafe-fl<= . zb) 1 2)))
-  (define t2 (if (z2 . unsafe-fl< . za) 0 (if (z2 . unsafe-fl<= . zb) 1 2)))
-  (define t3 (if (z3 . unsafe-fl< . za) 0 (if (z3 . unsafe-fl<= . zb) 1 2)))
-  (define t4 (if (z4 . unsafe-fl< . za) 0 (if (z4 . unsafe-fl<= . zb) 1 2)))
+(define (do-heights->polys za zb z1 z2 z3 z4)
+  (define t1 (if (z1 . < . za) 0 (if (z1 . <= . zb) 1 2)))
+  (define t2 (if (z2 . < . za) 0 (if (z2 . <= . zb) 1 2)))
+  (define t3 (if (z3 . < . za) 0 (if (z3 . <= . zb) 1 2)))
+  (define t4 (if (z4 . < . za) 0 (if (z4 . <= . zb) 1 2)))
   (define facet-num
     (unsafe-fx+ (unsafe-fx+ (unsafe-fx+ (unsafe-fx* (unsafe-fx* (unsafe-fx* t1 3) 3) 3)
                                         (unsafe-fx* (unsafe-fx* t2 3) 3))
@@ -600,36 +590,16 @@ above.
                          [za real?] [zb real?]
                          [z1 real?] [z2 real?] [z3 real?] [z4 real?]
                          ) (listof (listof (vector/c real? real? real?)))
-  (cond [(all inexact-real? xa xb ya yb za zb z1 z2 z3 z4)
-         (define polys (unsafe-heights->polys za zb z1 z2 z3 z4))
-         (for/list ([poly  (in-list polys)])
-           (cond [(eq? poly 'full)  (list (vector xa ya z1) (vector xb ya z2)
-                                          (vector xb yb z3) (vector xa yb z4))]
-                 [else  (for/list ([uv  (in-list poly)])
-                          (match-define (vector u v z) uv)
-                          (vector (unsafe-unsolve-t xa xb u) (unsafe-unsolve-t ya yb v) z))]))]
-        [(find-failure-index real? xa xb ya yb za zb z1 z2 z3 z4)
-         => (λ (i) (raise-type-error 'heights->polys "real number" i xa xb ya yb za zb z1 z2 z3 z4))]
-        [(= za zb z1 z2 z3 z4)  (list (list (vector xa ya z1) (vector xb ya z2)
-                                            (vector xb yb z3) (vector xa yb z4)))]
-        [else
-         (let-map
-          (za zb z1 z2 z3 z4) inexact->exact
-          (define z-min (min za zb z1 z2 z3 z4))
-          (define z-max (max za zb z1 z2 z3 z4))
-          (define z-scale (- z-max z-min))
-          (define polys
-            (unsafe-heights->polys (exact->inexact (/ (- za z-min) z-scale))
-                                   (exact->inexact (/ (- zb z-min) z-scale))
-                                   (exact->inexact (/ (- z1 z-min) z-scale))
-                                   (exact->inexact (/ (- z2 z-min) z-scale))
-                                   (exact->inexact (/ (- z3 z-min) z-scale))
-                                   (exact->inexact (/ (- z4 z-min) z-scale))))
-          (for/list ([poly  (in-list polys)])
-            (cond [(eq? poly 'full)  (list (vector xa ya z1) (vector xb ya z2)
-                                           (vector xb yb z3) (vector xa yb z4))]
-                  [else  (for/list ([uv  (in-list poly)])
-                           (match-define (vector u v z) uv)
-                           (vector (unsolve-t xa xb u)
-                                   (unsolve-t ya yb v)
-                                   (+ z-min (* (inexact->exact z) z-scale))))])))]))
+  (cond
+    [(and (zb . < . z1) (zb . < . z2) (zb . < . z3) (zb . < . z4))  empty]
+    [(and (za . > . z1) (za . > . z2) (za . > . z3) (za . > . z4))  empty]
+    [(= za zb z1 z2 z3 z4)  (list (list (vector xa ya z1) (vector xb ya z2)
+                                        (vector xb yb z3) (vector xa yb z4)))]
+    [else
+     (define polys (do-heights->polys za zb z1 z2 z3 z4))
+     (for/list ([poly  (in-list polys)])
+       (cond [(eq? poly 'full)  (list (vector xa ya z1) (vector xb ya z2)
+                                      (vector xb yb z3) (vector xa yb z4))]
+             [else  (for/list ([uvz  (in-list poly)])
+                      (match-define (vector u v z) uvz)
+                      (vector (unsolve-t xa xb u) (unsolve-t ya yb v) z))]))]))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt
index 6aeace039a..0b02734e9e 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/marching-utils.rkt
@@ -10,6 +10,9 @@
 ;;           if z = zb, this returns 1.0
 ;;           if z = (za + zb) / 2, this returns 0.5
 ;; Intuitively, regard a use (solve-t z za zb) as "the point between za and zb".
+(define-syntax-rule (solve-t z za zb)
+  (/ (- z za) (- zb za)))
+
 (define-syntax-rule (unsafe-solve-t z za zb)
   (unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za)))
 
@@ -17,8 +20,9 @@
   (unsafe-fl+ (unsafe-fl* t zb) (unsafe-fl* (unsafe-fl- 1.0 t) za)))
 
 (define-syntax-rule (unsolve-t za zb t)
-  (let ([t  (inexact->exact t)])
-    (+ (* t zb) (* (- 1 t) za))))
+  (cond [(eq? t 0)  za]
+        [(eq? t 1)  zb]
+        [else  (+ (* t zb) (* (- 1 t) za))]))
 
 (define-syntax-rule (unsafe-flavg4 z1 z2 z3 z4)
   (unsafe-fl* 0.25 (unsafe-fl+ (unsafe-fl+ (unsafe-fl+ z1 z2) z3) z4)))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt
index b1ccb7eef9..1ef5248931 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/math.rkt
@@ -42,29 +42,6 @@
     [xs  (cond [(not (andmap flonum? xs))  (raise-type-error 'fldistance "flonums" xs)]
                [else  (unsafe-flsqrt (flsum (λ (x) (unsafe-fl* x x)) xs))])]))
 
-(defproc (flonum-ok-for-range? [x-min rational?] [x-max rational?]
-                               [size exact-positive-integer?]) boolean?
-  (let/ec return
-    (let ([x-min  (inexact->exact (min x-min x-max))]
-          [x-max  (inexact->exact (max x-min x-max))])
-      (define step-size (/ (- x-max x-min) size))
-      
-      (define inexact-x-min (exact->inexact x-min))
-      (unless (rational? inexact-x-min) (return #f))
-      
-      (define inexact-x-max (exact->inexact x-max))
-      (unless (rational? inexact-x-max) (return #f))
-      
-      (define inexact-x-max-prev (flprev inexact-x-max))
-      (unless (rational? inexact-x-max-prev) (return #f))
-      
-      (define inexact-x-min-next (flnext inexact-x-min))
-      (unless (rational? inexact-x-min-next) (return #f))
-      
-      (define max-diff (- x-max (inexact->exact inexact-x-max-prev)))
-      (define min-diff (- (inexact->exact inexact-x-min-next) x-min))
-      (and (max-diff . < . step-size) (min-diff . < . step-size)))))
-
 ;; ===================================================================================================
 ;; Reals
 
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt
index d6d406cf2e..29771dac12 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/common/sample.rkt
@@ -142,11 +142,15 @@
       (hash-ref! memo (vector outer-ivl num tx)
                  (λ ()
                    (define xs (sample-points outer-ivl inner-ivl num tx))
-                   (define ys (map* f xs))
-                   (define rys (filter rational? ys))
-                   (define-values (y-min y-max)
-                     (cond [(empty? rys)  (values #f #f)]
-                           [else  (values (apply min* rys) (apply max* rys))]))
+                   (define y-min #f)
+                   (define y-max #f)
+                   (define ys (for/list ([x  (in-list xs)])
+                                (define y (f x))
+                                (cond [(rational? y)
+                                       (unless (and y-min (y . >= . y-min)) (set! y-min y))
+                                       (unless (and y-max (y . <= . y-max)) (set! y-max y))
+                                       (inexact->exact y)]
+                                      [else  y])))
                    (sample xs ys
                            (maybe-inexact->exact y-min)
                            (maybe-inexact->exact y-max)))))))
@@ -173,11 +177,12 @@
                    (define z-max #f)
                    (define zss (for/vector #:length (length ys) ([y  (in-list ys)])
                                  (for/vector #:length (length xs) ([x  (in-list xs)])
-                                   (let ([z  (f x y)])
-                                     (when (rational? z)
-                                       (unless (and z-min (z . >= . z-min)) (set! z-min z))
-                                       (unless (and z-max (z . <= . z-max)) (set! z-max z)))
-                                     z))))
+                                   (define z (f x y))
+                                   (cond [(rational? z)
+                                          (unless (and z-min (z . >= . z-min)) (set! z-min z))
+                                          (unless (and z-max (z . <= . z-max)) (set! z-max z))
+                                          (inexact->exact z)]
+                                         [else  z]))))
                    (2d-sample xs ys zss
                               (maybe-inexact->exact z-min)
                               (maybe-inexact->exact z-max)))))))
@@ -208,11 +213,12 @@
                    (define dsss (for/vector #:length (length zs) ([z  (in-list zs)])
                                   (for/vector #:length (length ys) ([y  (in-list ys)])
                                     (for/vector #:length (length xs) ([x  (in-list xs)])
-                                      (let ([d  (f x y z)])
-                                        (when (rational? d)
-                                          (unless (and d-min (d . >= . d-min)) (set! d-min d))
-                                          (unless (and d-max (d . <= . d-max)) (set! d-max d)))
-                                        d)))))
+                                      (define d (f x y z))
+                                      (cond [(rational? d)
+                                             (unless (and d-min (d . >= . d-min)) (set! d-min d))
+                                             (unless (and d-max (d . <= . d-max)) (set! d-max d))
+                                             (inexact->exact d)]
+                                            [else  d])))))
                    (3d-sample xs ys zs dsss
                               (maybe-inexact->exact d-min)
                               (maybe-inexact->exact d-max)))))))
@@ -262,45 +268,3 @@
         (values yb ds01 ds11))
       (values zb dss1))))
 
-(defproc (sample-exact->inexact [s sample?]) sample?
-  (match-define (sample xs ys y-min y-max) s)
-  (sample (map exact->inexact xs) (map exact->inexact ys)
-          y-min y-max))
-
-(defproc (2d-sample-exact->inexact [s 2d-sample?]) 2d-sample?
-  (match-define (2d-sample xs ys zss z-min z-max) s)
-  (2d-sample (map exact->inexact xs) (map exact->inexact ys)
-             (for/vector #:length (vector-length zss) ([zs  (in-vector zss)])
-               (for/vector #:length (vector-length zs) ([z  (in-vector zs)])
-                 (exact->inexact z)))
-             z-min z-max))
-
-(defproc (3d-sample-exact->inexact [s 3d-sample?]) 3d-sample?
-  (match-define (3d-sample xs ys zs dsss d-min d-max) s)
-  (3d-sample (map exact->inexact xs) (map exact->inexact ys) (map exact->inexact zs)
-             (for/vector #:length (vector-length dsss) ([dss  (in-vector dsss)])
-               (for/vector #:length (vector-length dss) ([ds  (in-vector dss)])
-                 (for/vector #:length (vector-length ds) ([d  (in-vector ds)])
-                   (exact->inexact d))))
-             d-min d-max))
-
-(defproc (flonum-ok-for-2d? [x-min rational?] [x-max rational?]
-                            [y-min rational?] [y-max rational?]) boolean?
-  (and (flonum-ok-for-range? x-min x-max 10000)
-       (flonum-ok-for-range? y-min y-max 10000)))
-
-(defproc (flonum-ok-for-3d? [x-min rational?] [x-max rational?]
-                            [y-min rational?] [y-max rational?]
-                            [z-min rational?] [z-max rational?]) boolean?
-  (and (flonum-ok-for-range? x-min x-max 10000)
-       (flonum-ok-for-range? y-min y-max 10000)
-       (flonum-ok-for-range? z-min z-max 10000)))
-
-(defproc (flonum-ok-for-4d? [x-min rational?] [x-max rational?]
-                            [y-min rational?] [y-max rational?]
-                            [z-min rational?] [z-max rational?]
-                            [d-min rational?] [d-max rational?]) boolean?
-  (and (flonum-ok-for-range? x-min x-max 10000)
-       (flonum-ok-for-range? y-min y-max 10000)
-       (flonum-ok-for-range? z-min z-max 10000)
-       (flonum-ok-for-range? d-min d-max 10000)))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt
index 8ecafffeef..4cf0305255 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/math.rkt
@@ -5,7 +5,7 @@
 (require "../common/math.rkt")
 (provide equal?*
          ;; Flonums
-         flblend flsum fldistance (activate-contract-out flonum-ok-for-range?)
+         flblend flsum fldistance
          ;; Reals
          maybe-inexact->exact
          min* max* blend atan2 sum real-modulo distance
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt
index 98b86146d0..34c466cdb5 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/contracted/sample.rkt
@@ -24,12 +24,6 @@
                                 make-function->sampler
                                 make-2d-function->sampler
                                 make-3d-function->sampler
-                                sample-exact->inexact
-                                2d-sample-exact->inexact
-                                3d-sample-exact->inexact
-                                flonum-ok-for-2d?
-                                flonum-ok-for-3d?
-                                flonum-ok-for-4d?
                                 )
          (contract-out (struct mapped-function ([f (any/c . -> . any/c)]
                                                 [fmap ((listof any/c) . -> . (listof any/c))])))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt
index 8117af9e79..9ab9a7c737 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/clip.rkt
@@ -1,57 +1,89 @@
 #lang racket/base
 
-;; Small library for clipping points, rectangles, lines and polygons against axial planes.
+(require (for-syntax racket/base)
+         racket/match
+         racket/list)
 
-(require racket/match racket/list
-         "../common/utils.rkt")
-
-(provide point-in-bounds? clip-line clip-lines clip-polygon)
+(provide point-in-bounds?
+         point-inside-plane?
+         clip-line/bounds
+         clip-line/plane
+         clip-lines/bounds
+         clip-lines/plane
+         clip-polygon/bounds
+         clip-polygon/plane)
 
 ;; ===================================================================================================
-;; Point clipping
+;; Basic plane intersection and distance tests
+
+;; Need to use this instead of `*' to make axial plane tests as fast as hand-written because the
+;; optimizer doesn't always figure out (* 0 e) = 0, etc.
+(define-syntax (times stx)
+  (syntax-case stx ()
+    [(_ 0 e)  (syntax/loc stx 0)]
+    [(_ 1 e)  (syntax/loc stx e)]
+    [(_ 0.0 e)  (syntax/loc stx 0.0)]
+    [(_ 1.0 e)  (syntax/loc stx (real->double-flonum e))]
+    [(_ . es)  (syntax/loc stx (* . es))]))
+
+(define-syntax-rule (plane-line-intersect a b c x1-stx y1-stx x2-stx y2-stx)
+  (let ([x1 x1-stx] [y1 y1-stx] [x2 x2-stx] [y2 y2-stx])
+    (let ([dot1  (+ (times a x1) (times b y1))]
+          [dot2  (+ (times a x2) (times b y2))])
+      (let ([denom  (- dot1 dot2)])
+        (if (zero? denom)
+            (values x2 y2)
+            (let ([t  (/ (+ dot1 c) denom)])
+              (values (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1))))))))))
+
+(define-syntax-rule (plane-point-dist a b c x y)
+  (+ (times a x) (times b y) c))
+
+;; ===================================================================================================
+;; Points
 
 (define (point-in-bounds? v x-min x-max y-min y-max)
   (match-define (vector x y) v)
   (and (<= x-min x x-max) (<= y-min y y-max)))
 
+(define (point-inside-plane? v plane)
+  (match-define (vector x y) v)
+  (match-define (vector a b c) plane)
+  (>= (plane-point-dist a b c x y) 0))
+
 ;; ===================================================================================================
 ;; Line clipping
 
-(define (clip-line-x start-in-bounds? x x1 y1 x2 y2)
-  (let-map
-   (x x1 y1 x2 y2) inexact->exact
-   (define t (/ (- x x1) (- x2 x1)))
-   (cond [start-in-bounds?  (values x1 y1 x (+ y1 (* t (- y2 y1))))]
-         [else              (values x (+ y1 (* t (- y2 y1))) x2 y2)])))
+(define (clip-line/plane v1 v2 plane)
+  (match-define (vector x1 y1) v1)
+  (match-define (vector x2 y2) v2)
+  (match-define (vector a b c) plane)
+  (define v1? (>= (plane-point-dist a b c x1 y1) 0))
+  (define v2? (>= (plane-point-dist a b c x2 y2) 0))
+  (cond [(and v1? v2?)       (values v1 v2)]
+        [(not (or v1? v2?))  (values #f #f)]
+        [else
+         (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2))
+         (if v1? (values v1 (vector x y)) (values (vector x y) v2))]))
 
-(define (clip-line-y start-in-bounds? y x1 y1 x2 y2)
-  (let-map
-   (y x1 y1 x2 y2) inexact->exact
-   (define t (/ (- y y1) (- y2 y1)))
-   (cond [start-in-bounds?  (values x1 y1 (+ x1 (* t (- x2 x1))) y)]
-         [else              (values (+ x1 (* t (- x2 x1))) y x2 y2)])))
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
 
-(define (clip-line-x-min x-min x1 y1 x2 y2)
-  (cond [(and (x1 . >= . x-min) (x2 . < . x-min))  (clip-line-x #t x-min x1 y1 x2 y2)]
-        [(and (x2 . >= . x-min) (x1 . < . x-min))  (clip-line-x #f x-min x1 y1 x2 y2)]
-        [else                                      (values x1 y1 x2 y2)]))
+(define-syntax-rule (make-clip-line/axis a b gte?)
+  (λ (val x1 y1 x2 y2)
+    (define v1? (gte? (plane-point-dist a b (- val) x1 y1) 0))
+    (define v2? (gte? (plane-point-dist a b (- val) x2 y2) 0))
+    (cond [(or (and v1? v2?) (not (or v1? v2?))) (values x1 y1 x2 y2)]
+          [else
+           (define-values (x y) (plane-line-intersect a b (- val) x1 y1 x2 y2))
+           (if v1? (values x1 y1 x y) (values x y x2 y2))])))
 
-(define (clip-line-x-max x-max x1 y1 x2 y2)
-  (cond [(and (x1 . <= . x-max) (x2 . > . x-max))  (clip-line-x #t x-max x1 y1 x2 y2)]
-        [(and (x2 . <= . x-max) (x1 . > . x-max))  (clip-line-x #f x-max x1 y1 x2 y2)]
-        [else                                      (values x1 y1 x2 y2)]))
+(define clip-line-x-min (make-clip-line/axis 1 0 >=))
+(define clip-line-x-max (make-clip-line/axis 1 0 <=))
+(define clip-line-y-min (make-clip-line/axis 0 1 >=))
+(define clip-line-y-max (make-clip-line/axis 0 1 <=))
 
-(define (clip-line-y-min y-min x1 y1 x2 y2)
-  (cond [(and (y1 . >= . y-min) (y2 . < . y-min))  (clip-line-y #t y-min x1 y1 x2 y2)]
-        [(and (y2 . >= . y-min) (y1 . < . y-min))  (clip-line-y #f y-min x1 y1 x2 y2)]
-        [else                                      (values x1 y1 x2 y2)]))
-
-(define (clip-line-y-max y-max x1 y1 x2 y2)
-  (cond [(and (y1 . <= . y-max) (y2 . > . y-max))  (clip-line-y #t y-max x1 y1 x2 y2)]
-        [(and (y2 . <= . y-max) (y1 . > . y-max))  (clip-line-y #f y-max x1 y1 x2 y2)]
-        [else                                      (values x1 y1 x2 y2)]))
-
-(define (clip-line v1 v2 x-min x-max y-min y-max)
+(define (clip-line/bounds v1 v2 x-min x-max y-min y-max)
   (let/ec return
     (match-define (vector x1 y1) v1)
     (match-define (vector x2 y2) v2)
@@ -70,90 +102,127 @@
       (values (vector x1 y1) (vector x2 y2)))))
 
 ;; ===================================================================================================
-;; Polygon clipping
+;; Connected lines clipping
 
-(define-syntax-rule (make-clip-polygon in-bounds? clip-line)
-  (λ (val xs ys)
-    (cond [(empty? xs)  (values empty empty)]
-          [else
-           (for/fold ([res-xs empty] [res-ys empty]) ([x1  (in-list (cons (last xs) xs))]
-                                                      [x2  (in-list xs)]
-                                                      [y1  (in-list (cons (last ys) ys))]
-                                                      [y2  (in-list ys)])
-             (define v1-in-bounds? (in-bounds? x1 y1 val))
-             (define v2-in-bounds? (in-bounds? x2 y2 val))
-             (cond [(and v1-in-bounds? v2-in-bounds?)              (values (cons x2 res-xs)
-                                                                           (cons y2 res-ys))]
-                   [(and (not v1-in-bounds?) (not v2-in-bounds?))  (values res-xs res-ys)]
-                   [else  (let-values ([(x1 y1 x2 y2)  (clip-line v1-in-bounds? val x1 y1 x2 y2)])
-                            (cond [v2-in-bounds?  (values (list* x2 x1 res-xs)
-                                                          (list* y2 y1 res-ys))]
-                                  [else           (values (cons x2 res-xs) (cons y2 res-ys))]))]))])))
+(define-syntax-rule (make-clip-lines/plane a b c gte?)
+  (λ (vs)
+    (cond
+      [(empty? vs)  empty]
+      [else
+       (define v1 (first vs))
+       (match-define (vector x1 y1) v1)
+       (define v1? (gte? (plane-point-dist a b c x1 y1) 0))
+       
+       (define init-vss (if v1? (list (list v1)) (list empty)))
+       
+       (define-values (vss _x1 _y1 _v1?)
+         (for/fold ([vss init-vss] [x1 x1] [y1 y1] [v1? v1?]) ([v2  (in-list (rest vs))])
+           (define x2 (vector-ref v2 0))
+           (define y2 (vector-ref v2 1))
+           (define v2? (gte? (plane-point-dist a b c x2 y2) 0))
+           (cond [(and v1? v2?)       (values (cons (cons v2 (first vss)) (rest vss)) x2 y2 v2?)]
+                 [(not (or v1? v2?))  (values vss x2 y2 v2?)]
+                 [else
+                  (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2))
+                  (if v1?
+                      (values (cons (cons (vector x y) (first vss)) (rest vss)) x2 y2 v2?)
+                      (values (cons (list v2 (vector x y)) vss) x2 y2 v2?))])))
+       
+       (filter (compose not empty?) vss)])))
 
-(define-syntax-rule (x-min-in-bounds? x y x-min) (x . >= . x-min)) 
-(define-syntax-rule (x-max-in-bounds? x y x-max) (x . <= . x-max))
-(define-syntax-rule (y-min-in-bounds? x y y-min) (y . >= . y-min))
-(define-syntax-rule (y-max-in-bounds? x y y-max) (y . <= . y-max))
+(define (clip-lines/plane vs plane)
+  (match-define (vector a b c) plane)
+  ((make-clip-lines/plane a b c >=) vs))
 
-(define clip-polygon-x-min (make-clip-polygon x-min-in-bounds? clip-line-x))
-(define clip-polygon-x-max (make-clip-polygon x-max-in-bounds? clip-line-x))
-(define clip-polygon-y-min (make-clip-polygon y-min-in-bounds? clip-line-y))
-(define clip-polygon-y-max (make-clip-polygon y-max-in-bounds? clip-line-y))
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
 
-(define (clip-polygon vs x-min x-max y-min y-max)
+;; Early accept: all endpoints in bounds (or empty vs)
+(define (early-accept? vs x-min x-max y-min y-max)
+  (andmap (λ (v) (and (<= x-min (vector-ref v 0) x-max)
+                      (<= y-min (vector-ref v 1) y-max)))
+          vs))
+
+;; Early reject: all endpoints on the outside of the same plane
+(define (early-reject? vs x-min x-max y-min y-max)
+  (or (andmap (λ (v) ((vector-ref v 0) . < . x-min)) vs)
+      (andmap (λ (v) ((vector-ref v 0) . > . x-max)) vs)
+      (andmap (λ (v) ((vector-ref v 1) . < . y-min)) vs)
+      (andmap (λ (v) ((vector-ref v 1) . > . y-max)) vs)))
+
+(define-syntax-rule (make-clip-lines/axis a b gte?)
+  (λ (val vs)
+    ((make-clip-lines/plane a b (- val) gte?) vs)))
+
+(define clip-lines-x-min (make-clip-lines/axis 1 0 >=))
+(define clip-lines-x-max (make-clip-lines/axis 1 0 <=))
+(define clip-lines-y-min (make-clip-lines/axis 0 1 >=))
+(define clip-lines-y-max (make-clip-lines/axis 0 1 <=))
+
+(define (clip-lines/bounds vs x-min x-max y-min y-max)
   (let/ec return
-    ;; early reject: no polygon
-    (when (empty? vs) (return empty))
-    (match-define (list (vector xs ys) ...) vs)
-    ;; early accept: all endpoints in bounds
-    (when (and (andmap (λ (x) (<= x-min x x-max)) xs)
-               (andmap (λ (y) (<= y-min y y-max)) ys))
-      (return vs))
-    ;; early reject: all endpoints on the outside of the same plane
-    (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs)
-              (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys))
-      (return empty))
-    (let*-values ([(xs ys)  (clip-polygon-x-min x-min xs ys)]
-                  [(xs ys)  (clip-polygon-x-max x-max xs ys)]
-                  [(xs ys)  (clip-polygon-y-min y-min xs ys)]
-                  [(xs ys)  (clip-polygon-y-max y-max xs ys)])
-      (map vector xs ys))))
+    (when (early-accept? vs x-min x-max y-min y-max) (return (list vs)))
+    (when (early-reject? vs x-min x-max y-min y-max) (return empty))
+    (let* ([vss  (clip-lines-x-min x-min vs)]
+           [_    (when (empty? vss) (return empty))]
+           [vss  (append* (map (λ (vs) (clip-lines-x-max x-max vs)) vss))]
+           [_    (when (empty? vss) (return empty))]
+           [vss  (append* (map (λ (vs) (clip-lines-y-min y-min vs)) vss))]
+           [_    (when (empty? vss) (return empty))]
+           [vss  (append* (map (λ (vs) (clip-lines-y-max y-max vs)) vss))])
+      vss)))
 
 ;; ===================================================================================================
-;; Lines clipping
+;; Polygon clipping
 
-(define (join-lines lines)
-  (let loop ([lines lines] [current-line empty])
-    (cond [(empty? lines)         (list (reverse current-line))]
-          [(empty? current-line)  (loop (rest lines) (reverse (first lines)))]
-          [else
-           (match-define v (first current-line))
-           (match-define (list v1 v2) (first lines))
-           (cond [(equal? v v1)  (loop (rest lines) (cons v2 current-line))]
-                 [else           (cons (reverse current-line) (loop lines empty))])])))
+(define-syntax-rule (make-clip-polygon/plane a b c gte?)
+  (λ (vs ls)
+    (define v1 (last vs))
+    (define x1 (vector-ref v1 0))
+    (define y1 (vector-ref v1 1))
+    (define v1? (gte? (plane-point-dist a b c x1 y1) 0))
+    
+    (define-values (new-vs new-ls _x1 _y1 _v1?)
+      (for/fold ([vs empty] [ls empty] [x1 x1] [y1 y1] [v1? v1?]) ([v2  (in-list vs)]
+                                                                   [l   (in-list ls)])
+        (define x2 (vector-ref v2 0))
+        (define y2 (vector-ref v2 1))
+        (define v2? (gte? (plane-point-dist a b c x2 y2) 0))
+        (cond [(and v1? v2?)       (values (cons v2 vs) (cons l ls) x2 y2 v2?)]
+              [(not (or v1? v2?))  (values vs ls x2 y2 v2?)]
+              [else
+               (define-values (x y) (plane-line-intersect a b c x1 y1 x2 y2))
+               (if v1?
+                   (values (cons (vector x y) vs) (cons l ls) x2 y2 v2?)
+                   (values (list* v2 (vector x y) vs) (list* l #t ls) x2 y2 v2?))])))
+    
+    (values (reverse new-vs) (reverse new-ls))))
 
-#;
-(join-lines
- '((#(1 2) #(3 4))
-   (#(3 4) #(5 6))
-   (#(5 7) #(7 8))
-   (#(7 8) #(9 10))))
+(define (clip-polygon/plane vs ls plane)
+  (match-define (vector a b c) plane)
+  ((make-clip-polygon/plane a b c >=) vs ls))
 
-(define (clip-lines vs x-min x-max y-min y-max)
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
+
+(define-syntax-rule (make-clip-polygon/axis a b gte?)
+  (λ (vs ls val)
+    ((make-clip-polygon/plane a b (- val) gte?) vs ls)))
+
+(define clip-polygon-x-min (make-clip-polygon/axis 1 0 >=))
+(define clip-polygon-x-max (make-clip-polygon/axis 1 0 <=))
+(define clip-polygon-y-min (make-clip-polygon/axis 0 1 >=))
+(define clip-polygon-y-max (make-clip-polygon/axis 0 1 <=))
+
+(define (clip-polygon/bounds vs ls x-min x-max y-min y-max)
   (let/ec return
-    ;; early reject: no lines
-    (when (empty? vs) (return empty))
-    (match-define (list (vector xs ys) ...) vs)
-    ;; early accept: all endpoints in bounds
-    (when (and (andmap (λ (x) (<= x-min x x-max)) xs) (andmap (λ (y) (<= y-min y y-max)) ys))
-      (return (list vs)))
-    ;; early reject: all endpoints on the outside of the same plane
-    (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs)
-              (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys))
-      (return empty))
-    (join-lines
-     (reverse
-      (for/fold ([res empty]) ([v1  (in-list vs)] [v2  (in-list (rest vs))])
-        (let-values ([(v1 v2)  (clip-line v1 v2 x-min x-max y-min y-max)])
-          (cond [(and v1 v2)  (cons (list v1 v2) res)]
-                [else  res])))))))
+    (when (early-accept? vs x-min x-max y-min y-max) (return vs ls))
+    (when (early-reject? vs x-min x-max y-min y-max) (return empty empty))
+    (let*-values ([(vs ls)  (clip-polygon-x-min vs ls x-min)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-x-max vs ls x-max)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-y-min vs ls y-min)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-y-max vs ls y-max)])
+      (values vs ls))))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt
index 2a0169906c..d35e04673e 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/contour.rkt
@@ -2,7 +2,7 @@
 
 ;; Renderers for contour lines and contour intervals
 
-(require racket/contract racket/class racket/match racket/list racket/flonum racket/vector racket/math
+(require racket/contract racket/class racket/match racket/list racket/vector
          unstable/latent-contract/defthing
          plot/utils
          "../common/utils.rkt")
@@ -23,13 +23,10 @@
   (when (<= z-min z z-max)
     (send area put-alpha alpha)
     (send area put-pen color width style)
-    (let* ([flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-           [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-           [z           (if flonum-ok? (exact->inexact z) z)])
-      (for-2d-sample
-       (xa xb ya yb z1 z2 z3 z4) sample
-       (for/list ([line  (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))])
-         (send/apply area put-line (map (λ (v) (vector-take v 2)) line))))))
+    (for-2d-sample
+     (xa xb ya yb z1 z2 z3 z4) sample
+     (for/list ([line  (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))])
+       (send/apply area put-line (map (λ (v) (vector-take v 2)) line)))))
   
   (cond [label  (line-legend-entry label color width style)]
         [else   empty]))
@@ -72,10 +69,7 @@
       (let* ([colors  (maybe-apply colors zs)]
              [widths  (maybe-apply widths zs)]
              [styles  (maybe-apply styles zs)]
-             [alphas  (maybe-apply alphas zs)]
-             [flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-             [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-             [zs          (if flonum-ok? (map exact->inexact zs) zs)])
+             [alphas  (maybe-apply alphas zs)])
         (for ([z      (in-list zs)]
               [color  (in-cycle colors)]
               [width  (in-cycle widths)]
@@ -138,10 +132,7 @@
     (send area put-pen 0 1 'transparent)
     (let* ([colors  (map ->brush-color (maybe-apply colors z-ivls))]
            [styles  (map ->brush-style (maybe-apply styles z-ivls))]
-           [alphas  (maybe-apply alphas z-ivls)]
-           [flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-           [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-           [zs          (if flonum-ok? (map exact->inexact zs) zs)])
+           [alphas  (maybe-apply alphas z-ivls)])
       (for ([za     (in-list zs)]
             [zb     (in-list (rest zs))]
             [color  (in-cycle colors)]
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt
index 40c6726f15..56435b9da7 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/decoration.rkt
@@ -7,9 +7,6 @@
          unstable/contract
          plot/utils
          "../common/utils.rkt"
-         "line.rkt"
-         "interval.rkt"
-         "point.rkt"
          "clip.rkt")
 
 (provide (all-defined-out))
@@ -98,8 +95,8 @@
   (define-values (r-mins r-maxs)
     (for/lists (r-mins r-maxs) ([θ  (in-list θs)])
       (define-values (v1 v2)
-        (clip-line (vector 0 0) (vector (* max-r (cos θ)) (* max-r (sin θ)))
-                   x-min x-max y-min y-max))
+        (clip-line/bounds (vector 0 0) (vector (* max-r (cos θ)) (* max-r (sin θ)))
+                          x-min x-max y-min y-max))
       (values (if v1 (vmag v1) #f)
               (if v2 (vmag v2) #f))))
   (for/lists (θs r-mins r-maxs) ([θ  (in-list θs)] [r-min  (in-list r-mins)] [r-max  (in-list r-maxs)]
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt
index 59e4a41cf4..590b41cd80 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/interval.rkt
@@ -2,7 +2,7 @@
 
 ;; Renderers for intervals between functions.
 
-(require racket/contract racket/class racket/match racket/math racket/list racket/sequence
+(require racket/contract racket/class racket/match racket/math racket/list
          unstable/latent-contract/defthing
          unstable/contract
          plot/utils
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt
index e86eee5b41..77ff41e35e 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/plot-area.rkt
@@ -1,18 +1,16 @@
 #lang racket/base
 
-(require racket/draw racket/class racket/contract racket/match racket/math racket/list racket/string
-         racket/flonum
+(require racket/class racket/match racket/math racket/list racket/flonum
+         (only-in math/flonum fl)
          "../common/plot-device.rkt"
          "../common/ticks.rkt"
-         "../common/contract.rkt"
          "../common/math.rkt"
          "../common/draw.rkt"
          "../common/axis-transform.rkt"
-         "../common/sample.rkt"
-         "../common/legend.rkt"
          "../common/parameters.rkt"
          "../common/utils.rkt"
-         "clip.rkt")
+         "clip.rkt"
+         "vector.rkt")
 
 (provide (all-defined-out))
 
@@ -87,45 +85,20 @@
       (and (equal? (plot-x-transform) id-transform)
            (equal? (plot-y-transform) id-transform)))
     
-    (define flonum-ok? (flonum-ok-for-2d? x-min x-max y-min y-max))
-    (when flonum-ok?
-      (set! x-min (exact->inexact x-min))
-      (set! x-max (exact->inexact x-max))
-      (set! y-min (exact->inexact y-min))
-      (set! y-max (exact->inexact y-max))
-      (set! x-size (exact->inexact x-size))
-      (set! y-size (exact->inexact y-size))
-      (set! x-mid (exact->inexact x-mid))
-      (set! y-mid (exact->inexact y-mid)))
-    
     (define plot->view
-      (if flonum-ok?
-          (let-map
-           (x-min y-min x-size y-size) exact->inexact
-           (if identity-transforms?
-               (match-lambda
-                 [(vector x y)
-                  (vector (fl/ (fl- (exact->inexact x) x-min) x-size)
-                          (fl/ (fl- (exact->inexact y) y-min) y-size))])
-               (match-lambda
-                 [(vector (? rational? x) (? rational? y))
-                  (vector (fl/ (fl- (exact->inexact (fx x)) x-min) x-size)
-                          (fl/ (fl- (exact->inexact (fy y)) y-min) y-size))]
-                 [(vector x y)
-                  (vector +nan.0 +nan.0)])))
-          (if identity-transforms?
-              (match-lambda
-                [(vector (? rational? x) (? rational? y))
-                 (vector (exact->inexact (/ (- (inexact->exact x) x-min) x-size))
-                         (exact->inexact (/ (- (inexact->exact y) y-min) y-size)))]
-                [(vector x y)
-                 (vector +nan.0 +nan.0)])
-              (match-lambda
-                [(vector (? rational? x) (? rational? y))
-                 (vector (exact->inexact (/ (- (inexact->exact (fx x)) x-min) x-size))
-                         (exact->inexact (/ (- (inexact->exact (fy y)) y-min) y-size)))]
-                [(vector x y)
-                 (vector +nan.0 +nan.0)]))))
+      (if identity-transforms?
+          (match-lambda
+            [(vector (? rational? x) (? rational? y))
+             (vector (fl (/ (- (inexact->exact x) x-min) x-size))
+                     (fl (/ (- (inexact->exact y) y-min) y-size)))]
+            [(vector x y)
+             (vector +nan.0 +nan.0)])
+          (match-lambda
+            [(vector (? rational? x) (? rational? y))
+             (vector (fl (/ (- (inexact->exact (fx x)) x-min) x-size))
+                     (fl (/ (- (inexact->exact (fy y)) y-min) y-size)))]
+            [(vector x y)
+             (vector +nan.0 +nan.0)])))
     
     (define view->dc #f)
     (define (plot->dc v) (view->dc (plot->view v)))
@@ -144,7 +117,7 @@
       (define area-per-view-x (/ (- area-x-max area-x-min) view-x-size))
       (define area-per-view-y (/ (- area-y-max area-y-min) view-y-size))
       (let-map
-       (area-x-min area-y-max area-per-view-x area-per-view-y) exact->inexact
+       (area-x-min area-y-max area-per-view-x area-per-view-y) fl
        (λ (v)
          (match-define (vector x y) v)
          (vector (fl+ area-x-min (fl* x area-per-view-x))
@@ -503,61 +476,69 @@
     ;; Shapes
     
     (define/public (put-lines vs)
-      (for ([vs  (vrational-sublists vs)])
-        (for ([vs  (if clipping?
-                       (in-list (clip-lines vs clip-x-min clip-x-max
-                                            clip-y-min clip-y-max))
-                       (in-value vs))])
-          (when (not (empty? vs))
-            (let* ([vs  (if identity-transforms? vs (subdivide-lines plot->dc vs))]
-                   [vs  (map (λ (v) (plot->dc v)) vs)])
-              (send pd draw-lines vs))))))
+      (for ([vs  (in-list (exact-vector2d-sublists vs))])
+        (let ([vss  (if clipping?
+                        (clip-lines/bounds vs
+                                           clip-x-min clip-x-max
+                                           clip-y-min clip-y-max)
+                        (list vs))])
+          (for ([vs  (in-list vss)])
+            (unless (empty? vs)
+              (let* ([vs  (if identity-transforms? vs (subdivide-lines plot->dc vs))]
+                     [vs  (map plot->dc vs)])
+                (send pd draw-lines vs)))))))
     
     (define/public (put-line v1 v2)
-      (when (and (vrational? v1) (vrational? v2))
-        (let-values ([(v1 v2)  (if clipping?
-                                   (clip-line v1 v2 clip-x-min clip-x-max
-                                              clip-y-min clip-y-max)
-                                   (values v1 v2))])
-          (when (and v1 v2)
-            (if identity-transforms?
-                (send pd draw-line (plot->dc v1) (plot->dc v2))
-                (send pd draw-lines (map (λ (v) (plot->dc v))
-                                         (subdivide-line plot->dc v1 v2))))))))
+      (let ([v1  (exact-vector2d v1)]
+            [v2  (exact-vector2d v2)])
+        (when (and v1 v2)
+          (let-values ([(v1 v2)  (if clipping?
+                                     (clip-line/bounds v1 v2
+                                                       clip-x-min clip-x-max
+                                                       clip-y-min clip-y-max)
+                                     (values v1 v2))])
+            (when (and v1 v2)
+              (if identity-transforms?
+                  (send pd draw-line (plot->dc v1) (plot->dc v2))
+                  (send pd draw-lines (map plot->dc (subdivide-line plot->dc v1 v2)))))))))
     
     (define/public (put-polygon vs)
-      (when (andmap vrational? vs)
-        (let* ([vs  (if clipping?
-                        (clip-polygon vs clip-x-min clip-x-max
-                                      clip-y-min clip-y-max)
-                        vs)])
-          (when (not (empty? vs))
-            (if identity-transforms?
-                (send pd draw-polygon (map (λ (v) (plot->dc v)) vs))
-                (send pd draw-polygon (map (λ (v) (plot->dc v))
-                                           (subdivide-polygon plot->dc vs))))))))
+      (define ls (make-list (length vs) #t))
+      (let-values ([(vs ls)  (exact-polygon2d vs ls)])
+        (unless (empty? vs)
+          (let-values ([(vs ls)  (if clipping?
+                                     (clip-polygon/bounds vs ls
+                                                          clip-x-min clip-x-max
+                                                          clip-y-min clip-y-max)
+                                     (values vs ls))])
+            (unless (empty? vs)
+              (if identity-transforms?
+                  (send pd draw-polygon (map plot->dc vs))
+                  (send pd draw-polygon (map plot->dc (subdivide-polygon plot->dc vs)))))))))
     
     (define/public (put-rect r)
       (when (rect-rational? r)
         (match-define (vector (ivl x1 x2) (ivl y1 y2)) r)
         (put-polygon (list (vector x1 y1) (vector x2 y1) (vector x2 y2) (vector x1 y2)))))
     
-    (define/public (put-text str v [anchor 'top-left] [angle 0] [dist 0]
-                             #:outline? [outline? #f])
-      (when (and (vrational? v) (in-bounds? v))
-        (send pd draw-text str (plot->dc v) anchor angle dist #:outline? outline?)))
+    (define/public (put-text str v [anchor 'top-left] [angle 0] [dist 0] #:outline? [outline? #f])
+      (let ([v  (exact-vector2d v)])
+        (when (and v (in-bounds? v))
+          (send pd draw-text str (plot->dc v) anchor angle dist #:outline? outline?))))
     
     (define/public (put-glyphs vs symbol size)
-      (send pd draw-glyphs (map (λ (v) (plot->dc v))
-                                (filter (λ (v) (and (vrational? v) (in-bounds? v)))
-                                        vs))
-            symbol size))
-    
+      (let ([vs  (filter (λ (v) (and v (in-bounds? v))) (map exact-vector2d vs))])
+        (unless (empty? vs)
+          (send pd draw-glyphs (map plot->dc vs) symbol size))))
+      
     (define/public (put-arrow v1 v2)
-      (when (and (vrational? v1) (vrational? v2) (in-bounds? v1))
-        (send pd draw-arrow (plot->dc v1) (plot->dc v2))))
+      (let ([v1  (exact-vector2d v1)]
+            [v2  (exact-vector2d v2)])
+        (when (and v1 v2 (in-bounds? v1))
+          (send pd draw-arrow (plot->dc v1) (plot->dc v2)))))
     
     (define/public (put-tick v r angle)
-      (when (and (vrational? v) (in-bounds? v))
-        (send pd draw-tick (plot->dc v) r angle)))
+      (let ([v  (exact-vector2d v)])
+        (when (and v (in-bounds? v))
+          (send pd draw-tick (plot->dc v) r angle))))
     ))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt
index daf87e5261..2db80ac1ce 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/point.rkt
@@ -2,7 +2,7 @@
 
 ;; Renderers for points and other point-like things.
 
-(require racket/contract racket/class racket/match racket/math racket/list racket/sequence
+(require racket/contract racket/class racket/match racket/math racket/list
          unstable/latent-contract/defthing
          unstable/contract
          plot/utils
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt
new file mode 100644
index 0000000000..aed08755e3
--- /dev/null
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot2d/vector.rkt
@@ -0,0 +1,46 @@
+#lang racket/base
+
+(require racket/list
+         (only-in racket/unsafe/ops unsafe-vector-ref))
+
+(provide exact-vector2d
+         exact-vector2d-sublists
+         exact-polygon2d)
+
+(define (exact-vector2d v)
+  (define n (vector-length v))
+  (cond [(= n 2)
+         (define v1 (unsafe-vector-ref v 0))
+         (define v2 (unsafe-vector-ref v 1))
+         (cond [(and (exact? v1) (exact? v2))  v]
+               [(and (rational? v1) (rational? v2))
+                (vector (inexact->exact v1) (inexact->exact v2))]
+               [else  #f])]
+        [else  #f]))
+
+(define (sublists vs)
+  (define vss
+    (for/fold ([vss  (list null)]) ([v  (in-list vs)])
+      (cond [v  (cons (cons v (car vss)) (cdr vss))]
+            [(null? (car vss))  vss]
+            [else  (cons null vss)])))
+  (cond [(null? (car vss))  (cdr vss)]
+        [else  vss]))
+
+(define (exact-vector2d-sublists vs)
+  (sublists (map exact-vector2d vs)))
+
+(define (exact-polygon2d vs ls)
+  (cond
+    [(null? vs)  (values null null)]
+    [else
+     (define-values (new-vs new-ls _)
+       (for/fold ([vs null] [ls null] [v1 (last vs)]) ([v2  (in-list vs)]
+                                                       [l   (in-list ls)])
+         (cond [(equal? v1 v2)  (values vs ls v2)]
+               [else
+                (define exact-v2 (exact-vector2d v2))
+                (if exact-v2
+                    (values (cons exact-v2 vs) (cons l ls) v2)
+                    (values vs ls v2))])))
+     (values (reverse new-vs) (reverse new-ls))]))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt
index 9c29b76f67..ff136c1e0f 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/bsp.rkt
@@ -3,9 +3,7 @@
 (require racket/list
          racket/match
          racket/promise
-         racket/math
          math/flonum
-         math/statistics
          "split.rkt")
 
 (provide
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt
index 9406621c41..fe75b7e0b7 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/clip.rkt
@@ -1,14 +1,45 @@
 #lang racket/base
 
-;; Small library for clipping points, lines and polygons against axial planes.
+(require (for-syntax racket/base)
+         racket/match
+         racket/list)
 
-(require racket/match racket/list racket/unsafe/ops
-         "../common/utils.rkt")
+(provide point-in-bounds?
+         point-inside-plane?
+         clip-line/bounds
+         clip-line/plane
+         clip-lines/bounds
+         clip-lines/plane
+         clip-polygon/bounds
+         clip-polygon/plane)
 
-(provide point-in-bounds? clip-line clip-polygon clip-lines
-         clip-polygon-x-min clip-polygon-x-max
-         clip-polygon-y-min clip-polygon-y-max
-         clip-polygon-z-min clip-polygon-z-max)
+;; ===================================================================================================
+;; Basic plane intersection and distance tests
+
+;; Need to use this instead of `*' to make axial plane tests as fast as hand-written because the
+;; optimizer doesn't always figure out (* 0 e) = 0, etc.
+(define-syntax (times stx)
+  (syntax-case stx ()
+    [(_ 0 e)  (syntax/loc stx 0)]
+    [(_ 1 e)  (syntax/loc stx e)]
+    [(_ 0.0 e)  (syntax/loc stx 0.0)]
+    [(_ 1.0 e)  (syntax/loc stx (real->double-flonum e))]
+    [(_ . es)  (syntax/loc stx (* . es))]))
+
+(define-syntax-rule (plane-line-intersect a b c d x1-stx y1-stx z1-stx x2-stx y2-stx z2-stx)
+  (let ([x1 x1-stx] [y1 y1-stx] [z1 z1-stx] [x2 x2-stx] [y2 y2-stx] [z2 z2-stx])
+    (let ([dot1  (+ (times a x1) (times b y1) (times c z1))]
+          [dot2  (+ (times a x2) (times b y2) (times c z2))])
+      (let ([denom  (- dot1 dot2)])
+        (if (zero? denom)
+            (values x2 y2 z2)
+            (let ([t  (/ (+ dot1 d) denom)])
+              (values (+ x1 (* t (- x2 x1)))
+                      (+ y1 (* t (- y2 y1)))
+                      (+ z1 (* t (- z2 z1))))))))))
+
+(define-syntax-rule (plane-point-dist a b c d x y z)
+  (+ (times a x) (times b y) (times c z) d))
 
 ;; ===================================================================================================
 ;; Points
@@ -17,61 +48,46 @@
   (match-define (vector x y z) v)
   (and (<= x-min x x-max) (<= y-min y y-max) (<= z-min z z-max)))
 
+(define (point-inside-plane? v plane)
+  (match-define (vector x y z) v)
+  (match-define (vector a b c d) plane)
+  (>= (plane-point-dist a b c d x y z) 0))
+
 ;; ===================================================================================================
 ;; Line clipping
 
-(define (clip-line-x start-in-bounds? x x1 y1 z1 x2 y2 z2)
-  (let-map
-   (x x1 y1 z1 x2 y2 z2) inexact->exact
-   (define t (/ (- x x1) (- x2 x1)))
-   (cond [start-in-bounds?  (values x1 y1 z1 x (+ y1 (* t (- y2 y1))) (+ z1 (* t (- z2 z1))))]
-         [else              (values x (+ y1 (* t (- y2 y1))) (+ z1 (* t (- z2 z1))) x2 y2 z2)])))
+(define (clip-line/plane v1 v2 plane)
+  (match-define (vector x1 y1 z1) v1)
+  (match-define (vector x2 y2 z2) v2)
+  (match-define (vector a b c d) plane)
+  (define v1? (>= (plane-point-dist a b c d x1 y1 z1) 0))
+  (define v2? (>= (plane-point-dist a b c d x2 y2 z2) 0))
+  (cond [(and v1? v2?)       (values v1 v2)]
+        [(not (or v1? v2?))  (values #f #f)]
+        [else
+         (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2))
+         (if v1? (values v1 (vector x y z)) (values (vector x y z) v2))]))
 
-(define (clip-line-y start-in-bounds? y x1 y1 z1 x2 y2 z2)
-  (let-map
-   (y x1 y1 z1 x2 y2 z2) inexact->exact
-   (define t (/ (- y y1) (- y2 y1)))
-   (cond [start-in-bounds?  (values x1 y1 z1 (+ x1 (* t (- x2 x1))) y (+ z1 (* t (- z2 z1))))]
-         [else              (values (+ x1 (* t (- x2 x1))) y (+ z1 (* t (- z2 z1))) x2 y2 z2)])))
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
 
-(define (clip-line-z start-in-bounds? z x1 y1 z1 x2 y2 z2)
-  (let-map
-   (z x1 y1 z1 x2 y2 z2) inexact->exact
-   (define t (/ (- z z1) (- z2 z1)))
-   (cond [start-in-bounds?  (values x1 y1 z1 (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1))) z)]
-         [else              (values (+ x1 (* t (- x2 x1))) (+ y1 (* t (- y2 y1))) z x2 y2 z2)])))
+(define-syntax-rule (make-clip-line/axis a b c gte?)
+  (λ (val x1 y1 z1 x2 y2 z2)
+    (define v1? (gte? (plane-point-dist a b c (- val) x1 y1 z1) 0))
+    (define v2? (gte? (plane-point-dist a b c (- val) x2 y2 z2) 0))
+    (cond [(or (and v1? v2?) (not (or v1? v2?))) (values x1 y1 z1 x2 y2 z2)]
+          [else
+           (define-values (x y z) (plane-line-intersect a b c (- val) x1 y1 z1 x2 y2 z2))
+           (if v1? (values x1 y1 z1 x y z) (values x y z x2 y2 z2))])))
 
-(define (clip-line-x-min x-min x1 y1 z1 x2 y2 z2)
-  (cond [(and (x1 . >= . x-min) (x2 . < . x-min))  (clip-line-x #t x-min x1 y1 z1 x2 y2 z2)]
-        [(and (x2 . >= . x-min) (x1 . < . x-min))  (clip-line-x #f x-min x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
+(define clip-line-x-min (make-clip-line/axis 1 0 0 >=))
+(define clip-line-x-max (make-clip-line/axis 1 0 0 <=))
+(define clip-line-y-min (make-clip-line/axis 0 1 0 >=))
+(define clip-line-y-max (make-clip-line/axis 0 1 0 <=))
+(define clip-line-z-min (make-clip-line/axis 0 0 1 >=))
+(define clip-line-z-max (make-clip-line/axis 0 0 1 <=))
 
-(define (clip-line-x-max x-max x1 y1 z1 x2 y2 z2)
-  (cond [(and (x1 . <= . x-max) (x2 . > . x-max))  (clip-line-x #t x-max x1 y1 z1 x2 y2 z2)]
-        [(and (x2 . <= . x-max) (x1 . > . x-max))  (clip-line-x #f x-max x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
-
-(define (clip-line-y-min y-min x1 y1 z1 x2 y2 z2)
-  (cond [(and (y1 . >= . y-min) (y2 . < . y-min))  (clip-line-y #t y-min x1 y1 z1 x2 y2 z2)]
-        [(and (y2 . >= . y-min) (y1 . < . y-min))  (clip-line-y #f y-min x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
-
-(define (clip-line-y-max y-max x1 y1 z1 x2 y2 z2)
-  (cond [(and (y1 . <= . y-max) (y2 . > . y-max))  (clip-line-y #t y-max x1 y1 z1 x2 y2 z2)]
-        [(and (y2 . <= . y-max) (y1 . > . y-max))  (clip-line-y #f y-max x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
-
-(define (clip-line-z-min z-min x1 y1 z1 x2 y2 z2)
-  (cond [(and (z1 . >= . z-min) (z2 . < . z-min))  (clip-line-z #t z-min x1 y1 z1 x2 y2 z2)]
-        [(and (z2 . >= . z-min) (z1 . < . z-min))  (clip-line-z #f z-min x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
-
-(define (clip-line-z-max z-max x1 y1 z1 x2 y2 z2)
-  (cond [(and (z1 . <= . z-max) (z2 . > . z-max))  (clip-line-z #t z-max x1 y1 z1 x2 y2 z2)]
-        [(and (z2 . <= . z-max) (z1 . > . z-max))  (clip-line-z #f z-max x1 y1 z1 x2 y2 z2)]
-        [else                                      (values x1 y1 z1 x2 y2 z2)]))
-
-(define (clip-line v1 v2 x-min x-max y-min y-max z-min z-max)
+(define (clip-line/bounds v1 v2 x-min x-max y-min y-max z-min z-max)
   (let/ec return
     (match-define (vector x1 y1 z1) v1)
     (match-define (vector x2 y2 z2) v2)
@@ -92,117 +108,74 @@
                   [(x1 y1 z1 x2 y2 z2)  (clip-line-z-max z-max x1 y1 z1 x2 y2 z2)])
       (values (vector x1 y1 z1) (vector x2 y2 z2)))))
 
-;; ===================================================================================================
-;; Polygon clipping
-
-(define-syntax-rule (make-clip-polygon idx test? clip-line)
-  (λ (val vs ls)
-    (define-values (new-vs new-ls)
-      (for/fold ([vs empty] [ls empty]) ([v1  (in-list (cons (last vs) vs))]
-                                         [v2  (in-list vs)]
-                                         [l   (in-list ls)])
-        (define v1-in-bounds? (test? (unsafe-vector-ref v1 idx) val))
-        (define v2-in-bounds? (test? (unsafe-vector-ref v2 idx) val))
-        (cond [(and v1-in-bounds? v2-in-bounds?)
-               (values (cons v2 vs) (cons l ls))]
-              [(and (not v1-in-bounds?) (not v2-in-bounds?))
-               (values vs ls)]
-              [else
-               (match-define (vector x1 y1 z1) v1)
-               (match-define (vector x2 y2 z2) v2)
-               (let-values ([(x1 y1 z1 x2 y2 z2) (clip-line v1-in-bounds? val x1 y1 z1 x2 y2 z2)])
-                 (cond [v2-in-bounds?
-                        (values (list* (vector x2 y2 z2) (vector x1 y1 z1) vs)
-                                (list* l #t ls))]
-                       [else
-                        (values (cons (vector x2 y2 z2) vs)
-                                (cons l ls))]))])))
-    (values (reverse new-vs) (reverse new-ls))))
-
-(define clip-polygon-x-min (make-clip-polygon 0 >= clip-line-x))
-(define clip-polygon-x-max (make-clip-polygon 0 <= clip-line-x))
-(define clip-polygon-y-min (make-clip-polygon 1 >= clip-line-y))
-(define clip-polygon-y-max (make-clip-polygon 1 <= clip-line-y))
-(define clip-polygon-z-min (make-clip-polygon 2 >= clip-line-z))
-(define clip-polygon-z-max (make-clip-polygon 2 <= clip-line-z))
-
-(define (clip-polygon vs ls x-min x-max y-min y-max z-min z-max)
-  (let/ec return
-    ;; early reject: no polygon
-    (when (empty? vs) (return empty empty))
-    (match-define (list (vector xs ys zs) ...) vs)
-    ;; early accept: all endpoints in bounds
-    (when (and (andmap (λ (x) (<= x-min x x-max)) xs)
-               (andmap (λ (y) (<= y-min y y-max)) ys)
-               (andmap (λ (z) (<= z-min z z-max)) zs))
-      (return vs ls))
-    ;; early reject: all endpoints on the outside of the same plane
-    (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs)
-              (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys)
-              (andmap (λ (z) (z . < . z-min)) zs) (andmap (λ (z) (z . > . z-max)) zs))
-      (return empty empty))
-    (let*-values ([(vs ls)  (clip-polygon-x-min x-min vs ls)]
-                  [(_)      (when (empty? vs) (return empty empty))]
-                  [(vs ls)  (clip-polygon-x-max x-max vs ls)]
-                  [(_)      (when (empty? vs) (return empty empty))]
-                  [(vs ls)  (clip-polygon-y-min y-min vs ls)]
-                  [(_)      (when (empty? vs) (return empty empty))]
-                  [(vs ls)  (clip-polygon-y-max y-max vs ls)]
-                  [(_)      (when (empty? vs) (return empty empty))]
-                  [(vs ls)  (clip-polygon-z-min z-min vs ls)]
-                  [(_)      (when (empty? vs) (return empty empty))]
-                  [(vs ls)  (clip-polygon-z-max z-max vs ls)])
-      (values vs ls))))
-
 ;; ===================================================================================================
 ;; Connected lines clipping
 
-(define-syntax-rule (make-clip-lines idx test? clip-line)
+(define-syntax-rule (make-clip-lines/plane a b c d gte?)
+  (λ (vs)
+    (cond
+      [(empty? vs)  empty]
+      [else
+       (define v1 (first vs))
+       (match-define (vector x1 y1 z1) v1)
+       (define v1? (gte? (plane-point-dist a b c d x1 y1 z1) 0))
+       
+       (define init-vss (if v1? (list (list v1)) (list empty)))
+       
+       (define-values (vss _x1 _y1 _z1 _v1?)
+         (for/fold ([vss init-vss] [x1 x1] [y1 y1] [z1 z1] [v1? v1?]) ([v2  (in-list (rest vs))])
+           (define x2 (vector-ref v2 0))
+           (define y2 (vector-ref v2 1))
+           (define z2 (vector-ref v2 2))
+           (define v2? (gte? (plane-point-dist a b c d x2 y2 z2) 0))
+           (cond [(and v1? v2?)       (values (cons (cons v2 (first vss)) (rest vss)) x2 y2 z2 v2?)]
+                 [(not (or v1? v2?))  (values vss x2 y2 z2 v2?)]
+                 [else
+                  (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2))
+                  (if v1?
+                      (values (cons (cons (vector x y z) (first vss)) (rest vss)) x2 y2 z2 v2?)
+                      (values (cons (list v2 (vector x y z)) vss) x2 y2 z2 v2?))])))
+       
+       (filter (compose not empty?) vss)])))
+
+(define (clip-lines/plane vs plane)
+  (match-define (vector a b c d) plane)
+  ((make-clip-lines/plane a b c d >=) vs))
+
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
+
+;; Early accept: all endpoints in bounds (or empty vs)
+(define (early-accept? vs x-min x-max y-min y-max z-min z-max)
+  (andmap (λ (v) (and (<= x-min (vector-ref v 0) x-max)
+                      (<= y-min (vector-ref v 1) y-max)
+                      (<= z-min (vector-ref v 2) z-max)))
+          vs))
+
+;; Early reject: all endpoints on the outside of the same plane
+(define (early-reject? vs x-min x-max y-min y-max z-min z-max)
+  (or (andmap (λ (v) ((vector-ref v 0) . < . x-min)) vs)
+      (andmap (λ (v) ((vector-ref v 0) . > . x-max)) vs)
+      (andmap (λ (v) ((vector-ref v 1) . < . y-min)) vs)
+      (andmap (λ (v) ((vector-ref v 1) . > . y-max)) vs)
+      (andmap (λ (v) ((vector-ref v 2) . < . z-min)) vs)
+      (andmap (λ (v) ((vector-ref v 2) . > . z-max)) vs)))
+
+(define-syntax-rule (make-clip-lines/axis a b c gte?)
   (λ (val vs)
-    (define v1 (first vs))
-    (define v1? (test? (unsafe-vector-ref v1 idx) val))
-    
-    (define-values (vss _v1 _v1?)
-      (for/fold ([vss (if v1? (list (list v1)) (list empty))]
-                 [v1 v1]
-                 [v1? v1?]
-                 ) ([v2  (in-list (rest vs))])
-        (define v2? (test? (unsafe-vector-ref v2 idx) val))
-        (cond [(and v1? v2?)       (values (cons (cons v2 (first vss)) (rest vss)) v2 v2?)]
-              [(not (or v1? v2?))  (values vss v2 v2?)]
-              [else
-               (match-define (vector x1 y1 z1) v1)
-               (match-define (vector x2 y2 z2) v2)
-               (let-values ([(x1 y1 z1 x2 y2 z2) (clip-line v1? val x1 y1 z1 x2 y2 z2)])
-                 (cond [v2?
-                        (values (cons (list (vector x2 y2 z2) (vector x1 y1 z1)) vss) v2 v2?)]
-                       [else
-                        (values (cons (cons (vector x2 y2 z2) (first vss)) (rest vss)) v2 v2?)]))])))
-    
-    (filter (compose not empty?) vss)))
+    ((make-clip-lines/plane a b c (- val) gte?) vs)))
 
-(define clip-lines-x-min (make-clip-lines 0 >= clip-line-x))
-(define clip-lines-x-max (make-clip-lines 0 <= clip-line-x))
-(define clip-lines-y-min (make-clip-lines 1 >= clip-line-y))
-(define clip-lines-y-max (make-clip-lines 1 <= clip-line-y))
-(define clip-lines-z-min (make-clip-lines 2 >= clip-line-z))
-(define clip-lines-z-max (make-clip-lines 2 <= clip-line-z))
+(define clip-lines-x-min (make-clip-lines/axis 1 0 0 >=))
+(define clip-lines-x-max (make-clip-lines/axis 1 0 0 <=))
+(define clip-lines-y-min (make-clip-lines/axis 0 1 0 >=))
+(define clip-lines-y-max (make-clip-lines/axis 0 1 0 <=))
+(define clip-lines-z-min (make-clip-lines/axis 0 0 1 >=))
+(define clip-lines-z-max (make-clip-lines/axis 0 0 1 <=))
 
-(define (clip-lines vs x-min x-max y-min y-max z-min z-max)
+(define (clip-lines/bounds vs x-min x-max y-min y-max z-min z-max)
   (let/ec return
-    ;; early reject: no polygon
-    (when (empty? vs) (return empty))
-    (match-define (list (vector xs ys zs) ...) vs)
-    ;; early accept: all endpoints in bounds
-    (when (and (andmap (λ (x) (<= x-min x x-max)) xs)
-               (andmap (λ (y) (<= y-min y y-max)) ys)
-               (andmap (λ (z) (<= z-min z z-max)) zs))
-      (return (list vs)))
-    ;; early reject: all endpoints on the outside of the same plane
-    (when (or (andmap (λ (x) (x . < . x-min)) xs) (andmap (λ (x) (x . > . x-max)) xs)
-              (andmap (λ (y) (y . < . y-min)) ys) (andmap (λ (y) (y . > . y-max)) ys)
-              (andmap (λ (z) (z . < . z-min)) zs) (andmap (λ (z) (z . > . z-max)) zs))
-      (return empty))
+    (when (early-accept? vs x-min x-max y-min y-max z-min z-max) (return (list vs)))
+    (when (early-reject? vs x-min x-max y-min y-max z-min z-max) (return empty))
     (let* ([vss  (clip-lines-x-min x-min vs)]
            [_    (when (empty? vss) (return empty))]
            [vss  (append* (map (λ (vs) (clip-lines-x-max x-max vs)) vss))]
@@ -215,3 +188,66 @@
            [_    (when (empty? vss) (return empty))]
            [vss  (append* (map (λ (vs) (clip-lines-z-max z-max vs)) vss))])
       vss)))
+
+;; ===================================================================================================
+;; Polygon clipping
+
+(define-syntax-rule (make-clip-polygon/plane a b c d gte?)
+  (λ (vs ls)
+    (define v1 (last vs))
+    (define x1 (vector-ref v1 0))
+    (define y1 (vector-ref v1 1))
+    (define z1 (vector-ref v1 2))
+    (define v1? (gte? (plane-point-dist a b c d x1 y1 z1) 0))
+    
+    (define-values (new-vs new-ls _x1 _y1 _z1 _v1?)
+      (for/fold ([vs empty] [ls empty] [x1 x1] [y1 y1] [z1 z1] [v1? v1?]) ([v2  (in-list vs)]
+                                                                           [l   (in-list ls)])
+        (define x2 (vector-ref v2 0))
+        (define y2 (vector-ref v2 1))
+        (define z2 (vector-ref v2 2))
+        (define v2? (gte? (plane-point-dist a b c d x2 y2 z2) 0))
+        (cond [(and v1? v2?)       (values (cons v2 vs) (cons l ls) x2 y2 z2 v2?)]
+              [(not (or v1? v2?))  (values vs ls x2 y2 z2 v2?)]
+              [else
+               (define-values (x y z) (plane-line-intersect a b c d x1 y1 z1 x2 y2 z2))
+               (if v1?
+                   (values (cons (vector x y z) vs) (cons l ls) x2 y2 z2 v2?)
+                   (values (list* v2 (vector x y z) vs) (list* l #t ls) x2 y2 z2 v2?))])))
+    
+    (values (reverse new-vs) (reverse new-ls))))
+
+(define (clip-polygon/plane vs ls plane)
+  (match-define (vector a b c d) plane)
+  ((make-clip-polygon/plane a b c d >=) vs ls))
+
+;; ---------------------------------------------------------------------------------------------------
+;; Clipping inside axial bounding box
+
+(define-syntax-rule (make-clip-polygon/axis a b c gte?)
+  (λ (vs ls val)
+    ((make-clip-polygon/plane a b c (- val) gte?) vs ls)))
+
+(define clip-polygon-x-min (make-clip-polygon/axis 1 0 0 >=))
+(define clip-polygon-x-max (make-clip-polygon/axis 1 0 0 <=))
+(define clip-polygon-y-min (make-clip-polygon/axis 0 1 0 >=))
+(define clip-polygon-y-max (make-clip-polygon/axis 0 1 0 <=))
+(define clip-polygon-z-min (make-clip-polygon/axis 0 0 1 >=))
+(define clip-polygon-z-max (make-clip-polygon/axis 0 0 1 <=))
+
+(define (clip-polygon/bounds vs ls x-min x-max y-min y-max z-min z-max)
+  (let/ec return
+    (when (early-accept? vs x-min x-max y-min y-max z-min z-max) (return vs ls))
+    (when (early-reject? vs x-min x-max y-min y-max z-min z-max) (return empty empty))
+    (let*-values ([(vs ls)  (clip-polygon-x-min vs ls x-min)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-x-max vs ls x-max)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-y-min vs ls y-min)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-y-max vs ls y-max)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-z-min vs ls z-min)]
+                  [(_)      (when (empty? vs) (return empty empty))]
+                  [(vs ls)  (clip-polygon-z-max vs ls z-max)])
+      (values vs ls))))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt
index 9040a1f717..66944701ec 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/contour.rkt
@@ -1,6 +1,6 @@
 #lang racket/base
 
-(require racket/class racket/match racket/list racket/flonum racket/contract
+(require racket/class racket/match racket/list racket/contract
          unstable/latent-contract/defthing
          plot/utils
          "../common/utils.rkt")
@@ -21,14 +21,11 @@
   (when (<= z-min z z-max)
     (send area put-alpha alpha)
     (send area put-pen color width style)
-    (let* ([flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-           [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-           [z           (if flonum-ok? (exact->inexact z) z)])
-      (for-2d-sample
-       (xa xb ya yb z1 z2 z3 z4) sample
-       (for ([line  (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))])
-         (match-define (list v1 v2) line)
-         (send area put-line v1 v2)))))
+    (for-2d-sample
+     (xa xb ya yb z1 z2 z3 z4) sample
+     (for ([line  (in-list (heights->lines xa xb ya yb z z1 z2 z3 z4))])
+       (match-define (list v1 v2) line)
+       (send area put-line v1 v2))))
   
   (cond [label  (line-legend-entry label color width style)]
         [else   empty]))
@@ -71,10 +68,7 @@
   (let* ([colors  (maybe-apply colors zs)]
          [widths  (maybe-apply widths zs)]
          [styles  (maybe-apply styles zs)]
-         [alphas  (maybe-apply alphas zs)]
-         [flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-         [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-         [zs          (if flonum-ok? (map exact->inexact zs) zs)])
+         [alphas  (maybe-apply alphas zs)])
     (for ([z  (in-list zs)]
           [color  (in-cycle colors)]
           [width  (in-cycle widths)]
@@ -141,10 +135,7 @@
          [alphas  (maybe-apply alphas z-ivls)]
          [line-colors  (maybe-apply line-colors z-ivls)]
          [line-widths  (maybe-apply line-widths z-ivls)]
-         [line-styles  (maybe-apply line-styles z-ivls)]
-         [flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-         [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)]
-         [zs          (if flonum-ok? (map exact->inexact zs) zs)])
+         [line-styles  (maybe-apply line-styles z-ivls)])
     (for ([za  (in-list zs)]
           [zb  (in-list (rest zs))]
           [color  (in-cycle colors)]
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt
index c08993af1e..7739d32c1a 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/isosurface.rkt
@@ -2,7 +2,6 @@
 
 (require racket/class racket/match racket/list racket/flonum racket/contract racket/math
          unstable/latent-contract/defthing
-         unstable/flonum
          plot/utils)
 
 (provide (all-defined-out))
@@ -24,13 +23,10 @@
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (let* ([flonum-ok?  (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)]
-         [sample      (if flonum-ok? (3d-sample-exact->inexact sample) sample)]
-         [d           (if flonum-ok? (exact->inexact d) d)])
-    (for-3d-sample
-     (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
-     (for ([vs  (in-list (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))])
-       (send area put-polygon vs))))
+  (for-3d-sample
+   (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
+   (for ([vs  (in-list (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))])
+     (send area put-polygon vs)))
   
   (cond [label  (rectangle-legend-entry
                  label color style line-color line-width line-style)]
@@ -85,10 +81,7 @@
            [alphas  (maybe-apply alphas ds)]
            [line-colors  (maybe-apply line-colors ds)]
            [line-widths  (maybe-apply line-widths ds)]
-           [line-styles  (maybe-apply line-styles ds)]
-           [flonum-ok?  (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)]
-           [sample      (if flonum-ok? (3d-sample-exact->inexact sample) sample)]
-           [ds          (if flonum-ok? (map exact->inexact ds) ds)])
+           [line-styles  (maybe-apply line-styles ds)])
       (for ([d      (in-list ds)]
             [color  (in-cycle colors)]
             [style  (in-cycle styles)]
@@ -153,26 +146,24 @@
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (let* ([flonum-ok?  (flonum-ok-for-4d? x-min x-max y-min y-max z-min z-max d-min d-max)]
-         [sample      (if flonum-ok? (3d-sample-exact->inexact sample) sample)]
-         [d           (if flonum-ok? 0.0 0)])
-    (for-3d-sample
-     (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
-     (cond [(and (xb . > . 0) (ya . < . 0) (yb . > . 0))
-            (let* ([yb  -0.00001]
-                   [d3  (f xb yb za)]
-                   [d4  (f xa yb za)]
-                   [d7  (f xb yb zb)]
-                   [d8  (f xa yb zb)])
-              (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
-            (let* ([ya  0.00001]
-                   [d1  (f xa ya za)]
-                   [d2  (f xb ya za)]
-                   [d5  (f xa ya zb)]
-                   [d6  (f xb ya zb)])
-              (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))]
-           [else
-            (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)])))
+  (define d 0)
+  (for-3d-sample
+   (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
+   (cond [(and (xb . > . 0) (ya . < . 0) (yb . > . 0))
+          (let* ([yb  -0.00001]
+                 [d3  (f xb yb za)]
+                 [d4  (f xa yb za)]
+                 [d7  (f xb yb zb)]
+                 [d8  (f xa yb zb)])
+            (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
+          (let* ([ya  0.00001]
+                 [d1  (f xa ya za)]
+                 [d2  (f xb ya za)]
+                 [d5  (f xa ya zb)]
+                 [d6  (f xb ya zb)])
+            (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))]
+         [else
+          (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)]))
   
   (cond [label  (rectangle-legend-entry
                  label color style line-color line-width line-style)]
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt
deleted file mode 100644
index 5e861fe768..0000000000
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/matrix.rkt
+++ /dev/null
@@ -1,61 +0,0 @@
-#lang racket/base
-
-;; A small rotation matrix library, used to transform plot coordinates into view coordinates.
-
-(require racket/match
-         (only-in math/flonum fl)
-         racket/flonum)
-
-(provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x m3-scale
-         fl3-dot)
-
-(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (fl+ (fl+ (fl* x1 x2) (fl* y1 y2)) (fl* z1 z2)))
-
-(define (m3-apply m v)
-  (match-define (vector v1 v2 v3) m)
-  (define x (flvector-ref v 0))
-  (define y (flvector-ref v 1))
-  (define z (flvector-ref v 2))
-  (flvector (dot x y z (flvector-ref v1 0) (flvector-ref v1 1) (flvector-ref v1 2))
-            (dot x y z (flvector-ref v2 0) (flvector-ref v2 1) (flvector-ref v2 2))
-            (dot x y z (flvector-ref v3 0) (flvector-ref v3 1) (flvector-ref v3 2))))
-
-(define (m3-transpose m)
-  (match-define (vector v1 v2 v3) m)
-  (vector (flvector (flvector-ref v1 0) (flvector-ref v2 0) (flvector-ref v3 0))
-          (flvector (flvector-ref v1 1) (flvector-ref v2 1) (flvector-ref v3 1))
-          (flvector (flvector-ref v1 2) (flvector-ref v2 2) (flvector-ref v3 2))))
-
-(define (m3* m1 m2)
-  (match-define (vector v1 v2 v3) m1)
-  (define m (m3-transpose m2))
-  (vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3)))
-
-(define (m3-rotate-z theta)
-  (let ([theta  (fl theta)])
-    (define cos-theta (flcos theta))
-    (define sin-theta (flsin theta))
-    (vector (flvector cos-theta (- sin-theta) 0.0)
-            (flvector sin-theta cos-theta 0.0)
-            (flvector 0.0 0.0 1.0))))
-
-(define (m3-rotate-x rho)
-  (let ([rho  (fl rho)])
-    (define cos-rho (flcos rho))
-    (define sin-rho (flsin rho))
-    (vector (flvector 1.0 0.0 0.0)
-            (flvector 0.0 cos-rho (- sin-rho))
-            (flvector 0.0 sin-rho cos-rho))))
-
-(define (m3-scale x-scale y-scale z-scale)
-  (vector (flvector (fl x-scale) 0.0 0.0)
-          (flvector 0.0 (fl y-scale) 0.0)
-          (flvector 0.0 0.0 (fl z-scale))))
-
-(define (fl3-dot v1 v2)
-  (fl+ (fl* (flvector-ref v1 0)
-            (flvector-ref v2 0))
-       (fl+ (fl* (flvector-ref v1 1)
-                 (flvector-ref v2 1))
-            (fl* (flvector-ref v1 2)
-                 (flvector-ref v2 2)))))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt
index 3caa75f0fe..2799c7c159 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/plot-area.rkt
@@ -1,20 +1,15 @@
 #lang racket/base
 
-(require racket/class racket/match racket/list racket/math racket/contract racket/vector
-         racket/flonum
+(require racket/class racket/match racket/list racket/math racket/flonum
          (only-in math fl flvector->vector vector->flvector)
-         racket/promise
-         unstable/flonum
          "../common/math.rkt"
          "../common/plot-device.rkt"
          "../common/ticks.rkt"
          "../common/draw.rkt"
-         "../common/contract.rkt"
          "../common/axis-transform.rkt"
          "../common/parameters.rkt"
-         "../common/sample.rkt"
          "../common/utils.rkt"
-         "matrix.rkt"
+         "vector.rkt"
          "clip.rkt"
          "bsp-trees.rkt"
          "bsp.rkt")
@@ -104,7 +99,8 @@
     (define/public (clear-clip-rect) (set! clipping? #f))
     
     (define (in-bounds? v)
-      (or (not clipping?) (point-in-bounds? v clip-x-min clip-x-max
+      (or (not clipping?) (point-in-bounds? v
+                                            clip-x-min clip-x-max
                                             clip-y-min clip-y-max
                                             clip-z-min clip-z-max)))
     
@@ -132,7 +128,7 @@
     ;; There are four coordinate systems:
     ;;  1. Plot coordinates (original, user-facing coordinate system)
     ;;  2. Normalized coordinates (from plot coordinates: for each axis: transform, center, and scale
-    ;;     to [-1/2,1/2]) - these are always vectors of flonum
+    ;;     to [-0.5,0.5]) - these are always flvectors
     ;;  3. View coordinates (from normalized coordinates: rotate)
     ;;  4. Device context coordinates (from view coordinates: project to 2D)
     
@@ -145,40 +141,22 @@
            (equal? (plot-y-transform) id-transform)
            (equal? (plot-z-transform) id-transform)))
     
-    (define flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max))
-    
     (define plot->norm
-      (if flonum-ok?
-          (let-map
-           (x-mid y-mid z-mid x-size y-size z-size) fl
-           (if identity-transforms?
-               (match-lambda
-                 [(vector x y z)
-                  (flvector (fl/ (fl- (fl x) x-mid) x-size)
-                            (fl/ (fl- (fl y) y-mid) y-size)
-                            (fl/ (fl- (fl z) z-mid) z-size))])
-               (match-lambda
-                 [(vector (? rational? x) (? rational? y) (? rational? z))
-                  (flvector (fl/ (fl- (fl (fx x)) x-mid) x-size)
-                            (fl/ (fl- (fl (fy y)) y-mid) y-size)
-                            (fl/ (fl- (fl (fz z)) z-mid) z-size))]
-                 [(vector x y z)
-                  (flvector +nan.0 +nan.0 +nan.0)])))
-          (if identity-transforms?
-              (match-lambda
-                [(vector (? rational? x) (? rational? y) (? rational? z))
-                 (flvector (fl (/ (- (inexact->exact x) x-mid) x-size))
-                           (fl (/ (- (inexact->exact y) y-mid) y-size))
-                           (fl (/ (- (inexact->exact z) z-mid) z-size)))]
-                [(vector x y z)
-                 (flvector +nan.0 +nan.0 +nan.0)])
-              (match-lambda
-                [(vector (? rational? x) (? rational? y) (? rational? z))
-                 (flvector (fl (/ (- (inexact->exact (fx x)) x-mid) x-size))
-                           (fl (/ (- (inexact->exact (fy y)) y-mid) y-size))
-                           (fl (/ (- (inexact->exact (fz z)) z-mid) z-size)))]
-                [(vector x y z)
-                 (flvector +nan.0 +nan.0 +nan.0)]))))
+      (if identity-transforms?
+          (match-lambda
+            [(vector (? rational? x) (? rational? y) (? rational? z))
+             (flvector (fl (/ (- x x-mid) x-size))
+                       (fl (/ (- y y-mid) y-size))
+                       (fl (/ (- z z-mid) z-size)))]
+            [(vector x y z)
+             (flvector +nan.0 +nan.0 +nan.0)])
+          (match-lambda
+            [(vector (? rational? x) (? rational? y) (? rational? z))
+             (flvector (fl (/ (- (inexact->exact (fx x)) x-mid) x-size))
+                       (fl (/ (- (inexact->exact (fy y)) y-mid) y-size))
+                       (fl (/ (- (inexact->exact (fz z)) z-mid) z-size)))]
+            [(vector x y z)
+             (flvector +nan.0 +nan.0 +nan.0)])))
     
     (define rotate-theta-matrix (m3-rotate-z theta))
     (define rotate-rho-matrix (m3-rotate-x rho))
@@ -789,11 +767,11 @@
            ;; Diffuse lighting: typical Lambertian surface model (using absolute value because we
            ;; can't expect surface normals to point the right direction)
            (define diff
-             (cond [diffuse-light?  (flabs (fl3-dot normal light-dir))]
+             (cond [diffuse-light?  (flabs (flv3-dot normal light-dir))]
                    [else  1.0]))
            ;; Specular highlighting: Blinn-Phong model
            (define spec
-             (cond [specular-light?  (fl* 32.0 (expt (flabs (fl3-dot normal half-dir)) 20.0))]
+             (cond [specular-light?  (fl* 32.0 (expt (flabs (flv3-dot normal half-dir)) 20.0))]
                    [else  0.0]))
            ;; Blend ambient light with diffuse light, return specular as it is
            ;; As ambient-light -> 1.0, contribution of diffuse -> 0.0
@@ -809,7 +787,7 @@
                           vs ls normal)
         s)
       (define view-normal (norm->view normal))
-      (define cos-view (fl3-dot view-dir view-normal))
+      (define cos-view (flv3-dot view-dir view-normal))
       (cond
         [(and (cos-view . < . 0.0) (eq? face 'front))  (void)]
         [(and (cos-view . > . 0.0) (eq? face 'back))   (void)]
@@ -819,7 +797,7 @@
          (let ([pen-color    (map (λ (v) (+ (* v diff) spec)) pen-color)]
                [brush-color  (map (λ (v) (+ (* v diff) spec)) brush-color)]
                [vs  (map norm->dc vs)])
-           ;(send pd set-pen "black" 0.5 'solid)
+           ;(send pd set-pen "black" 0.5 'solid)  ; for BSP debugging
            (send pd set-pen "black" 0 'transparent)
            (send pd set-brush brush-color brush-style)
            (send pd draw-polygon vs)
@@ -990,79 +968,79 @@
     ;; Drawing shapes
     
     (define/public (put-line v1 v2)
-      (let/ec return
-        (unless (and (vrational? v1) (vrational? v2)) (return (void)))
-        (let-values ([(v1 v2)  (if clipping?
-                                   (clip-line v1 v2 clip-x-min clip-x-max
-                                              clip-y-min clip-y-max
-                                              clip-z-min clip-z-max)
-                                   (values v1 v2))])
-          (unless (and v1 v2) (return (void)))
-          (cond [identity-transforms?
-                 (add-shape! plot3d-area-layer
-                             (line (line-data alpha pen-color pen-width pen-style)
-                                   (plot->norm v1)
-                                   (plot->norm v2)))]
-                [else
-                 (define vs (subdivide-line plot->dc v1 v2))
-                 (add-shape! plot3d-area-layer
-                             (lines (line-data alpha pen-color pen-width pen-style)
-                                    (map plot->norm vs)))]))))
-    
-    (define/public (put-lines vs)
-      (for ([vs  (vrational-sublists vs)])
-        (unless (empty? vs)
-          (let ([vss  (if clipping?
-                          (clip-lines vs clip-x-min clip-x-max
-                                      clip-y-min clip-y-max
-                                      clip-z-min clip-z-max)
-                          (list vs))])
-            (cond [identity-transforms?
-                   (for ([vs  (in-list vss)])
+      (let ([v1  (exact-vector3d v1)]
+            [v2  (exact-vector3d v2)])
+        (when (and v1 v2)
+          (let-values ([(v1 v2)  (if clipping?
+                                     (clip-line/bounds v1 v2
+                                                       clip-x-min clip-x-max
+                                                       clip-y-min clip-y-max
+                                                       clip-z-min clip-z-max)
+                                     (values v1 v2))])
+            (when (and v1 v2)
+              (cond [identity-transforms?
+                     (add-shape! plot3d-area-layer
+                                 (line (line-data alpha pen-color pen-width pen-style)
+                                       (plot->norm v1)
+                                       (plot->norm v2)))]
+                    [else
+                     (define vs (subdivide-line plot->dc v1 v2))
                      (add-shape! plot3d-area-layer
                                  (lines (line-data alpha pen-color pen-width pen-style)
-                                        (map plot->norm vs))))]
-                  [else
-                   (for ([vs  (in-list vss)])
-                     (let ([vs  (subdivide-lines plot->dc vs)])
-                       (add-shape! plot3d-area-layer
-                                   (lines (line-data alpha pen-color pen-width pen-style)
-                                          (map plot->norm vs)))))])))))
+                                        (map plot->norm vs)))]))))))
     
-    (define (add-polygon! vs face ls)
-      (let/ec return
-        (when (or (empty? vs) (not (andmap vrational? vs)))
-          (return empty))
-        
-        (define norm-vs (map (λ (v) (flvector->vector (plot->norm v))) vs))
-        (define normal (vector->flvector (vnormal norm-vs)))
-        (define center (vector->flvector (vcenter norm-vs)))
-        (let*-values
-            ([(vs ls)  (if clipping?
-                           (clip-polygon vs ls
-                                         clip-x-min clip-x-max
-                                         clip-y-min clip-y-max
-                                         clip-z-min clip-z-max)
-                           (values vs ls))]
-             [(vs ls)  (if identity-transforms?
-                           (values vs ls)
-                           (subdivide-polygon plot->dc vs ls))])
-          (when (empty? vs) (return empty))
-          (add-shape! plot3d-area-layer
-                      (poly (poly-data alpha center
-                                       pen-color pen-width pen-style
-                                       brush-color brush-style face)
-                            (map plot->norm vs) ls normal)))))
+    (define/public (put-lines vs)
+      (for ([vs  (in-list (exact-vector3d-sublists vs))])
+        (let ([vss  (if clipping?
+                        (clip-lines/bounds vs
+                                           clip-x-min clip-x-max
+                                           clip-y-min clip-y-max
+                                           clip-z-min clip-z-max)
+                        (list vs))])
+          (cond [identity-transforms?
+                 (for ([vs  (in-list vss)])
+                   (add-shape! plot3d-area-layer
+                               (lines (line-data alpha pen-color pen-width pen-style)
+                                      (map plot->norm vs))))]
+                [else
+                 (for ([vs  (in-list vss)])
+                   (let ([vs  (subdivide-lines plot->dc vs)])
+                     (add-shape! plot3d-area-layer
+                                 (lines (line-data alpha pen-color pen-width pen-style)
+                                        (map plot->norm vs)))))]))))
     
     (define/public (put-polygon vs [face 'both] [ls (make-list (length vs) #t)])
-      (add-polygon! vs face ls))
+      (let-values ([(vs ls)  (exact-polygon3d vs ls)])
+        (unless (empty? vs)
+          (let*-values ([(vs ls)  (if clipping?
+                                      (clip-polygon/bounds vs ls
+                                                           clip-x-min clip-x-max
+                                                           clip-y-min clip-y-max
+                                                           clip-z-min clip-z-max)
+                                      (values vs ls))]
+                        [(vs ls)  (if identity-transforms?
+                                      (values vs ls)
+                                      (subdivide-polygon plot->dc vs ls))])
+            (unless (empty? vs)
+              (define norm-vs (map plot->norm vs))
+              (define normal (flv3-normal norm-vs))
+              (define center (flv3-center norm-vs))
+              (add-shape! plot3d-area-layer
+                          (poly (poly-data alpha center
+                                           pen-color pen-width pen-style
+                                           brush-color brush-style face)
+                                norm-vs ls normal)))))))
     
     (define/public (put-rect r)
       (when (rect-rational? r)
         (let ([r  (rect-meet r bounds-rect)])
           (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r)
-          (define v-min (plot->norm (vector x-min y-min z-min)))
-          (define v-max (plot->norm (vector x-max y-max z-max)))
+          (define v-min (plot->norm (vector (inexact->exact x-min)
+                                            (inexact->exact y-min)
+                                            (inexact->exact z-min))))
+          (define v-max (plot->norm (vector (inexact->exact x-max)
+                                            (inexact->exact y-max)
+                                            (inexact->exact z-max))))
           (let ()
             (define x-min (flvector-ref v-min 0))
             (define y-min (flvector-ref v-min 1))
@@ -1112,13 +1090,14 @@
     (define/public (put-text str v [anchor 'center] [angle 0] [dist 0]
                              #:outline? [outline? #f]
                              #:layer [layer plot3d-area-layer])
-      (when (and (vrational? v) (in-bounds? v))
-        (add-shape! layer (points (text-data alpha anchor angle dist str
-                                             font-size font-family text-foreground outline?)
-                                  (list (plot->norm v))))))
+      (let ([v  (exact-vector3d v)])
+        (when (and v (in-bounds? v))
+          (add-shape! layer (points (text-data alpha anchor angle dist str
+                                               font-size font-family text-foreground outline?)
+                                    (list (plot->norm v)))))))
     
     (define/public (put-glyphs vs symbol size #:layer [layer plot3d-area-layer])
-      (let ([vs  (filter (λ (v) (and (vrational? v) (in-bounds? v))) vs)])
+      (let ([vs  (filter (λ (v) (and v (in-bounds? v))) (map exact-vector3d vs))])
         (unless (empty? vs)
           (add-shape! layer (points (glyph-data alpha symbol size
                                                 pen-color pen-width pen-style
@@ -1126,14 +1105,16 @@
                                     (map plot->norm vs))))))
     
     (define/public (put-arrow v1 v2)
-      (when (and (vrational? v1) (vrational? v2) (in-bounds? v1))
-        (cond [(in-bounds? v2)
-               (define c (v* (v+ v1 v2) 1/2))
-               (define outline-color (->brush-color (plot-background)))
-               (add-shape! plot3d-area-layer
-                           (points (arrow-data alpha (plot->norm v1) (plot->norm v2)
-                                               outline-color pen-color pen-width pen-style)
-                                   (list (plot->norm c))))]
-              [else
-               (put-line v1 v2)])))
+      (let ([v1  (exact-vector3d v1)]
+            [v2  (exact-vector3d v2)])
+        (when (and v1 v2 (in-bounds? v1))
+          (cond [(in-bounds? v2)
+                 (define c (v* (v+ v1 v2) 1/2))
+                 (define outline-color (->brush-color (plot-background)))
+                 (add-shape! plot3d-area-layer
+                             (points (arrow-data alpha (plot->norm v1) (plot->norm v2)
+                                                 outline-color pen-color pen-width pen-style)
+                                     (list (plot->norm c))))]
+                [else
+                 (put-line v1 v2)]))))
     )) ; end class
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt
index fc42973205..a785cc1a85 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/split.rkt
@@ -1,7 +1,6 @@
 #lang typed/racket/base
 
-(require racket/match
-         racket/list
+(require racket/list
          racket/flonum)
 
 (provide point3d-plane-dist
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt
index 7b0bf5469d..ca00ac0190 100644
--- a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/surface.rkt
@@ -1,6 +1,6 @@
 #lang racket/base
 
-(require racket/class racket/match racket/list racket/flonum racket/contract
+(require racket/class racket/match racket/list racket/contract
          unstable/latent-contract/defthing
          plot/utils)
 
@@ -22,12 +22,10 @@
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (let* ([flonum-ok?  (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max)]
-         [sample      (if flonum-ok? (2d-sample-exact->inexact sample) sample)])
-    (for-2d-sample
-     (xa xb ya yb z1 z2 z3 z4) sample
-     (define vs (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4)))
-     (send area put-polygon vs)))
+  (for-2d-sample
+   (xa xb ya yb z1 z2 z3 z4) sample
+   (define vs (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4)))
+   (send area put-polygon vs))
   
   (cond [label  (rectangle-legend-entry label color style line-color line-width line-style)]
         [else   empty]))
diff --git a/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt
new file mode 100644
index 0000000000..026a1b9f5b
--- /dev/null
+++ b/pkgs/plot-pkgs/plot-lib/plot/private/plot3d/vector.rkt
@@ -0,0 +1,144 @@
+#lang racket/base
+
+(require racket/match
+         racket/list
+         (only-in math/flonum fl flvector+ flvector-)
+         (only-in racket/unsafe/ops unsafe-vector-ref)
+         racket/flonum)
+
+(provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x
+         flv3-dot
+         flv3-normal
+         flv3-center
+         exact-vector3d
+         exact-vector3d-sublists
+         exact-polygon3d)
+
+(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (fl+ (fl+ (fl* x1 x2) (fl* y1 y2)) (fl* z1 z2)))
+
+(define (m3-apply m v)
+  (match-define (vector v1 v2 v3) m)
+  (define x (flvector-ref v 0))
+  (define y (flvector-ref v 1))
+  (define z (flvector-ref v 2))
+  (flvector (dot x y z (flvector-ref v1 0) (flvector-ref v1 1) (flvector-ref v1 2))
+            (dot x y z (flvector-ref v2 0) (flvector-ref v2 1) (flvector-ref v2 2))
+            (dot x y z (flvector-ref v3 0) (flvector-ref v3 1) (flvector-ref v3 2))))
+
+(define (m3-transpose m)
+  (match-define (vector v1 v2 v3) m)
+  (vector (flvector (flvector-ref v1 0) (flvector-ref v2 0) (flvector-ref v3 0))
+          (flvector (flvector-ref v1 1) (flvector-ref v2 1) (flvector-ref v3 1))
+          (flvector (flvector-ref v1 2) (flvector-ref v2 2) (flvector-ref v3 2))))
+
+(define (m3* m1 m2)
+  (match-define (vector v1 v2 v3) m1)
+  (define m (m3-transpose m2))
+  (vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3)))
+
+(define (m3-rotate-z theta)
+  (let ([theta  (fl theta)])
+    (define cos-theta (flcos theta))
+    (define sin-theta (flsin theta))
+    (vector (flvector cos-theta (- sin-theta) 0.0)
+            (flvector sin-theta cos-theta 0.0)
+            (flvector 0.0 0.0 1.0))))
+
+(define (m3-rotate-x rho)
+  (let ([rho  (fl rho)])
+    (define cos-rho (flcos rho))
+    (define sin-rho (flsin rho))
+    (vector (flvector 1.0 0.0 0.0)
+            (flvector 0.0 cos-rho (- sin-rho))
+            (flvector 0.0 sin-rho cos-rho))))
+
+(define (flv3-dot v1 v2)
+  (fl+ (fl* (flvector-ref v1 0)
+            (flvector-ref v2 0))
+       (fl+ (fl* (flvector-ref v1 1)
+                 (flvector-ref v2 1))
+            (fl* (flvector-ref v1 2)
+                 (flvector-ref v2 2)))))
+
+(define default-normal (flvector 0.0 -1.0 0.0))
+
+(define (flv3-normal vs)
+  (define n (length vs))
+  (cond
+    [(n . < . 3)  default-normal]
+    [else
+     (match-define (list v1 v2) (take-right vs 2))
+     (define x1 (flvector-ref v1 0))
+     (define y1 (flvector-ref v1 1))
+     (define z1 (flvector-ref v1 2))
+     (define x2 (flvector-ref v2 0))
+     (define y2 (flvector-ref v2 1))
+     (define z2 (flvector-ref v2 2))
+     (define-values (x y z _x1 _y1 _z1 _x2 _y2 _z2)
+       (for/fold ([x 0.0] [y 0.0] [z 0.0] [x1 x1] [y1 y1] [z1 z1] [x2 x2] [y2 y2] [z2 z2]
+                          ) ([v3  (in-list vs)])
+         (define x3 (flvector-ref v3 0))
+         (define y3 (flvector-ref v3 1))
+         (define z3 (flvector-ref v3 2))
+         (define x32 (fl- x3 x2))
+         (define y32 (fl- y3 y2))
+         (define z32 (fl- z3 z2))
+         (define x12 (fl- x1 x2))
+         (define y12 (fl- y1 y2))
+         (define z12 (fl- z1 z2))
+         (values (+ x (fl- (fl* y32 z12) (fl* z32 y12)))
+                 (+ y (fl- (fl* z32 x12) (fl* x32 z12)))
+                 (+ z (fl- (fl* x32 y12) (fl* y32 x12)))
+                 x2 y2 z2
+                 x3 y3 z3)))
+     (define m (flsqrt (fl+ (fl* x x) (fl+ (fl* y y) (fl* z z)))))
+     (if (fl> m 0.0)
+         (flvector (fl/ x m) (fl/ y m) (fl/ z m))
+         default-normal)]))
+
+(define (flv3-center vs)
+  (define xs (map (λ (v) (flvector-ref v 0)) vs))
+  (define ys (map (λ (v) (flvector-ref v 1)) vs))
+  (define zs (map (λ (v) (flvector-ref v 2)) vs))
+  (flvector (* 0.5 (+ (apply min xs) (apply max xs)))
+            (* 0.5 (+ (apply min ys) (apply max ys)))
+            (* 0.5 (+ (apply min zs) (apply max zs)))))
+
+(define (exact-vector3d v)
+  (define n (vector-length v))
+  (cond [(= n 3)
+         (define v1 (unsafe-vector-ref v 0))
+         (define v2 (unsafe-vector-ref v 1))
+         (define v3 (unsafe-vector-ref v 2))
+         (cond [(and (exact? v1) (exact? v2) (exact? v3))  v]
+               [(and (rational? v1) (rational? v2) (rational? v3))
+                (vector (inexact->exact v1) (inexact->exact v2) (inexact->exact v3))]
+               [else  #f])]
+        [else  #f]))
+
+(define (sublists vs)
+  (define vss
+    (for/fold ([vss  (list null)]) ([v  (in-list vs)])
+      (cond [v  (cons (cons v (car vss)) (cdr vss))]
+            [(null? (car vss))  vss]
+            [else  (cons null vss)])))
+  (cond [(null? (car vss))  (cdr vss)]
+        [else  vss]))
+
+(define (exact-vector3d-sublists vs)
+  (sublists (map exact-vector3d vs)))
+
+(define (exact-polygon3d vs ls)
+  (cond
+    [(null? vs)  (values null null)]
+    [else
+     (define-values (new-vs new-ls _)
+       (for/fold ([vs null] [ls null] [v1 (last vs)]) ([v2  (in-list vs)]
+                                                       [l   (in-list ls)])
+         (cond [(equal? v1 v2)  (values vs ls v2)]
+               [else
+                (define exact-v2 (exact-vector3d v2))
+                (if exact-v2
+                    (values (cons exact-v2 vs) (cons l ls) v2)
+                    (values vs ls v2))])))
+     (values (reverse new-vs) (reverse new-ls))]))
diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt
index 0d8bfe2e5d..2f2c1432ef 100644
--- a/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt
+++ b/pkgs/plot-pkgs/plot-test/plot/tests/extreme-bounds-tests.rkt
@@ -6,18 +6,20 @@
   (module config info
     (define timeout 150)))
 
+(printf "Point at 0,0:~n")
 (plot (points '(#(0 0)))
       #:x-min +min.0 #:x-max (flstep +min.0 1000)
       #:y-min 0 #:y-max 1)
 
+(printf "Point at 0,0,0:~n")
 (plot3d (points3d '(#(0 0 0)))
       #:x-min +min.0 #:x-max (flstep +min.0 1000)
       #:y-min 0 #:y-max 1 #:z-min 0 #:z-max 1)
 
-;; this should go all the way from the bottom to the top
+(printf "Steps should appear all the way from the bottom to the top:~n")
 (plot (function exact->inexact (flstep 0.0 -1) (flstep 0.0 2)))
 
-;; this should go all the way from the bottom to the top
+(printf "Steps should appear all the way from the bottom to the top:~n")
 (plot3d (surface3d (λ (x y) (exact->inexact x))
                    (flstep 0.0 -1) (flstep 0.0 2)
                    (flstep 0.0 -1) (flstep 0.0 2)))
diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt
index fa831585f5..5c535bb5c6 100644
--- a/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt
+++ b/pkgs/plot-pkgs/plot-test/plot/tests/plot2d-tests.rkt
@@ -89,6 +89,7 @@
  (plot (list (points (map vector ys xs) #:sym 'full6star #:color 1)
              (points (map vector xs ys) #:sym "B" #:color 3))))
 
+(printf "This plot should be blank:~n")
 (time
  (plot (vector-field (λ (x y) (vector +nan.0 +nan.0)))
        #:x-min -2 #:x-max 2 #:y-min -2 #:y-max 2))
diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt
index be8ed2a60c..4fbe061bc9 100644
--- a/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt
+++ b/pkgs/plot-pkgs/plot-test/plot/tests/plot3d-bsp-tests.rkt
@@ -12,18 +12,10 @@
   (plot-y-transform cbrt-transform)
   (plot-z-transform cbrt-transform))
 
-#;
-(plot3d (contour-intervals3d * -1 1 -1 1 #:samples 3 #:line-widths '(2)))
-
-(plot3d (list (surface3d * -1 1 -1 1 #:samples 2)
+(plot3d (list (contour-intervals3d * -1 1 -1 1 #:samples 2)
               (isosurface3d (λ (x y z) z)
                               0 -1 1 -1 1 -1 1
-                              #:samples 2 #:color 3)
-              #;
-              (isosurface3d (λ (x y z) (+ x y -0.5))
-                              0 -1 1 -1 1 -1 1
-                              #:samples 2 #:color 2)
-              ))
+                              #:samples 2 #:color 3)))
 
 (plot3d (list (parametric3d (λ (t)
                               (list (* 0.5 (cos (* 2 pi (- t))))
@@ -42,8 +34,7 @@
                             #:samples 53
                             #:color "red"
                             #:width 15
-                            #:alpha 0.6)
-              ))
+                            #:alpha 0.6)))
 
 (plot3d (list (isosurface3d (λ (x y z) (+ x y z)) 0 -1 1 -1 1 -1 1
                             #:samples 2
@@ -114,6 +105,10 @@
     ;#:out-file "test.pdf"
     )))
 
+(define (f2 x y)
+  (let ([x  (fl x)] [y  (fl y)])
+    (- (sqrt (+ (abs y) (abs x))))))
+
 (plot3d (list (surface3d * -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 1)
               (surface3d (λ (x y) (+ 0.1 (* x y))) -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 2)
               (surface3d (λ (x y) (+ 0.2 (* x y))) -1 1 -1 1 #:samples 6 #:alpha 0.75 #:color 3
diff --git a/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt b/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt
index d355e850b2..fe864caebd 100644
--- a/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt
+++ b/pkgs/plot-pkgs/plot-test/plot/tests/typed/parameter-regression-tests.rkt
@@ -14,7 +14,8 @@
                         (list (ivl 1 2) (ivl 1 2)))))
 
 ;; ---------------------------------------------------------------------------------------------------
-;; The following two plots should be the same
+
+(printf "The following two plots should be the same:~n")
 
 (contour-colors (λ (n) '((255 0 0) (0 0 255))))
 (contour-widths (λ (n) '(.25 .5 .75)))
@@ -33,7 +34,8 @@
 (plot (contour-intervals + 0 1 0 1))
 
 ;; ---------------------------------------------------------------------------------------------------
-;; The following two plots should be the same
+
+(printf "The following two plots should be the same:~n")
 
 (stacked-histogram-alphas (λ (n) '(0.25 0.5 0.75)))
 (stacked-histogram-colors (λ (n) '(1 2 3)))
@@ -50,7 +52,8 @@
                                (list 'c '()) (list 'd '(1/2)))))
 
 ;; ---------------------------------------------------------------------------------------------------
-;; The following two plots should be the same
+
+(printf "The following two plots should be the same:~n")
 
 (contour-interval-line-colors (λ (n) '(1 2 3)))
 (contour-interval-line-widths (λ (n) '(0.25 0.5 0.75)))
@@ -63,7 +66,8 @@
 (plot3d (contour-intervals3d + 0 1 0 1))
 
 ;; ---------------------------------------------------------------------------------------------------
-;; The following two plots should be the same
+
+(printf "The following two plots should be the same:~n")
 
 (isosurface-colors (λ (n) '(1 2 3)))
 (isosurface-styles (λ (n) '(0 1 2)))