diff --git a/collects/plot/common/draw.rkt b/collects/plot/common/draw.rkt
index a2576f3e25..2f75da637c 100644
--- a/collects/plot/common/draw.rkt
+++ b/collects/plot/common/draw.rkt
@@ -279,8 +279,6 @@
 ;; As long as this process is monotone and bounded, the distance off the dc is zero in the limit. In
 ;; practice, only a few iterations drives this distance to less than 1 drawing unit.
 
-(define (appx= x y) ((abs (- x y)) . < . 1/2))
-
 (define (margin-fixpoint x-min x-max y-min y-max
                          init-left init-right init-top init-bottom
                          get-vs)
@@ -288,21 +286,20 @@
     (for/fold ([left init-left] [right init-right] [top init-top] [bottom init-bottom]
                                 ) ([i  (in-range 3)])
       (match-define (list (vector xs ys) ...) (get-vs left right top bottom))
-      
       (define param-x-min (apply min x-min xs))
       (define param-x-max (apply max (sub1 x-max) xs))
       (define param-y-min (apply min y-min ys))
       (define param-y-max (apply max (sub1 y-max) ys))
       
-      (define new-left (+ left (- x-min param-x-min)))
-      (define new-right (- right (- (sub1 x-max) param-x-max)))
-      (define new-top (+ top (- y-min param-y-min)))
-      (define new-bottom (- bottom (- (sub1 y-max) param-y-max)))
+      (define new-left (round (+ left (- x-min param-x-min))))
+      (define new-right (round (- right (- (sub1 x-max) param-x-max))))
+      (define new-top (round (+ top (- y-min param-y-min))))
+      (define new-bottom (round (- bottom (- (sub1 y-max) param-y-max))))
       
       ;; Early out: if the margins haven't changed much, another iteration won't change them more
       ;; (hopefully)
-      (when (and (appx= left new-left) (appx= right new-right)
-                 (appx= top new-top) (appx= bottom new-bottom))
+      (when (and (= left new-left) (= right new-right)
+                 (= top new-top) (= bottom new-bottom))
         (return new-left new-right new-top new-bottom))
       
       (values new-left new-right new-top new-bottom))))
@@ -380,30 +377,30 @@
   (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r)
   (list 
    ;; Top (z-max) face
-   (list (vector 0 0 1)
+   (list (vector 0.0 0.0 1.0)
          (vector x-min y-min z-max) (vector x-max y-min z-max)
          (vector x-max y-max z-max) (vector x-min y-max z-max))
    ;; Front (y-min) face
-   (if ((cos theta) . > . 0)
-       (list (vector 0 -1 0)
+   (if ((cos theta) . > . 0.0)
+       (list (vector 0.0 -1.0 0.0)
              (vector x-min y-min z-min) (vector x-max y-min z-min)
              (vector x-max y-min z-max) (vector x-min y-min z-max))
        empty)
    ;; Back (y-max) face
-   (if ((cos theta) . < . 0)
-       (list (vector 0 1 0)
+   (if ((cos theta) . < . 0.0)
+       (list (vector 0.0 1.0 0.0)
              (vector x-min y-max z-min) (vector x-max y-max z-min)
              (vector x-max y-max z-max) (vector x-min y-max z-max))
        empty)
    ;; Left (x-min) face
-   (if ((sin theta) . > . 0)
-       (list (vector -1 0 0)
+   (if ((sin theta) . > . 0.0)
+       (list (vector -1.0 0.0 0.0)
              (vector x-min y-min z-min) (vector x-min y-max z-min)
              (vector x-min y-max z-max) (vector x-min y-min z-max))
        empty)
    ;; Right (x-max) face
-   (if ((sin theta) . < . 0)
-       (list (vector 1 0 0)
+   (if ((sin theta) . < . 0.0)
+       (list (vector 1.0 0.0 0.0)
              (vector x-max y-min z-min) (vector x-max y-max z-min)
              (vector x-max y-max z-max) (vector x-max y-min z-max))
        empty)))
diff --git a/collects/plot/common/marching-cubes.rkt b/collects/plot/common/marching-cubes.rkt
index ca1d64a38f..e9d8c4536e 100644
--- a/collects/plot/common/marching-cubes.rkt
+++ b/collects/plot/common/marching-cubes.rkt
@@ -3,6 +3,7 @@
 (require racket/contract racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops
          (for-syntax racket/base racket/syntax racket/match racket/list)
          "math.rkt"
+         "utils.rkt"
          "marching-utils.rkt"
          "contract-doc.rkt")
 
@@ -16,20 +17,20 @@
 
 ;; edge vertexes 
 
-(define-syntax-rule (edge-1-2 d d1 d2) (vector (solve-t d d1 d2) 0.0 0.0))
-(define-syntax-rule (edge-2-3 d d2 d3) (vector 1.0 (solve-t d d2 d3) 0.0))
-(define-syntax-rule (edge-3-4 d d3 d4) (vector (solve-t d d4 d3) 1.0 0.0))
-(define-syntax-rule (edge-1-4 d d1 d4) (vector 0.0 (solve-t d d1 d4) 0.0))
+(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-5-6 d d5 d6) (vector (solve-t d d5 d6) 0.0 1.0))
-(define-syntax-rule (edge-6-7 d d6 d7) (vector 1.0 (solve-t d d6 d7) 1.0))
-(define-syntax-rule (edge-7-8 d d7 d8) (vector (solve-t d d7 d8) 1.0 1.0))
-(define-syntax-rule (edge-5-8 d d5 d8) (vector 0.0 (solve-t d d5 d8) 1.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-1-5 d d1 d5) (vector 0.0 0.0 (solve-t d d1 d5)))
-(define-syntax-rule (edge-2-6 d d2 d6) (vector 1.0 0.0 (solve-t d d2 d6)))
-(define-syntax-rule (edge-3-7 d d3 d7) (vector 1.0 1.0 (solve-t d d3 d7)))
-(define-syntax-rule (edge-4-8 d d4 d8) (vector 0.0 1.0 (solve-t d d4 d8)))
+(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)))
 
 #|
 Cube vertex numbers:
@@ -74,10 +75,10 @@ Cube vertex numbers:
   (define da (unsafe-flavg4 d1 d2 d3 d4))
   (cond
     [(test? da d)
-      (list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3)
-                  (edge-3-7 d d3 d7) (edge-1-5 d d1 d5))
-            (list (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)
-                  (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)))]
+     (list (list (edge-1-2 d d1 d2) (edge-2-3 d d2 d3)
+                 (edge-3-7 d d3 d7) (edge-1-5 d d1 d5))
+           (list (edge-3-4 d d3 d4) (edge-1-4 d d1 d4)
+                 (edge-1-5 d d1 d5) (edge-3-7 d d3 d7)))]
     [else
      (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)))]))
@@ -104,7 +105,7 @@ Cube vertex numbers:
      (list (list (edge-1-5 d d1 d5) (edge-2-6 d d2 d6)
                  (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<))
 
@@ -510,32 +511,62 @@ 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 facet-num
+    (add-digit
+     j1 (add-digit
+         j2 (add-digit
+             j3 (add-digit
+                 j4 (add-digit
+                     j5 (add-digit
+                         j6 (add-digit j7 j8))))))))
+  (define f (vector-ref cube-dispatch-table facet-num))
+  (f d d1 d2 d3 d4 d5 d6 d7 d8))
+
 (defproc (heights->cube-polys [xa real?] [xb real?] [ya real?] [yb real?] [za real?] [zb real?]
                               [d real?]
                               [d1 real?] [d2 real?] [d3 real?] [d4 real?]
                               [d5 real?] [d6 real?] [d7 real?] [d8 real?]
                               ) (listof (vector/c real? real? real?))
-  (check-all-real! 'heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)
-  (let-exact->inexact
-   (xa xb ya yb za zb 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 facet-num
-     (add-digit
-      j1 (add-digit
-          j2 (add-digit
-              j3 (add-digit
-                  j4 (add-digit
-                      j5 (add-digit
-                          j6 (add-digit j7 j8))))))))
-   (define f (vector-ref cube-dispatch-table facet-num))
-   (for/list ([poly  (in-list (f d d1 d2 d3 d4 d5 d6 d7 d8))])
-     (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))))))
+  (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))
+         (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)))))]))
diff --git a/collects/plot/common/marching-squares.rkt b/collects/plot/common/marching-squares.rkt
index 3f517f7e0a..810c36e46f 100644
--- a/collects/plot/common/marching-squares.rkt
+++ b/collects/plot/common/marching-squares.rkt
@@ -3,6 +3,7 @@
 (require racket/flonum racket/fixnum racket/list racket/match racket/unsafe/ops racket/contract
          (for-syntax racket/base)
          "math.rkt"
+         "utils.rkt"
          "marching-utils.rkt"
          "contract-doc.rkt")
 
@@ -58,23 +59,23 @@ above.
 
 ;(: lines1000 (FType Lines))
 (define-syntax-rule (lines1000 z z1 z2 z3 z4)
-  (list (vector (solve-t z z1 z2) 0.0
-                0.0 (solve-t z z1 z4))))
+  (list (vector (unsafe-solve-t z z1 z2) 0.0
+                0.0 (unsafe-solve-t z z1 z4))))
 
 ;(: lines0100 (FType Lines))
 (define-syntax-rule (lines0100 z z1 z2 z3 z4)
-  (list (vector (solve-t z z1 z2) 0.0
-                1.0 (solve-t z z2 z3))))
+  (list (vector (unsafe-solve-t z z1 z2) 0.0
+                1.0 (unsafe-solve-t z z2 z3))))
 
 ;(: lines0010 (FType Lines))
 (define-syntax-rule (lines0010 z z1 z2 z3 z4)
-  (list (vector 1.0 (solve-t z z2 z3)
-                (solve-t z z4 z3) 1.0)))
+  (list (vector 1.0 (unsafe-solve-t z z2 z3)
+                (unsafe-solve-t z z4 z3) 1.0)))
 
 ;(: lines0001 (FType Lines))
 (define-syntax-rule (lines0001 z z1 z2 z3 z4)
-  (list (vector 0.0 (solve-t z z1 z4)
-                (solve-t z z4 z3) 1.0)))
+  (list (vector 0.0 (unsafe-solve-t z z1 z4)
+                (unsafe-solve-t z z4 z3) 1.0)))
 
 (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))
@@ -86,13 +87,13 @@ above.
 
 ;(: lines1100 (FType Lines))
 (define-syntax-rule (lines1100 z z1 z2 z3 z4)
-  (list (vector 0.0 (solve-t z z1 z4)
-                1.0 (solve-t z z2 z3))))
+  (list (vector 0.0 (unsafe-solve-t z z1 z4)
+                1.0 (unsafe-solve-t z z2 z3))))
 
 ;(: lines0110 (FType Lines))
 (define-syntax-rule (lines0110 z z1 z2 z3 z4)
-  (list (vector (solve-t z z1 z2) 0.0
-                (solve-t z z4 z3) 1.0)))
+  (list (vector (unsafe-solve-t z z1 z2) 0.0
+                (unsafe-solve-t z z4 z3) 1.0)))
 
 (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))
@@ -105,64 +106,83 @@ above.
   ; disambiguate using average of corners as guess for center value
   (let ([z5  (unsafe-flavg4 z1 z2 z3 z4)])
     (if (test? z5 z)
-        (list (vector (solve-t z z1 z2) 0.0
-                      1.0 (solve-t z z2 z3))
-              (vector 0.0 (solve-t z z1 z4)
-                      (solve-t z z4 z3) 1.0))
-        (list (vector (solve-t z z1 z2) 0.0
-                      0.0 (solve-t z z1 z4))
-              (vector 1.0 (solve-t z z2 z3)
-                      (solve-t z z4 z3) 1.0)))))
+        (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)))))
 
 (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 (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))
+  (if p1
+      (if p2
+          (if p3
+              (if p4
+                  (lines1111 z z1 z2 z3 z4)
+                  (lines1110 z z1 z2 z3 z4))
+              (if p4
+                  (lines1101 z z1 z2 z3 z4)
+                  (lines1100 z z1 z2 z3 z4)))
+          (if p3
+              (if p4
+                  (lines1011 z z1 z2 z3 z4)
+                  (lines1010 z z1 z2 z3 z4))
+              (if p4
+                  (lines1001 z z1 z2 z3 z4)
+                  (lines1000 z z1 z2 z3 z4))))
+      (if p2
+          (if p3
+              (if p4
+                  (lines0111 z z1 z2 z3 z4)
+                  (lines0110 z z1 z2 z3 z4))
+              (if p4
+                  (lines0101 z z1 z2 z3 z4)
+                  (lines0100 z z1 z2 z3 z4)))
+          (if p3
+              (if p4
+                  (lines0011 z z1 z2 z3 z4)
+                  (lines0010 z z1 z2 z3 z4))
+              (if p4
+                  (lines0001 z z1 z2 z3 z4)
+                  (lines0000 z z1 z2 z3 z4))))))
+
 (defproc (heights->lines [xa real?] [xb real?] [ya real?] [yb real?]
                          [z real?] [z1 real?] [z2 real?] [z3 real?] [z4 real?]
                          ) (list/c (vector/c real? real? real?) (vector/c real? real? real?))
-  (check-all-real! 'heights->lines xa xb ya yb z z1 z2 z3 z4)
-  (let-exact->inexact
-   (xa xb ya yb 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 lines
-     (if p1
-         (if p2
-             (if p3
-                 (if p4
-                     (lines1111 z z1 z2 z3 z4)
-                     (lines1110 z z1 z2 z3 z4))
-                 (if p4
-                     (lines1101 z z1 z2 z3 z4)
-                     (lines1100 z z1 z2 z3 z4)))
-             (if p3
-                 (if p4
-                     (lines1011 z z1 z2 z3 z4)
-                     (lines1010 z z1 z2 z3 z4))
-                 (if p4
-                     (lines1001 z z1 z2 z3 z4)
-                     (lines1000 z z1 z2 z3 z4))))
-         (if p2
-             (if p3
-                 (if p4
-                     (lines0111 z z1 z2 z3 z4)
-                     (lines0110 z z1 z2 z3 z4))
-                 (if p4
-                     (lines0101 z z1 z2 z3 z4)
-                     (lines0100 z z1 z2 z3 z4)))
-             (if p3
-                 (if p4
-                     (lines0011 z z1 z2 z3 z4)
-                     (lines0010 z z1 z2 z3 z4))
-                 (if p4
-                     (lines0001 z z1 z2 z3 z4)
-                     (lines0000 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)))))
+  (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))]
+        [(= 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))))]))
 
 ;; =============================================================================
 ;; Isoband marching squares: polygonizes contour between two isoline values
@@ -183,16 +203,16 @@ above.
 ;; all corners same
 
 (define (polys0000 za zb z1 z2 z3 z4) empty)
-(define (polys1111 za zb z1 z2 z3 z4) (list 'full))
 (define (polys2222 za zb z1 z2 z3 z4) empty)
+(define (polys1111 za zb z1 z2 z3 z4) (list 'full))
 
 ;; -----------------------------------------------------------------------------
 ;; single triangle
 
 (define (polys1000 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (solve-t za z1 z4) za))))
+              (vector (unsafe-solve-t za z1 z2) 0.0 za)
+              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define polys0100 (rotate-facet polys1000))
 (define polys0010 (rotate-facet polys0100))
@@ -200,8 +220,8 @@ above.
 
 (define (polys1222 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (vector (unsafe-solve-t zb z1 z2) 0.0 zb)
+              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
 
 (define polys2122 (rotate-facet polys1222))
 (define polys2212 (rotate-facet polys2122))
@@ -211,20 +231,20 @@ above.
 ;; single trapezoid
 
 (define (polys2000 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (solve-t za z1 z4) za)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+  (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))))
 
 (define polys0200 (rotate-facet polys2000))
 (define polys0020 (rotate-facet polys0200))
 (define polys0002 (rotate-facet polys0020))
 
 (define (polys0222 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) za))))
+  (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))))
 
 (define polys2022 (rotate-facet polys0222))
 (define polys2202 (rotate-facet polys2022))
@@ -236,8 +256,8 @@ above.
 (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 (solve-t za z2 z3) za)
-              (vector 0.0 (solve-t za z1 z4) za))))
+              (vector 1.0 (unsafe-solve-t za z2 z3) za)
+              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define polys0110 (rotate-facet polys1100))
 (define polys0011 (rotate-facet polys0110))
@@ -246,18 +266,18 @@ above.
 (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 (solve-t zb z2 z3) zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (vector 1.0 (unsafe-solve-t zb z2 z3) zb)
+              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
 
 (define polys2112 (rotate-facet polys1122))
 (define polys2211 (rotate-facet polys2112))
 (define polys1221 (rotate-facet polys2211))
 
 (define (polys0022 za zb z1 z2 z3 z4)
-  (list (list (vector 0.0 (solve-t za z1 z4) za)
-              (vector 1.0 (solve-t za z2 z3) za)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+  (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))))
 
 (define polys2002 (rotate-facet polys0022))
 (define polys2200 (rotate-facet polys2002))
@@ -267,22 +287,22 @@ above.
 ;; single pentagon
 
 (define (polys0111 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
+  (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 (solve-t za z1 z4) za))))
+              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define polys1011 (rotate-facet polys0111))
 (define polys1101 (rotate-facet polys1011))
 (define polys1110 (rotate-facet polys1101))
 
 (define (polys2111 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t zb z1 z2) 0.0 zb)
+  (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 (solve-t zb z1 z4) zb))))
+              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
 
 (define polys1211 (rotate-facet polys2111))
 (define polys1121 (rotate-facet polys1211))
@@ -290,10 +310,10 @@ above.
 
 (define (polys1002 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t za z4 z3) 1.0 za)
-              (vector (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (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))))
 
 (define polys2100 (rotate-facet polys1002))
 (define polys0210 (rotate-facet polys2100))
@@ -301,10 +321,10 @@ above.
 
 (define (polys1220 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector (solve-t zb z4 z3) 1.0 zb)
-              (vector (solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (solve-t za z1 z4) za))))
+              (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))))
 
 (define polys0122 (rotate-facet polys1220))
 (define polys2012 (rotate-facet polys0122))
@@ -312,10 +332,10 @@ above.
 
 (define (polys1200 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 1.0 (solve-t za z2 z3) za)
-              (vector 0.0 (solve-t za z1 z4) 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 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define polys0120 (rotate-facet polys1200))
 (define polys0012 (rotate-facet polys0120))
@@ -323,10 +343,10 @@ above.
 
 (define (polys1022 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (solve-t za z2 z3) za)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (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))))
 
 (define polys2102 (rotate-facet polys1022))
 (define polys2210 (rotate-facet polys2102))
@@ -336,36 +356,36 @@ above.
 ;; single hexagon
 
 (define (polys0112 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
+  (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 (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) 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))))
 
 (define polys2011 (rotate-facet polys0112))
 (define polys1201 (rotate-facet polys2011))
 (define polys1120 (rotate-facet polys1201))
 
 (define (polys2110 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t zb z1 z2) 0.0 zb)
+  (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 (solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (solve-t za z1 z4) za)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (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))))
 
 (define polys0211 (rotate-facet polys2110))
 (define polys1021 (rotate-facet polys0211))
 (define polys1102 (rotate-facet polys1021))
 
 (define (polys0121 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
+  (list (list (vector (unsafe-solve-t za z1 z2) 0.0 za)
               (vector 1.0 0.0 z2)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector (solve-t zb z4 z3) 1.0 zb)
+              (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 (solve-t za z1 z4) za))))
+              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define polys1012 (rotate-facet polys0121))
 (define polys2101 (rotate-facet polys1012))
@@ -376,19 +396,19 @@ above.
 
 (define (polys10100 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (solve-t za z1 z4) za))
+              (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 (solve-t za z4 z3) 1.0 za)
-              (vector 1.0 (solve-t za z2 z3) za))))
+              (vector (unsafe-solve-t za z4 z3) 1.0 za)
+              (vector 1.0 (unsafe-solve-t za z2 z3) za))))
 
 (define (polys10101 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (solve-t za z2 z3) za)
+              (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 (solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (solve-t za z1 z4) za))))
+              (vector (unsafe-solve-t za z4 z3) 1.0 za)
+              (vector 0.0 (unsafe-solve-t za z1 z4) za))))
 
 (define (polys1010 za zb z1 z2 z3 z4)
   (define z5 (unsafe-flavg4 z1 z2 z3 z4))
@@ -400,19 +420,19 @@ above.
 
 (define (polys1212-2 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))
+              (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 (solve-t zb z4 z3) 1.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb))))
+              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
+              (vector 1.0 (unsafe-solve-t zb z2 z3) zb))))
 
 (define (polys1212-1 za zb z1 z2 z3 z4)
   (list (list (vector 0.0 0.0 z1)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb)
+              (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 (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (vector (unsafe-solve-t zb z4 z3) 1.0 zb)
+              (vector 0.0 (unsafe-solve-t zb z1 z4) zb))))
 
 (define (polys1212 za zb z1 z2 z3 z4)
   (define z5 (unsafe-flavg4 z1 z2 z3 z4))
@@ -426,22 +446,22 @@ above.
 ;; 7-sided saddle
 
 (define (polys0212-1 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb)
+  (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 (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) 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))))
 
 (define (polys0212-2 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) za))
-        (list (vector 1.0 (solve-t zb z2 z3) zb)
+  (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 (solve-t zb z4 z3) 1.0 zb))))
+              (vector (unsafe-solve-t zb z4 z3) 1.0 zb))))
 
 (define (polys0212 za zb z1 z2 z3 z4)
   (define z5 (unsafe-flavg4 z1 z2 z3 z4))
@@ -454,22 +474,22 @@ above.
 (define polys2120 (rotate-facet polys1202))
 
 (define (polys2010-1 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 1.0 (solve-t za z2 z3) za)
+  (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 (solve-t za z4 z3) 1.0 za)
-              (vector 0.0 (solve-t za z1 z4) za)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+              (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))))
 
 (define (polys2010-0 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector (solve-t za z1 z2) 0.0 za)
-              (vector 0.0 (solve-t za z1 z4) za)
-              (vector 0.0 (solve-t zb z1 z4) zb))
-        (list (vector 1.0 (solve-t za z2 z3) za)
+  (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 (solve-t za z4 z3) 1.0 za))))
+              (vector (unsafe-solve-t za z4 z3) 1.0 za))))
 
 (define (polys2010 za zb z1 z2 z3 z4)
   (define z5 (unsafe-flavg4  z1 z2 z3 z4))
@@ -485,34 +505,34 @@ above.
 ;; 8-sided saddle
 
 (define (polys0202-0 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 1.0 (solve-t za z2 z3) za))
-        (list (vector 0.0 (solve-t za z1 z4) za)
-              (vector (solve-t za z4 z3) 1.0 za)
-              (vector (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb))))
+  (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))))
 
 (define (polys0202-1 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 1.0 (solve-t za z2 z3) za)
-              (vector (solve-t za z4 z3) 1.0 za)
-              (vector (solve-t zb z4 z3) 1.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) za))))
+  (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))))
 
 (define (polys0202-2 za zb z1 z2 z3 z4)
-  (list (list (vector (solve-t za z1 z2) 0.0 za)
-              (vector (solve-t zb z1 z2) 0.0 zb)
-              (vector 0.0 (solve-t zb z1 z4) zb)
-              (vector 0.0 (solve-t za z1 z4) za))
-        (list (vector 1.0 (solve-t zb z2 z3) zb)
-              (vector 1.0 (solve-t za z2 z3) za)
-              (vector (solve-t za z4 z3) 1.0 za)
-              (vector (solve-t zb z4 z3) 1.0 zb))))
+  (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))))
 
 (define (polys0202 za zb z1 z2 z3 z4)
   (define z5 (unsafe-flavg4 z1 z2 z3 z4))
@@ -563,26 +583,53 @@ 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 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))
+                            (unsafe-fx* t3 3))
+                t4))
+  (define f (vector-ref polys-dispatch-table facet-num))
+  (f za zb z1 z2 z3 z4))
+
 (defproc (heights->polys [xa real?] [xb real?] [ya real?] [yb real?]
                          [za real?] [zb real?]
                          [z1 real?] [z2 real?] [z3 real?] [z4 real?]
                          ) (listof (vector/c real? real? real?))
-  (check-all-real! 'heights->polys xa xb ya yb za zb z1 z2 z3 z4)
-  (let-exact->inexact
-   (xa xb ya yb 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 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))
-                             (unsafe-fx* t3 3))
-                 t4))
-   (define f (vector-ref polys-dispatch-table facet-num))
-   (for/list ([poly  (in-list (f za zb z1 z2 z3 z4))])
-     (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))]))))
+  (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 (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))))])))]))
diff --git a/collects/plot/common/marching-utils.rkt b/collects/plot/common/marching-utils.rkt
index 09a0af0b44..6aeace039a 100644
--- a/collects/plot/common/marching-utils.rkt
+++ b/collects/plot/common/marching-utils.rkt
@@ -1,6 +1,7 @@
 #lang racket/base
 
-(require (for-syntax racket/base) racket/unsafe/ops)
+(require (for-syntax racket/base) racket/unsafe/ops
+         "utils.rkt")
 
 (provide (all-defined-out))
 
@@ -9,24 +10,24 @@
 ;;           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)
+(define-syntax-rule (unsafe-solve-t z za zb)
   (unsafe-fl/ (unsafe-fl- z za) (unsafe-fl- zb za)))
 
-(define-syntax-rule (unsolve-t za zb t)
+(define-syntax-rule (unsafe-unsolve-t za zb t)
   (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))))
+
 (define-syntax-rule (unsafe-flavg4 z1 z2 z3 z4)
   (unsafe-fl* 0.25 (unsafe-fl+ (unsafe-fl+ (unsafe-fl+ z1 z2) z3) z4)))
 
-(define-syntax (check-all-real! stx)
-  (syntax-case stx ()
-    [(_ name id ...)
-     (with-syntax ([(i ...)  (build-list (length (syntax->list #'(id ...))) values)])
-       (define ids (syntax->list #'(id ...)))
-       #`(begin (unless (real? id)
-                  (raise-type-error name "real number" i #,@ids))
-                ...))]))
+(define-syntax-rule (all pred? obj ...) (and (pred? obj) ...))
 
-(define-syntax-rule (let-exact->inexact (id ...) body0 body ...)
-  (let ([id  (exact->inexact id)] ...)
-    body0 body ...))
+(define-syntax (find-failure-index stx)
+  (syntax-case stx ()
+    [(_ pred? id ...)
+     (with-syntax ([(i ...)   (build-list (length (syntax->list #'(id ...))) values)])
+       (syntax/loc stx
+         (or (and (not (pred? id)) i) ...)))]))
diff --git a/collects/plot/common/math.rkt b/collects/plot/common/math.rkt
index ec8b57e5c8..4d41b2a663 100644
--- a/collects/plot/common/math.rkt
+++ b/collects/plot/common/math.rkt
@@ -54,12 +54,79 @@
     [xs  (cond [(not (andmap flonum? xs))  (raise-type-error 'fldistance "flonums" xs)]
                [else  (unsafe-flsqrt (flsum (λ (x) (unsafe-fl* x x)) xs))])]))
 
+(define (flonum->bit-field x)
+  (integer-bytes->integer (real->floating-point-bytes x 8) #f))
+
+(define (bit-field->flonum i)
+  (floating-point-bytes->real (integer->integer-bytes i 8 #f)))
+
+(defproc (flonum->ordinal [x flonum?]) integer?
+  (cond [(x . < . 0)  (- (flonum->bit-field (- x)))]
+        [else         (flonum->bit-field (abs x))]))
+
+(defproc (ordinal->flonum [i (integer-in #x-7FFFFFFFFFFFFFFF #x7FFFFFFFFFFFFFFF)]) flonum?
+  (cond [(i . < . 0)  (- (bit-field->flonum (- i)))]
+        [else         (bit-field->flonum i)]))
+
+(define +inf-ordinal (flonum->ordinal +inf.0))
+(define -inf-ordinal (flonum->ordinal -inf.0))
+
+(defproc (flstep [x flonum?] [n exact-integer?]) flonum?
+  (cond [(eqv? x +nan.0)  +nan.0]
+        [(and (eqv? x +inf.0) (n . >= . 0))  +inf.0]
+        [(and (eqv? x -inf.0) (n . <= . 0))  -inf.0]
+        [else
+         (define i (+ n (flonum->ordinal x)))
+         (cond [(i . < . -inf-ordinal)  -inf.0]
+               [(i . > . +inf-ordinal)  +inf.0]
+               [else  (ordinal->flonum i)])]))
+
+(defproc (flnext [x flonum?]) flonum? #:document-body
+  (flstep x 1))
+
+(defproc (flprev [x flonum?]) flonum? #:document-body
+  (flstep x -1))
+
+(defthing +min.0 flonum? #:document-value (flnext 0.0))
+(defthing -min.0 flonum? #:document-value (flprev 0.0))
+(defthing +max.0 flonum? #:document-value (flprev +inf.0))
+(defthing -max.0 flonum? #:document-value (flnext -inf.0))
+
 ;; ===================================================================================================
 ;; Reals
 
 (defproc (regular-real? [x any/c]) boolean?
   (and (real? x) (not (special-real? x))))
 
+(defproc (maybe-inexact->exact [x (or/c regular-real? #f)]) (or/c regular-real? #f)
+  (cond [x  (unless (regular-real? x)
+              (raise-type-error 'maybe-inexact->exact "regular real or #f" x))
+            (inexact->exact x)]
+        [else  #f]))
+
+(defproc (flonum-ok-for-range? [x-min regular-real?] [x-max regular-real?]
+                               [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 (regular-real? inexact-x-min) (return #f))
+      
+      (define inexact-x-max (exact->inexact x-max))
+      (unless (regular-real? inexact-x-max) (return #f))
+      
+      (define inexact-x-max-prev (flprev inexact-x-max))
+      (unless (regular-real? inexact-x-max-prev) (return #f))
+      
+      (define inexact-x-min-next (flnext inexact-x-min))
+      (unless (regular-real? 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)))))
+
 (define equal?*
   (case-lambda
     [()   #t]
@@ -154,25 +221,23 @@
   (cond [(not (and (exact-integer? b) (b . >= . 2)))  (raise-type-error 'floor-log/base
                                                                         "exact integer >= 2" 0 b x)]
         [(not (and (real? x) (x . > . 0)))  (raise-type-error 'floor-log/base "real > 0" 1 b x)]
-        [else  (define y (inexact->exact (floor (/ (log x) (log b)))))
-               (cond [(exact? x)
-                      (let loop ([y y] [x  (/ x (expt b y))])
+        [else  (cond [(exact? x)
+                      (let loop ([y 0] [x x])
                         (cond [(x . >= . b)  (loop (add1 y) (/ x b))]
                               [(x . < . 1)   (loop (sub1 y) (* x b))]
                               [else  y]))]
-                     [else  y])]))
+                     [else  (inexact->exact (floor (/ (log x) (log b))))])]))
 
 (defproc (ceiling-log/base [b (and/c exact-integer? (>=/c 2))] [x (>/c 0)]) exact-integer?
   (cond [(not (and (exact-integer? b) (b . >= . 2)))  (raise-type-error 'floor-log/base
                                                                         "exact integer >= 2" 0 b x)]
         [(not (and (real? x) (x . > . 0)))  (raise-type-error 'floor-log/base "real > 0" 1 b x)]
-        [else  (define y (inexact->exact (ceiling (/ (log x) (log b)))))
-               (cond [(exact? x)
-                      (let loop ([y y] [x  (/ x (expt b y))])
+        [else  (cond [(exact? x)
+                      (let loop ([y 0] [x x])
                         (cond [(x . > . 1)         (loop (add1 y) (/ x b))]
                               [(x . <= . (/ 1 b))  (loop (sub1 y) (* x b))]
                               [else  y]))]
-                     [else  y])]))
+                     [else  (inexact->exact (ceiling (/ (log x) (log b))))])]))
 
 (defproc (polar->cartesian [θ real?] [r real?]) (vector/c real? real?)
   (cond [(not (real? θ))  (raise-type-error 'polar->cartesian "real number" 0 θ r)]
@@ -367,14 +432,12 @@
 (defproc (vregular? [v (vectorof real?)]) boolean?
   (match v
     [(vector (? real? x) (? real? y))
-     (cond [(flonum? x)  (unsafe-flregular? x)]
-           [(flonum? y)  (unsafe-flregular? y)]
-           [else  #t])]
+     (not (or (and (flonum? x) (unsafe-flspecial? x))
+              (and (flonum? y) (unsafe-flspecial? y))))]
     [(vector (? real? x) (? real? y) (? real? z))
-     (cond [(flonum? x)  (unsafe-flregular? x)]
-           [(flonum? y)  (unsafe-flregular? y)]
-           [(flonum? z)  (unsafe-flregular? z)]
-           [else  #t])]
+     (not (or (and (flonum? x) (unsafe-flspecial? x))
+              (and (flonum? y) (unsafe-flspecial? y))
+              (and (flonum? z) (unsafe-flspecial? z))))]
     [_  (cond [(vector-andmap real? v)  (let/ec break
                                           (for ([x  (in-vector v)])
                                             (when (and (flonum? x) (unsafe-flspecial? x))
@@ -435,8 +498,6 @@
   (cond [(and (not (null? res)) (null? (car res)))  (cdr res)]
         [else  res]))
 
-(define default-normal (vector 0 -1 0))
-
 (define (remove-degenerate-edges vs)
   (cond
     [(empty? vs)  empty]
@@ -450,6 +511,8 @@
        (cond [(v= last (first vs))  (rest vs)]
              [else  vs]))]))
 
+(define default-normal (vector 0.0 -1.0 0.0))
+
 (define (vnormal vs)
   (let ([vs  (remove-degenerate-edges vs)])
     (cond
@@ -515,7 +578,8 @@
 
 (defproc (ivl-inexact->exact [i ivl?]) ivl?
   (match-define (ivl a b) i)
-  (ivl (inexact->exact a) (inexact->exact b)))
+  (ivl (and a (if (nan? a) a (inexact->exact a)))
+       (and b (if (nan? b) b (inexact->exact b)))))
 
 (defproc (ivl-contains? [i ivl?] [x real?]) boolean?
   (match-define (ivl a b) i)
diff --git a/collects/plot/common/plot-device.rkt b/collects/plot/common/plot-device.rkt
index c79747a6cf..a352868427 100644
--- a/collects/plot/common/plot-device.rkt
+++ b/collects/plot/common/plot-device.rkt
@@ -7,6 +7,7 @@
 ;; It is up to callers to transform view or plot coordinates into dc coordinates.
 
 (require racket/draw racket/class racket/match racket/math racket/bool racket/list racket/contract
+         racket/vector
          "contract.rkt"
          "draw.rkt"
          "math.rkt"
@@ -308,9 +309,11 @@
         (draw-text/anchor dc str x y anchor #t 0 angle)))
     
     (define/public (get-text-corners str v [anchor 'top-left] [angle 0])
-      (when (vregular? v)
-        (match-define (vector x y) v)
-        (get-text-corners/anchor dc str x y anchor #t 0 angle)))
+      (cond [(vregular? v)
+             (match-define (vector x y) v)
+             (map (λ (v) (vector-map inexact->exact v))
+                  (get-text-corners/anchor dc str x y anchor #t 0 angle))]
+            [else  empty]))
     
     (define/public (draw-arrow v1 v2)
       (when (and (vregular? v1) (vregular? v2))
@@ -370,8 +373,8 @@
     
     (define/public (get-tick-endpoints v r angle)
       (match-define (vector x y) v)
-      (define dx (* (cos angle) r))
-      (define dy (* (sin angle) r))
+      (define dx (* (inexact->exact (cos angle)) r))
+      (define dy (* (inexact->exact (sin angle)) r))
       (list (vector (- x dx) (- y dy)) (vector (+ x dx) (+ y dy))))
     
     (define/public (make-draw-tick r angle)
diff --git a/collects/plot/common/plot-element.rkt b/collects/plot/common/plot-element.rkt
index 1460f5cee6..4ad489d797 100644
--- a/collects/plot/common/plot-element.rkt
+++ b/collects/plot/common/plot-element.rkt
@@ -84,30 +84,34 @@
 ;; More precisely, starting with the given plot bounds, it attempts to compute a fixpoint of
 ;; (apply-bounds* elems), overridden at every iteration by the plot bounds (if given). Because a
 ;; fixpoint doesn't always exist, or only exists in the limit, it stops after max-iters.
-(define (bounds-fixpoint elems plot-bounds-rect [max-iters 4])
+(define (bounds-fixpoint elems given-bounds-rect [max-iters 4])
+  (define res
   (let/ec break
-    ;; Shortcut eval: if the plot bounds are all known, the code below just returns them anyway
-    (when (rect-known? plot-bounds-rect) (break plot-bounds-rect))
-    ;; Objective: find the fixpoint of F starting at plot-bounds-rect
-    (define (F bounds-rect) (rect-meet plot-bounds-rect (apply-bounds* elems bounds-rect)))
-    ;; Iterate joint bounds to (hopefully) a fixpoint
-    (define-values (bounds-rect area delta-area)
-      (for/fold ([bounds-rect  plot-bounds-rect]
-                 [area  (rect-area plot-bounds-rect)] [delta-area  #f]
-                 ) ([n  (in-range max-iters)])
-        ;(printf "bounds-rect = ~v~n" bounds-rect)
-        ;; Get new bounds from the elements' bounds functions
-        (define new-bounds-rect (F bounds-rect))
-        (define new-area (rect-area new-bounds-rect))
-        (define new-delta-area (and area new-area (- new-area area)))
-        (cond
-          ;; Shortcut eval: if the bounds haven't changed, we have a fixpoint
-          [(equal? bounds-rect new-bounds-rect)  (break bounds-rect)]
-          ;; If the area grew more this iteration than last, it may not converge, so stop now
-          [(and delta-area new-delta-area (new-delta-area . > . delta-area))  (break bounds-rect)]
-          ;; All good - one more iteration
-          [else  (values new-bounds-rect new-area new-delta-area)])))
-    bounds-rect))
+    ;; Shortcut eval: if the plot bounds are all given, the code below just returns them anyway
+    (when (rect-known? given-bounds-rect) (break given-bounds-rect))
+    (let ([given-bounds-rect  (rect-inexact->exact given-bounds-rect)])
+      ;; Objective: find the fixpoint of F starting at given-bounds-rect
+      (define (F bounds-rect) (rect-meet given-bounds-rect (apply-bounds* elems bounds-rect)))
+      ;; Iterate joint bounds to (hopefully) a fixpoint
+      (define-values (bounds-rect area delta-area)
+        (for/fold ([bounds-rect  given-bounds-rect]
+                   [area  (rect-area given-bounds-rect)] [delta-area  #f]
+                   ) ([n  (in-range max-iters)])
+          ;(printf "bounds-rect = ~v~n" bounds-rect)
+          ;; Get new bounds from the elements' bounds functions
+          (define new-bounds-rect (F bounds-rect))
+          (define new-area (rect-area new-bounds-rect))
+          (define new-delta-area (and area new-area (- new-area area)))
+          (cond
+            ;; Shortcut eval: if the bounds haven't changed, we have a fixpoint
+            [(equal? bounds-rect new-bounds-rect)  (break bounds-rect)]
+            ;; If the area grew more this iteration than last, it may not converge, so stop now
+            [(and delta-area new-delta-area (new-delta-area . > . delta-area))  (break bounds-rect)]
+            ;; All good - one more iteration
+            [else  (values new-bounds-rect new-area new-delta-area)])))
+      bounds-rect)))
+  ;(printf "fixpoint bounds-rect = ~v~n" res)
+  res)
 
 ;; Applies the bounds functions of multiple plot elements, in parallel, and returns the smallest
 ;; bounds containing all the new bounds. This function is monotone and increasing regardless of
@@ -120,7 +124,12 @@
 ;; what bounds will you try to use?
 (define (apply-bounds elem bounds-rect)
   (match-define (plot-element elem-bounds-rect elem-bounds-fun _) elem)
-  (let ([elem-bounds-rect  (cond [elem-bounds-rect  (rect-meet bounds-rect elem-bounds-rect)]
-                                 [else  bounds-rect])])
-    (cond [elem-bounds-fun  (elem-bounds-fun elem-bounds-rect)]
-          [else  elem-bounds-rect])))
+  ;(printf "elem-bounds-rect = ~v~n" elem-bounds-rect)
+  (let ([elem-bounds-rect  (if elem-bounds-rect
+                               (rect-meet bounds-rect (rect-inexact->exact elem-bounds-rect))
+                               bounds-rect)])
+    (if elem-bounds-fun
+        (let ([new-elem-bounds-rect  (elem-bounds-fun elem-bounds-rect)])
+          ;(printf "new-elem-bounds-rect = ~v~n" new-elem-bounds-rect)
+          (rect-inexact->exact new-elem-bounds-rect))
+        elem-bounds-rect)))
diff --git a/collects/plot/common/sample.rkt b/collects/plot/common/sample.rkt
index d97ea1da93..b896718360 100644
--- a/collects/plot/common/sample.rkt
+++ b/collects/plot/common/sample.rkt
@@ -2,7 +2,7 @@
 
 ;; Functions that sample from functions, and functions that create memoized samplers.
 
-(require racket/match racket/flonum racket/math racket/contract racket/list
+(require racket/match racket/flonum racket/math racket/contract racket/list racket/vector
          "contract-doc.rkt"
          "math.rkt"
          "axis-transform.rkt")
@@ -75,7 +75,6 @@
 (struct 2d-sample (xs ys zss z-min z-max) #:transparent)
 (struct 3d-sample (xs ys zs dsss d-min d-max) #:transparent)
 
-(defcontract sample/c (list/c (listof real?) (listof real?)))
 (defcontract sampler/c (real? real? exact-nonnegative-integer? . -> . sample?))
 
 (defcontract 2d-sampler/c (real? real? exact-nonnegative-integer?
@@ -101,7 +100,9 @@
                    (define-values (y-min y-max)
                      (cond [(empty? rys)  (values #f #f)]
                            [else  (values (apply min* rys) (apply max* rys))]))
-                   (sample xs ys y-min y-max))))))
+                   (sample xs ys
+                           (maybe-inexact->exact y-min)
+                           (maybe-inexact->exact y-max)))))))
 
 (defproc (make-2d-function->sampler [transform-x-thnk (-> axis-transform/c)]
                                     [transform-y-thnk (-> axis-transform/c)]
@@ -124,7 +125,9 @@
                                        (unless (and z-min (z . >= . z-min)) (set! z-min z))
                                        (unless (and z-max (z . <= . z-max)) (set! z-max z)))
                                      z))))
-                   (2d-sample xs ys zss z-min z-max))))))
+                   (2d-sample xs ys zss
+                              (maybe-inexact->exact z-min)
+                              (maybe-inexact->exact z-max)))))))
 
 (defproc (make-3d-function->sampler [transform-x-thnk (-> axis-transform/c)]
                                     [transform-y-thnk (-> axis-transform/c)]
@@ -153,7 +156,9 @@
                                           (unless (and d-min (d . >= . d-min)) (set! d-min d))
                                           (unless (and d-max (d . <= . d-max)) (set! d-max d)))
                                         d)))))
-                   (3d-sample xs ys zs dsss d-min d-max))))))
+                   (3d-sample xs ys zs dsss
+                              (maybe-inexact->exact d-min)
+                              (maybe-inexact->exact d-max)))))))
 
 (define-syntax-rule (for-2d-sample (xa xb ya yb z1 z2 z3 z4) sample expr ...)
   (let ()
@@ -199,3 +204,46 @@
           (values xb d2 d3 d6 d7))
         (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 regular-real?] [x-max regular-real?]
+                            [y-min regular-real?] [y-max regular-real?]) 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 regular-real?] [x-max regular-real?]
+                            [y-min regular-real?] [y-max regular-real?]
+                            [z-min regular-real?] [z-max regular-real?]) 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 regular-real?] [x-max regular-real?]
+                            [y-min regular-real?] [y-max regular-real?]
+                            [z-min regular-real?] [z-max regular-real?]
+                            [d-min regular-real?] [d-max regular-real?]) 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/collects/plot/common/utils.rkt b/collects/plot/common/utils.rkt
index de21c39d57..92335ea682 100644
--- a/collects/plot/common/utils.rkt
+++ b/collects/plot/common/utils.rkt
@@ -81,3 +81,7 @@
         (values (cons (reverse lst) res)
                 rest-xs)))
     (reverse res)))
+
+(define-syntax-rule (let-map (id ...) fun body0 body ...)
+  (let ([id  (fun id)] ...)
+    body0 body ...))
diff --git a/collects/plot/contracted/math.rkt b/collects/plot/contracted/math.rkt
index 40e104d18f..c185bc116a 100644
--- a/collects/plot/contracted/math.rkt
+++ b/collects/plot/contracted/math.rkt
@@ -7,8 +7,11 @@
          ;; Flonums
          nan? infinite? special-real?
          flblend flatan2 flsum flmodulo fldistance
+         (activate-contract-out flonum->ordinal ordinal->flonum flstep flnext flprev
+                                flonum-ok-for-range?)
+         -max.0 -min.0 +min.0 +max.0
          ;; Reals
-         regular-real?
+         regular-real? maybe-inexact->exact
          min* max* degrees->radians radians->degrees blend atan2 sum real-modulo distance
          floor-log/base ceiling-log/base
          polar->cartesian 3d-polar->3d-cartesian
diff --git a/collects/plot/contracted/sample.rkt b/collects/plot/contracted/sample.rkt
index 6a28a7fea4..93668689b6 100644
--- a/collects/plot/contracted/sample.rkt
+++ b/collects/plot/contracted/sample.rkt
@@ -23,7 +23,14 @@
                                 sampler/c 2d-sampler/c 3d-sampler/c
                                 make-function->sampler
                                 make-2d-function->sampler
-                                make-3d-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))])))
          map*
diff --git a/collects/plot/plot2d/contour.rkt b/collects/plot/plot2d/contour.rkt
index 9c8a584ce0..a4a7180261 100644
--- a/collects/plot/plot2d/contour.rkt
+++ b/collects/plot/plot2d/contour.rkt
@@ -21,10 +21,13 @@
   (when (<= z-min z z-max)
     (send area put-alpha alpha)
     (send area put-pen color width style)
-    (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)))))
+    (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))))))
   
   (cond [label  (line-legend-entry label color width style)]
         [else   empty]))
@@ -57,10 +60,13 @@
     (match-define (2d-sample xs ys zss z-min z-max) sample)
     (match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #f))
     
-    (let ([colors  (maybe-apply colors zs)]
-          [widths  (maybe-apply widths zs)]
-          [styles  (maybe-apply styles zs)]
-          [alphas  (maybe-apply alphas zs)])
+    (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)])
       (for ([z      (in-list zs)]
             [color  (in-cycle colors)]
             [width  (in-cycle widths)]
@@ -116,9 +122,12 @@
         (values (ivl za zb) (format "[~a,~a]" la lb))))
     
     (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)])
+    (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)])
       (for ([za     (in-list zs)]
             [zb     (in-list (rest zs))]
             [color  (in-cycle colors)]
diff --git a/collects/plot/plot2d/plot.rkt b/collects/plot/plot2d/plot.rkt
index 510b1f2e65..c7bf219e8b 100644
--- a/collects/plot/plot2d/plot.rkt
+++ b/collects/plot/plot2d/plot.rkt
@@ -35,14 +35,12 @@
 (define (get-bounds-rect renderer-list x-min x-max y-min y-max)
   (define given-bounds-rect (vector (ivl x-min x-max) (ivl y-min y-max)))
   (define plot-bounds-rect (bounds-fixpoint renderer-list given-bounds-rect))
-  
   (when (or (not (rect-regular? plot-bounds-rect))
             (rect-zero-area? plot-bounds-rect))
     (match-define (vector (ivl x-min x-max) (ivl y-min y-max)) plot-bounds-rect)
     (error 'plot "could not determine sensible plot bounds; got x ∈ [~a,~a], y ∈ [~a,~a]"
            x-min x-max y-min y-max))
-  
-  (rect-inexact->exact plot-bounds-rect))
+  plot-bounds-rect)
 
 (define (get-ticks renderer-list bounds-rect)
   (define-values (all-x-ticks all-x-far-ticks all-y-ticks all-y-far-ticks)
@@ -51,7 +49,6 @@
       (define ticks-fun (plot-element-ticks-fun r))
       (cond [ticks-fun  (ticks-fun bounds-rect)]
             [else       (values empty empty empty empty)])))
-  
   (values (remove-duplicates (append* all-x-ticks))
           (remove-duplicates (append* all-x-far-ticks))
           (remove-duplicates (append* all-y-ticks))
diff --git a/collects/plot/plot3d/contour.rkt b/collects/plot/plot3d/contour.rkt
index 3ae3782fcd..60b4ce8b0d 100644
--- a/collects/plot/plot3d/contour.rkt
+++ b/collects/plot/plot3d/contour.rkt
@@ -19,13 +19,16 @@
   (when (<= z-min z z-max)
     (send area put-alpha alpha)
     (send area put-pen color width style)
-    (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)
-       (define center (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb))
-                              (* 1/2 (+ (min z1 z2 z3 z4) (max z1 z2 z3 z4)))))
-       (send area put-line v1 v2 center))))
+    (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)
+         (define center (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb))
+                                (* 1/2 (+ (min z1 z2 z3 z4) (max z1 z2 z3 z4)))))
+         (send area put-line v1 v2 center)))))
   
   (cond [label  (line-legend-entry label color width style)]
         [else   empty]))
@@ -62,10 +65,13 @@
                     y-min y-max (animated-samples samples)))
   (match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #f))
   
-  (let ([colors  (maybe-apply colors zs)]
-        [widths  (maybe-apply widths zs)]
-        [styles  (maybe-apply styles zs)]
-        [alphas  (maybe-apply alphas zs)])
+  (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)])
     (for ([z  (in-list zs)]
           [color  (in-cycle colors)]
           [width  (in-cycle widths)]
@@ -115,9 +121,8 @@
          area)
   (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
     (send area get-bounds-rect))
-  (define sample (f x-min x-max (animated-samples samples) y-min y-max (animated-samples samples)))
-  (match-define (2d-sample xs ys zss fz-min fz-max) sample)
-  
+  (define sample (f x-min x-max (animated-samples samples)
+                    y-min y-max (animated-samples samples)))
   (match-define (list (tick zs _ labels) ...) (contour-ticks z-min z-max levels #t))
   
   (define-values (z-ivls ivl-labels)
@@ -127,12 +132,15 @@
                                     [lb  (in-list (rest labels))])
       (values (ivl za zb) (format "[~a,~a]" la lb))))
   
-  (let ([colors  (maybe-apply colors z-ivls)]
-        [styles  (maybe-apply styles z-ivls)]
-        [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)])
+  (let* ([colors  (maybe-apply colors z-ivls)]
+         [styles  (maybe-apply styles z-ivls)]
+         [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)])
     (for ([za  (in-list zs)]
           [zb  (in-list (rest zs))]
           [color  (in-cycle colors)]
diff --git a/collects/plot/plot3d/isosurface.rkt b/collects/plot/plot3d/isosurface.rkt
index 3cca673655..19df425861 100644
--- a/collects/plot/plot3d/isosurface.rkt
+++ b/collects/plot/plot3d/isosurface.rkt
@@ -22,12 +22,15 @@
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (for-3d-sample
-   (xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8) sample
-   (define polys (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
-   (when (not (empty? polys))
-     (send area put-polygons polys
-           (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
+  (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
+     (define polys (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
+     (when (not (empty? polys))
+       (send area put-polygons polys
+             (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb)))))))
   
   (cond [label  (rectangle-legend-entry
                  label color style line-color line-width line-style)]
@@ -76,13 +79,16 @@
      (match-define (list (tick ds _ labels) ...) (isosurface-ticks d-min d-max levels))
      #;(define ds (linear-seq d-min d-max levels #:start? (and rd-min #t) #:end? (and rd-max #t)))
      
-     (let ([colors  (maybe-apply colors ds)]
-           [styles  (maybe-apply styles ds)]
-           [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)])
-       (for ([d           (in-list ds)]
+     (let* ([colors  (maybe-apply colors ds)]
+            [styles  (maybe-apply styles ds)]
+            [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)])
+       (for ([d      (in-list ds)]
              [color  (in-cycle colors)]
              [style  (in-cycle styles)]
              [alpha  (in-cycle alphas)]
@@ -137,8 +143,8 @@
                     z-min z-max (animated-samples samples)))
   (match-define (3d-sample xs ys zs dsss d-min d-max) sample)
   
-  (define (draw-cube xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8)
-    (define polys (heights->cube-polys xa xb ya yb za zb 0.0 d1 d2 d3 d4 d5 d6 d7 d8))
+  (define (draw-cube xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8)
+    (define polys (heights->cube-polys xa xb ya yb za zb d d1 d2 d3 d4 d5 d6 d7 d8))
     (when (not (empty? polys))
       (send area put-polygons polys
             (vector (* 1/2 (+ xa xb)) (* 1/2 (+ ya yb)) (* 1/2 (+ za zb))))))
@@ -146,23 +152,26 @@
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (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 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 d1 d2 d3 d4 d5 d6 d7 d8))]
-         [else
-          (draw-cube xa xb ya yb za zb d1 d2 d3 d4 d5 d6 d7 d8)]))
+  (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)])))
   
   (cond [label  (rectangle-legend-entry
                  label color style line-color line-width line-style)]
diff --git a/collects/plot/plot3d/matrix.rkt b/collects/plot/plot3d/matrix.rkt
index 5b28403880..f0f1076079 100644
--- a/collects/plot/plot3d/matrix.rkt
+++ b/collects/plot/plot3d/matrix.rkt
@@ -7,7 +7,7 @@
 
 (provide m3-apply m3-transpose m3* m3-rotate-z m3-rotate-x m3-scale)
 
-(define-syntax-rule (dot x1 y1 z1 x2 y2 z2) (+ (* x1 x2) (* y1 y2) (* z1 z2)))
+(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 (vector v11 v12 v13) (vector v21 v22 v23) (vector v31 v32 v33)) m)
@@ -28,20 +28,22 @@
   (vector (m3-apply m v1) (m3-apply m v2) (m3-apply m v3)))
 
 (define (m3-rotate-z theta)
-  (define cos-theta (cos theta))
-  (define sin-theta (sin theta))
-  (vector (vector cos-theta (- sin-theta) 0.0)
-          (vector sin-theta cos-theta 0.0)
-          (vector 0.0 0.0 1.0)))
+  (let ([theta  (exact->inexact theta)])
+    (define cos-theta (flcos theta))
+    (define sin-theta (flsin theta))
+    (vector (vector cos-theta (- sin-theta) 0.0)
+            (vector sin-theta cos-theta 0.0)
+            (vector 0.0 0.0 1.0))))
 
 (define (m3-rotate-x rho)
-  (define cos-rho (cos rho))
-  (define sin-rho (sin rho))
-  (vector (vector 1.0 0.0 0.0)
-          (vector 0.0 cos-rho (- sin-rho))
-          (vector 0.0 sin-rho cos-rho)))
+  (let ([rho  (exact->inexact rho)])
+    (define cos-rho (flcos rho))
+    (define sin-rho (flsin rho))
+    (vector (vector 1.0 0.0 0.0)
+            (vector 0.0 cos-rho (- sin-rho))
+            (vector 0.0 sin-rho cos-rho))))
 
 (define (m3-scale x-scale y-scale z-scale)
-  (vector (vector x-scale 0.0 0.0)
-          (vector 0.0 y-scale 0.0)
-          (vector 0.0 0.0 z-scale)))
+  (vector (vector (exact->inexact x-scale) 0.0 0.0)
+          (vector 0.0 (exact->inexact y-scale) 0.0)
+          (vector 0.0 0.0 (exact->inexact z-scale))))
diff --git a/collects/plot/plot3d/plot-area.rkt b/collects/plot/plot3d/plot-area.rkt
index f0fc8b90a2..bf835194c1 100644
--- a/collects/plot/plot3d/plot-area.rkt
+++ b/collects/plot/plot3d/plot-area.rkt
@@ -1,6 +1,6 @@
 #lang racket/base
 
-(require racket/class racket/match racket/list racket/math racket/contract racket/vector
+(require racket/class racket/match racket/list racket/math racket/contract racket/vector racket/flonum
          "../common/math.rkt"
          "../common/plot-device.rkt"
          "../common/ticks.rkt"
@@ -9,6 +9,7 @@
          "../common/axis-transform.rkt"
          "../common/parameters.rkt"
          "../common/sample.rkt"
+         "../common/utils.rkt"
          "matrix.rkt"
          "shape.rkt"
          "clip.rkt")
@@ -98,53 +99,60 @@
     (define theta (+ (degrees->radians angle) 0.00001))
     (define rho (degrees->radians altitude))
     
-    (define identity-transforms?
-      (and (equal? (plot-x-transform) id-transform)
-           (equal? (plot-y-transform) id-transform)
-           (equal? (plot-z-transform) id-transform)))
+    ;; 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
+    ;;  3. View coordinates (from normalized coordinates: rotate)
+    ;;  4. Device context coordinates (from view coordinates: project to 2D)
     
     (match-define (invertible-function fx _) (apply-axis-transform (plot-x-transform) x-min x-max))
     (match-define (invertible-function fy _) (apply-axis-transform (plot-y-transform) y-min y-max))
     (match-define (invertible-function fz _) (apply-axis-transform (plot-z-transform) z-min z-max))
     
-    (define axis-transform
-      (cond
-        [identity-transforms?  (λ (v) v)]
-        [else  (λ (v)
-                 (match-define (vector x y z) v)
-                 (vector (fx x) (fy y) (fz z)))]))
+    (define identity-transforms?
+      (and (equal? (plot-x-transform) id-transform)
+           (equal? (plot-y-transform) id-transform)
+           (equal? (plot-z-transform) id-transform)))
     
-    (define (center v)
-      (match-define (vector x y z) v)
-      (vector (- x x-mid) (- y y-mid) (- z z-mid)))
+    (define flonum-ok? (flonum-ok-for-3d? x-min x-max y-min y-max z-min z-max))
     
-    ;; There are four coordinate systems:
-    ;;  1. Plot coordinates (original, user-facing coordinate system)
-    ;;  2. Normalized coordinates (axis-transformed, centered, and scaled to [-1,1] on each axis)
-    ;;  3. View coordinates (normalized coordinates, rotated)
-    ;;  4. DC coordinates (view coordinates, projected to 2D)
+    (define plot->norm
+      (cond [identity-transforms?
+             (cond [flonum-ok?
+                    (let-map
+                     (x-mid y-mid z-mid x-size y-size z-size) exact->inexact
+                     (λ (v)
+                       (match-define (vector x y z) v)
+                       (vector (fl/ (fl- (exact->inexact x) x-mid) x-size)
+                               (fl/ (fl- (exact->inexact y) y-mid) y-size)
+                               (fl/ (fl- (exact->inexact z) z-mid) z-size))))]
+                   [else
+                    (λ (v)
+                      (match-define (vector x y z) v)
+                      (vector (exact->inexact (/ (- x x-mid) x-size))
+                              (exact->inexact (/ (- y y-mid) y-size))
+                              (exact->inexact (/ (- z z-mid) z-size))))])]
+            [else
+             (λ (v)
+               (match-define (vector x y z) v)
+               (vector (exact->inexact (/ (- (fx x) x-mid) x-size))
+                       (exact->inexact (/ (- (fy y) y-mid) y-size))
+                       (exact->inexact (/ (- (fz z) z-mid) z-size))))]))
     
-    ;; View coordinates show up mostly in tick decorations. Normalized coordinates are only used for
-    ;; surface normals. Most user vertexes get transformed from plot coordinates directly to DC
-    ;; coordinates.
-    
-    (define scale-matrix (m3-scale (/ x-size) (/ y-size) (/ z-size)))
     (define rotate-theta-matrix (m3-rotate-z theta))
     (define rotate-rho-matrix (m3-rotate-x rho))
     (define rotation-matrix (m3* rotate-rho-matrix rotate-theta-matrix))
-    (define transform-matrix (m3* rotation-matrix scale-matrix))
     
-    (define (plot->norm v) (m3-apply scale-matrix (center (axis-transform v))))
     (define (norm->view v) (m3-apply rotation-matrix v))
-    (define (plot->view v) (m3-apply transform-matrix (center (axis-transform v))))
-    
-    (define transform-matrix/no-rho (m3* (m3-rotate-z theta) scale-matrix))
-    (define (plot->view/no-rho v) (m3-apply transform-matrix/no-rho (center (axis-transform v))))
+    (define (plot->view v) (norm->view (plot->norm v)))
+    (define (plot->view/no-rho v) (m3-apply rotate-theta-matrix (plot->norm v)))
+    (define (norm->view/no-rho v) (m3-apply rotate-theta-matrix v))
     (define (rotate/rho v) (m3-apply rotate-rho-matrix v))
     
     (define view->dc #f)
-    (define (plot->dc/no-axis-trans v) (view->dc (m3-apply transform-matrix (center v))))
     (define (plot->dc v) (view->dc (plot->view v)))
+    (define (norm->dc v) (view->dc (norm->view v)))
     
     (define-values (view-x-size view-y-size view-z-size)
       (match-let ([(vector view-x-ivl view-y-ivl view-z-ivl)
@@ -164,30 +172,34 @@
       (define area-y-mid (* 1/2 (+ area-y-min area-y-max)))
       (define area-per-view-x (/ (- area-x-max area-x-min) view-x-size))
       (define area-per-view-z (/ (- area-y-max area-y-min) view-z-size))
-      (λ (v)
-        (match-define (vector x _ z) v)
-        (vector (+ area-x-mid (* x area-per-view-x))
-                (- area-y-mid (* z area-per-view-z)))))
+      (let-map
+       (area-x-mid area-y-mid area-per-view-x area-per-view-z) exact->inexact
+       (λ (v)
+         (match-define (vector x _ z) v)
+         (vector (fl+ area-x-mid (fl* x area-per-view-x))
+                 (fl- area-y-mid (fl* z area-per-view-z))))))
     
     ;; Initial view->dc
     (define init-top-margin (if (and (plot-decorations?) (plot-title)) (* 3/2 char-height) 0))
     (set! view->dc (make-view->dc 0 0 init-top-margin 0))
     
-    (define (plot-dir->dc-angle v)
-      (match-define (vector dx dy)
-        (v- (plot->dc/no-axis-trans (v+ v (vector x-mid y-mid z-mid)))
-            (plot->dc/no-axis-trans (vector x-mid y-mid z-mid))))
+    (define (x-axis-angle)
+      (match-define (vector dx dy) (v- (norm->dc (vector 0.5 0.0 0.0))
+                                       (norm->dc (vector -0.5 0.0 0.0))))
       (- (atan2 (- dy) dx)))
     
-    (define (x-axis-angle) (plot-dir->dc-angle #(1 0 0)))
-    (define (y-axis-angle) (plot-dir->dc-angle #(0 1 0)))
+    (define (y-axis-angle)
+      (match-define (vector dx dy) (v- (norm->dc (vector 0.0 0.5 0.0))
+                                       (norm->dc (vector 0.0 -0.5 0.0))))
+      (- (atan2 (- dy) dx)))
     
-    (define/public (plot-dir->dc-dir v)
-      (vnormalize (v- (plot->dc/no-axis-trans (v+ v (vector x-mid y-mid z-mid)))
-                      (plot->dc/no-axis-trans (vector x-mid y-mid z-mid)))))
+    (define (x-axis-dir)
+      (vnormalize (v- (norm->dc (vector 0.5 0.0 0.0))
+                      (norm->dc (vector -0.5 0.0 0.0)))))
     
-    (define (x-axis-dir) (plot-dir->dc-dir #(1 0 0)))
-    (define (y-axis-dir) (plot-dir->dc-dir #(0 1 0)))
+    (define (y-axis-dir)
+      (vnormalize (v- (norm->dc (vector 0.0 0.5 0.0))
+                      (norm->dc (vector 0.0 -0.5 0.0)))))
     
     ;; ===============================================================================================
     ;; Tick and label constants
@@ -208,6 +220,16 @@
     (define z-far-axis-x (if x-axis-y-min? x-max x-min))
     (define z-far-axis-y (if y-axis-x-min? y-min y-max))
     
+    (define x-axis-norm-y (if x-axis-y-min? -0.5 0.5))
+    (define y-axis-norm-x (if y-axis-x-min? -0.5 0.5))
+    (define z-axis-norm-x (if x-axis-y-min? -0.5 0.5))
+    (define z-axis-norm-y (if y-axis-x-min? 0.5 -0.5))
+    
+    (define x-far-axis-norm-y (if x-axis-y-min? 0.5 -0.5))
+    (define y-far-axis-norm-x (if y-axis-x-min? 0.5 -0.5))
+    (define z-far-axis-norm-x (if x-axis-y-min? 0.5 -0.5))
+    (define z-far-axis-norm-y (if y-axis-x-min? -0.5 0.5))
+    
     (define near-dist^2 (sqr (* 3 (plot-line-width))))
     (define (vnear? v1 v2)
       ((vmag^2 (v- (plot->dc v1) (plot->dc v2))) . <= . near-dist^2))
@@ -432,13 +454,13 @@
           0))
     
     (define (get-x-label-params)
-      (define v0 (plot->dc/no-axis-trans (vector x-mid x-axis-y z-min)))
+      (define v0 (norm->dc (vector 0.0 x-axis-norm-y -0.5)))
       (define dist (+ max-x-tick-offset (max-x-tick-label-diag) half-char-height))
       (list (plot-x-label) (v+ v0 (v* (y-axis-dir) (if x-axis-y-min? (- dist) dist)))
             'top (- (if x-axis-y-min? 0 pi) (x-axis-angle))))
     
     (define (get-y-label-params)
-      (define v0 (plot->dc/no-axis-trans (vector y-axis-x y-mid z-min)))
+      (define v0 (norm->dc (vector y-axis-norm-x 0.0 -0.5)))
       (define dist (+ max-y-tick-offset (max-y-tick-label-diag) half-char-height))
       (list (plot-y-label) (v+ v0 (v* (x-axis-dir) (if y-axis-x-min? (- dist) dist)))
             'top (- (if y-axis-x-min? pi 0) (y-axis-angle))))
@@ -449,13 +471,13 @@
             'bottom-left 0))
     
     (define (get-x-far-label-params)
-      (define v0 (plot->dc/no-axis-trans (vector x-mid x-far-axis-y z-min)))
+      (define v0 (norm->dc (vector 0.0 x-far-axis-norm-y -0.5)))
       (define dist (+ max-x-far-tick-offset (max-x-far-tick-label-diag) half-char-height))
       (list (plot-x-far-label) (v+ v0 (v* (y-axis-dir) (if x-axis-y-min? dist (- dist))))
             'bottom (- (if x-axis-y-min? 0 pi) (x-axis-angle))))
     
     (define (get-y-far-label-params)
-      (define v0 (plot->dc/no-axis-trans (vector y-far-axis-x y-mid z-min)))
+      (define v0 (norm->dc (vector y-far-axis-norm-x 0.0 -0.5)))
       (define dist (+ max-y-far-tick-offset (max-y-far-tick-label-diag) half-char-height))
       (list (plot-y-far-label) (v+ v0 (v* (x-axis-dir) (if y-axis-x-min? dist (- dist))))
             'bottom (- (if y-axis-x-min? pi 0) (y-axis-angle))))
@@ -519,6 +541,7 @@
       ;(printf "label params = ~v~n" (get-all-label-params))
       ;(printf "tick params = ~v~n" (get-all-tick-params))
       (set! view->dc (make-view->dc left right top bottom))
+      ;(printf "~v~n" (get-all-tick-params))
       (append (append* (map (λ (params) (send/apply pd get-text-corners params))
                             (get-all-label-params)))
               (append* (map (λ (params) (send/apply pd get-tick-endpoints (rest params)))
@@ -545,32 +568,32 @@
         (send pd set-minor-pen)
         (when (plot-x-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector x-min x-axis-y z-min))
-                (plot->dc/no-axis-trans (vector x-max x-axis-y z-min))))
+                (norm->dc (vector -0.5 x-axis-norm-y -0.5))
+                (norm->dc (vector 0.5 x-axis-norm-y -0.5))))
         (when (plot-x-far-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector x-min x-far-axis-y z-min))
-                (plot->dc/no-axis-trans (vector x-max x-far-axis-y z-min))))
+                (norm->dc (vector -0.5 x-far-axis-norm-y -0.5))
+                (norm->dc (vector 0.5 x-far-axis-norm-y -0.5))))
         (when (plot-y-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector y-axis-x y-min z-min))
-                (plot->dc/no-axis-trans (vector y-axis-x y-max z-min))))
+                (norm->dc (vector y-axis-norm-x -0.5 -0.5))
+                (norm->dc (vector y-axis-norm-x 0.5 -0.5))))
         (when (plot-y-far-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector y-far-axis-x y-min z-min))
-                (plot->dc/no-axis-trans (vector y-far-axis-x y-max z-min))))))
+                (norm->dc (vector y-far-axis-norm-x -0.5 -0.5))
+                (norm->dc (vector y-far-axis-norm-x 0.5 -0.5))))))
     
     (define (draw-front-axes)
       (when (plot-decorations?)
         (send pd set-minor-pen)
         (when (plot-z-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector z-axis-x z-axis-y z-min))
-                (plot->dc/no-axis-trans (vector z-axis-x z-axis-y z-max))))
+                (norm->dc (vector z-axis-norm-x z-axis-norm-y -0.5))
+                (norm->dc (vector z-axis-norm-x z-axis-norm-y 0.5))))
         (when (plot-z-far-axis?)
           (send pd draw-line
-                (plot->dc/no-axis-trans (vector z-far-axis-x z-far-axis-y z-min))
-                (plot->dc/no-axis-trans (vector z-far-axis-x z-far-axis-y z-max))))))
+                (norm->dc (vector z-far-axis-norm-x z-far-axis-norm-y -0.5))
+                (norm->dc (vector z-far-axis-norm-x z-far-axis-norm-y 0.5))))))
     
     (define (draw-ticks tick-params)
       (for ([params  (in-list tick-params)])
@@ -590,18 +613,19 @@
     (define (add-shapes! shapes) (set! render-list (append shapes render-list)))
     
     (define (draw-shapes ss)
-      (define s+cs (map (λ (s) (cons s (plot->view/no-rho (shape-center s)))) ss))
+      (define s+cs (map (λ (s) (cons s (norm->view/no-rho (shape-center s)))) ss))
       (for ([s+c  (in-list (depth-sort (reverse s+cs)))])
         (match-define (cons s c) s+c)
         (draw-shape s (rotate/rho c))))
     
-    (define (draw-polygon alpha center vs norm pen-color pen-width pen-style brush-color brush-style)
-      (define-values (diff spec) (get-light-values center (norm->view norm)))
+    (define (draw-polygon alpha center vs normal
+                          pen-color pen-width pen-style brush-color brush-style)
+      (define-values (diff spec) (get-light-values center (norm->view normal)))
       (let ([pen-color  (map (λ (v) (+ (* v diff) spec)) pen-color)]
             [brush-color  (map (λ (v) (+ (* v diff) spec)) brush-color)])
         (send pd set-pen pen-color pen-width pen-style)
         (send pd set-brush brush-color brush-style)
-        (send pd draw-polygon (map (λ (v) (plot->dc v)) vs))))
+        (send pd draw-polygon (map (λ (v) (norm->dc v)) vs))))
     
     (define (draw-shape s center)
       (send pd set-alpha (shape-alpha s))
@@ -609,19 +633,20 @@
         ;; shapes
         [(shapes alpha _ ss)  (draw-shapes ss)]
         ;; polygon
-        [(polygon alpha _ vs norm pen-color pen-width pen-style brush-color brush-style)
-         (draw-polygon alpha center vs norm pen-color pen-width pen-style brush-color brush-style)]
+        [(polygon alpha _ vs normal pen-color pen-width pen-style brush-color brush-style)
+         (draw-polygon alpha center vs normal pen-color pen-width pen-style brush-color brush-style)]
         ;; rectangle
         [(rectangle alpha _ r pen-color pen-width pen-style brush-color brush-style)
          (for ([face  (in-list (rect-visible-faces r theta))])
            (match face
-             [(list norm vs ...)  (draw-polygon alpha center vs norm pen-color pen-width pen-style
-                                                brush-color brush-style)]
+             [(list normal vs ...)  (draw-polygon alpha center vs
+                                                  normal pen-color pen-width pen-style
+                                                  brush-color brush-style)]
              [_  (void)]))]
         ;; line
         [(line alpha _ v1 v2 pen-color pen-width pen-style)
          (send pd set-pen pen-color pen-width pen-style)
-         (send pd draw-line (plot->dc v1) (plot->dc v2))]
+         (send pd draw-line (norm->dc v1) (norm->dc v2))]
         ;; text
         [(text alpha _ anchor angle str font-size font-family color)
          (send pd set-font font-size font-family)
@@ -639,21 +664,14 @@
         ;; arrow glyph
         [(arrow-glyph alpha _ v1 v2 pen-color pen-width pen-style)
          (send pd set-pen pen-color pen-width pen-style)
-         (send pd draw-arrow (plot->dc v1) (plot->dc v2))]
+         (send pd draw-arrow (norm->dc v1) (norm->dc v2))]
         [_  (error 'draw-shapes "shape not implemented: ~e" s)]))
     
-    ;; Use a special view transform for the light so that the light angle is always the same
-    ;; regardless of theta (but rotates rho). This also doesn't do any axis transforms, which could
-    ;; fail; e.g. log transform when the light is at a negative position.
-    (define light-transform-matrix (m3* (m3-rotate-x rho) scale-matrix))
-    (define (plot->light-view v) (m3-apply light-transform-matrix (center v)))
-    
     ;; Light position, in normalized view coordinates: 5 units up, ~3 units back and to the left
     ;; (simulates non-noon daylight conditions)
-    (define light (vector-map exact->inexact
-                              (plot->light-view (vector (- x-min (* 2 x-size))
-                                                        (- y-min (* 2 y-size))
-                                                        (+ z-max (* 5 z-size))))))
+    (define light (m3-apply rotate-rho-matrix (vector (- -0.5 2.0)
+                                                      (- -0.5 2.0)
+                                                      (+ 0.5 5.0))))
     ;; View direction, in normalized view coordinates: many graph widths backward
     (define view-dir (vector 0.0 -50.0 0.0))
     
@@ -663,17 +681,17 @@
     
     (define get-light-values
       (cond
-        [(not (or diffuse-light? specular-light?))  (λ (v norm) (values 1.0 0.0))]
+        [(not (or diffuse-light? specular-light?))  (λ (v normal) (values 1.0 0.0))]
         [else
-         (λ (v norm)
+         (λ (v normal)
            ; common lighting values
            (define light-dir (vnormalize (v- light v)))
            ; diffuse lighting: typical Lambertian surface model
-           (define diff (if diffuse-light? (abs (vdot norm light-dir)) 1.0))
+           (define diff (if diffuse-light? (abs (vdot normal light-dir)) 1.0))
            ; specular highlighting: Blinn-Phong model
            (define spec (cond [specular-light?
                                (define lv (v* (v+ light-dir view-dir) 0.5))
-                               (define cos-angle (/ (abs (vdot norm lv)) (vmag lv)))
+                               (define cos-angle (/ (abs (vdot normal lv)) (vmag lv)))
                                (* 32.0 (expt cos-angle 10.0))]
                               [else  0.0]))
            ; put it all together
@@ -793,12 +811,13 @@
                                    (values v1 v2))])
           (unless (and v1 v2) (return (void)))
           (cond [identity-transforms?
-                 (add-shape! (line alpha c v1 v2
+                 (add-shape! (line alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
                                    pen-color pen-width pen-style))]
                 [else
                  (define vs (subdivide-line plot->dc v1 v2))
                  (for ([v1  (in-list vs)] [v2  (in-list (rest vs))])
-                   (add-shape! (line alpha c v1 v2 pen-color pen-width pen-style)))]))))
+                   (add-shape! (line alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
+                                     pen-color pen-width pen-style)))]))))
     
     (define/public (put-lines vs)
       (for ([vs  (vregular-sublists vs)])
@@ -811,7 +830,7 @@
         (when (or (empty? vs) (not (and (andmap vregular? vs) (vregular? c))))
           (return lst))
         
-        (define norm (vnormal (map plot->norm vs)))
+        (define normal (vnormal (map plot->norm vs)))
         (let* ([vs  (if clipping?
                         (clip-polygon vs clip-x-min clip-x-max
                                       clip-y-min clip-y-max
@@ -819,7 +838,8 @@
                         vs)]
                [vs  (if identity-transforms? vs (subdivide-polygon plot->dc vs))])
           (when (empty? vs) (return lst))
-          (cons (polygon alpha c vs norm pen-color pen-width pen-style brush-color brush-style)
+          (cons (polygon alpha (plot->norm c) (map plot->norm vs)
+                         normal pen-color pen-width pen-style brush-color brush-style)
                 lst))))
     
     (define/public (put-polygon vs [c (vcenter vs)])
@@ -830,33 +850,43 @@
                                            #:when (not (empty? vs)))
                     (add-polygon lst vs (vcenter vs))))
       (when (not (empty? lst))
-        (add-shape! (shapes alpha c lst))))
+        (add-shape! (shapes alpha (plot->norm c) lst))))
     
     (define/public (put-rect r [c (rect-center r)])
       (when (rect-regular? r)
         (let ([r  (rect-meet r bounds-rect)])
-          (add-shape! (rectangle alpha c r pen-color pen-width pen-style brush-color brush-style)))))
+          (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) r)
+          (match-let ([(vector x-min y-min z-min)  (plot->norm (vector x-min y-min z-min))]
+                      [(vector x-max y-max z-max)  (plot->norm (vector x-max y-max z-max))])
+            (add-shape! (rectangle alpha (plot->norm c)
+                                   (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
+                                   pen-color pen-width pen-style brush-color brush-style))))))
     
     (define/public (put-text str v [anchor 'center] [angle 0])
       (when (and (vregular? v) (in-bounds? v))
-        (add-shape! (text alpha v anchor angle str font-size font-family text-foreground))))
+        (add-shape! (text alpha (plot->norm v) anchor angle str
+                          font-size font-family text-foreground))))
     
     (define/public (put-glyphs vs symbol size)
       (for ([v  (in-list vs)])
         (when (and (vregular? v) (in-bounds? v))
           (add-shape!
-           (glyph alpha v symbol size pen-color pen-width pen-style brush-color brush-style)))))
+           (glyph alpha (plot->norm v) symbol size
+                  pen-color pen-width pen-style brush-color brush-style)))))
     
     (define/public (put-arrow v1 v2 [c (v* (v+ v1 v2) 1/2)])
       (when (and (vregular? v1) (vregular? v2) (in-bounds? v1))
         (cond [(in-bounds? v2)
                (add-shape!
-                (arrow-glyph alpha c v1 v2 (->brush-color (plot-background)) (+ 2 pen-width) 'solid))
+                (arrow-glyph alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
+                             (->brush-color (plot-background)) (+ 2 pen-width) 'solid))
                (add-shape!
-                (arrow-glyph alpha c v1 v2 pen-color pen-width pen-style))]
+                (arrow-glyph alpha (plot->norm c) (plot->norm v1) (plot->norm v2)
+                             pen-color pen-width pen-style))]
               [else  (put-line v1 v2)])))
     
     (define/public (put-tick v radius angle)
       (when (and (vregular? v) (in-bounds? v))
-        (add-shape! (tick-glyph alpha v radius angle pen-color pen-width pen-style))))
+        (add-shape! (tick-glyph alpha (plot->norm v) radius angle
+                                pen-color pen-width pen-style))))
     )) ; end class
diff --git a/collects/plot/plot3d/plot.rkt b/collects/plot/plot3d/plot.rkt
index 5617724200..3e660596a0 100644
--- a/collects/plot/plot3d/plot.rkt
+++ b/collects/plot/plot3d/plot.rkt
@@ -36,14 +36,12 @@
 (define (get-bounds-rect renderer-list x-min x-max y-min y-max z-min z-max)
   (define given-bounds-rect (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)))
   (define plot-bounds-rect (bounds-fixpoint renderer-list given-bounds-rect))
-  
   (when (or (not (rect-regular? plot-bounds-rect))
             (rect-zero-area? plot-bounds-rect))
     (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max)) plot-bounds-rect)
     (error 'plot "could not determine sensible plot bounds; got x ∈ [~a,~a], y ∈ [~a,~a], z ∈ [~a,~a]"
            x-min x-max y-min y-max z-min z-max))
-  
-  (rect-inexact->exact plot-bounds-rect))
+  plot-bounds-rect)
 
 (define (get-ticks renderer-list bounds-rect)
   (define-values (all-x-ticks all-x-far-ticks all-y-ticks all-y-far-ticks all-z-ticks all-z-far-ticks)
@@ -56,7 +54,6 @@
       (define ticks-fun (plot-element-ticks-fun r))
       (cond [ticks-fun  (ticks-fun bounds-rect)]
             [else       (values empty empty empty empty empty empty)])))
-  
   (values (remove-duplicates (append* all-x-ticks))
           (remove-duplicates (append* all-x-far-ticks))
           (remove-duplicates (append* all-y-ticks))
diff --git a/collects/plot/plot3d/surface.rkt b/collects/plot/plot3d/surface.rkt
index f5ea70358f..28da7c5ba9 100644
--- a/collects/plot/plot3d/surface.rkt
+++ b/collects/plot/plot3d/surface.rkt
@@ -11,17 +11,20 @@
 
 (define ((surface3d-render-proc f samples color style line-color line-width line-style alpha label)
          area)
-  (match-define (vector (ivl x-min x-max) (ivl y-min y-max) z-ivl) (send area get-bounds-rect))
+  (match-define (vector (ivl x-min x-max) (ivl y-min y-max) (ivl z-min z-max))
+    (send area get-bounds-rect))
   (define sample (f x-min x-max (animated-samples samples)
                     y-min y-max (animated-samples samples)))
   
   (send area put-alpha alpha)
   (send area put-brush color style)
   (send area put-pen line-color line-width line-style)
-  (for-2d-sample
-   (xa xb ya yb z1 z2 z3 z4) sample
-   (send area put-polygon
-         (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4))))
+  (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
+     (send area put-polygon
+           (list (vector xa ya z1) (vector xb ya z2) (vector xb yb z3) (vector xa yb z4)))))
   
   (cond [label  (rectangle-legend-entry label color style line-color line-width line-style)]
         [else   empty]))
diff --git a/collects/plot/tests/extreme-bounds-tests.rkt b/collects/plot/tests/extreme-bounds-tests.rkt
new file mode 100644
index 0000000000..6b4b6b62b7
--- /dev/null
+++ b/collects/plot/tests/extreme-bounds-tests.rkt
@@ -0,0 +1,56 @@
+#lang racket
+
+(require plot plot/utils)
+
+(plot-x-label #f)
+(plot-y-label #f)
+
+(define stops (list (* 2 (inexact->exact -max.0))
+                    (inexact->exact -max.0)
+                    (inexact->exact -min.0)
+                    0
+                    (inexact->exact +min.0)
+                    (inexact->exact +max.0)
+                    (* 2 (inexact->exact +max.0))))
+
+(define (extreme-real->string x)
+  (real->plot-label x (digits-for-range 0 (abs x))))
+
+(define num-2d 0)
+(time
+ (for ([x-min  (in-list stops)]
+       [x-max  (in-list (rest stops))]
+       #:when #t
+       [y-min  (in-list stops)]
+       [y-max  (in-list (rest stops))])
+   (displayln
+    (list (format "[~a,~a] × [~a,~a]"
+                  (extreme-real->string x-min) (extreme-real->string x-max)
+                  (extreme-real->string y-min) (extreme-real->string y-max))
+          (time (plot (lines (list (vector x-min y-min) (vector x-max y-max)))))))
+   (set! num-2d (+ 1 num-2d))))
+
+(printf "Total 2D plots: ~a~n" num-2d)
+(printf "-----------------------------------------~n")
+
+(define num-3d 0)
+
+(time
+ (for ([x-min  (in-list stops)]
+       [x-max  (in-list (rest stops))]
+       #:when #t
+       [y-min  (in-list stops)]
+       [y-max  (in-list (rest stops))]
+       #:when #t
+       [z-min  (in-list stops)]
+       [z-max  (in-list (rest stops))])
+   (displayln
+    (list (format "[~a,~a] × [~a,~a] × [~a,~a]"
+                  (extreme-real->string x-min) (extreme-real->string x-max)
+                  (extreme-real->string y-min) (extreme-real->string y-max)
+                  (extreme-real->string z-min) (extreme-real->string z-max))
+          (time (plot3d (lines3d (list (vector x-min y-min z-min) (vector x-max y-max z-max)))))))
+   (set! num-3d (+ 1 num-3d))))
+
+(printf "Total plots: ~a~n" num-3d)
+(printf "-----------------------------------------~n")
diff --git a/collects/plot/tests/isosurface-tests.rkt b/collects/plot/tests/isosurface-tests.rkt
index ea93033e86..e403fff196 100644
--- a/collects/plot/tests/isosurface-tests.rkt
+++ b/collects/plot/tests/isosurface-tests.rkt
@@ -55,3 +55,9 @@
                              #:color "navajowhite" #:width 2)
                (lines3d '(#(0 0 2) #(0 0 -2)) #:color "navajowhite" #:width 2))
          #:title "A Seashell" #:x-label #f #:y-label #f #:angle 210 #:altitude 30))
+
+(define flepsilon (flnext 0.0))
+
+(time
+ (plot3d (isosurface3d (λ (x y z) (+ (sqr x) (sqr y) (sqr z))) (sqr (inexact->exact flepsilon))
+                       (- flepsilon) flepsilon (- flepsilon) flepsilon (- flepsilon) flepsilon)))
diff --git a/collects/plot/tests/plot2d-tests.rkt b/collects/plot/tests/plot2d-tests.rkt
index e278707766..287e55eab9 100644
--- a/collects/plot/tests/plot2d-tests.rkt
+++ b/collects/plot/tests/plot2d-tests.rkt
@@ -18,6 +18,10 @@
 
 (time (plot (function values 0 1000)))
 
+(time
+ (define flepsilon (ordinal->flonum 1))
+ (plot (lines (list (vector 0 0) (vector flepsilon flepsilon)))))
+
 (parameterize ([plot-x-transform  log-transform]
                [plot-x-ticks      (log-ticks #:base 4)]
                [plot-y-transform  log-transform]
diff --git a/collects/plot/tests/plot3d-tests.rkt b/collects/plot/tests/plot3d-tests.rkt
index 418471d766..44771f0cc5 100644
--- a/collects/plot/tests/plot3d-tests.rkt
+++ b/collects/plot/tests/plot3d-tests.rkt
@@ -10,6 +10,15 @@
 (time
  (plot3d (points3d empty) #:x-min -1 #:x-max 1 #:y-min -1 #:y-max 1 #:z-min -1 #:z-max 1))
 
+(define flepsilon (ordinal->flonum 1))
+
+(time
+ (plot3d (lines3d (list (vector 0 0 0) (vector flepsilon flepsilon flepsilon)))))
+
+(time
+ (plot3d (contour-intervals3d (λ (x y) (- (sqr x) (sqr y)))
+                              (- flepsilon) flepsilon (- flepsilon) flepsilon)))
+
 (parameterize ([plot-background  "black"]
                [plot-foreground  "white"]
                [plot-background-alpha  1/2]