diff --git a/collects/math/flonum.rkt b/collects/math/flonum.rkt
index 80d0db1d75..d603c0829d 100644
--- a/collects/math/flonum.rkt
+++ b/collects/math/flonum.rkt
@@ -11,7 +11,7 @@
          "private/flonum/flonum-factorial.rkt"
          "private/flonum/flonum-log1pmx.rkt"
          "private/flonum/flonum-polyfun.rkt"
-         "private/flonum/flonum-syntax.rkt"
+         "private/flonum/flonum-error.rkt"
          "private/flonum/expansion/expansion-base.rkt"
          "private/flonum/expansion/expansion-exp.rkt"
          "private/flonum/expansion/expansion-log.rkt"
@@ -28,7 +28,7 @@
           "private/flonum/flonum-factorial.rkt"
           "private/flonum/flonum-log1pmx.rkt"
           "private/flonum/flonum-polyfun.rkt"
-          "private/flonum/flonum-syntax.rkt"
+          "private/flonum/flonum-error.rkt"
           "private/flonum/expansion/expansion-base.rkt"
           "private/flonum/expansion/expansion-exp.rkt"
           "private/flonum/expansion/expansion-log.rkt"
diff --git a/collects/math/private/bigfloat/bigfloat-mpfr.rkt b/collects/math/private/bigfloat/bigfloat-mpfr.rkt
index 6c12131692..560b7d2dc9 100644
--- a/collects/math/private/bigfloat/bigfloat-mpfr.rkt
+++ b/collects/math/private/bigfloat/bigfloat-mpfr.rkt
@@ -1,6 +1,7 @@
 #lang typed/racket/base
 
-(require (only-in "mpfr.rkt" 1ary-funs 1ary-preds 1ary2-funs 2ary-funs)
+(require racket/flonum
+         (only-in "mpfr.rkt" 1ary-funs 1ary-preds 1ary2-funs 2ary-funs)
          "../base/base-random.rkt"
          "utils.rkt")
 
@@ -139,11 +140,18 @@
 (: bigfloat->fl2 (Bigfloat -> (Values Flonum Flonum)))
 (define (bigfloat->fl2 x)
   (define x2 (bigfloat->flonum x))
-  (values x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2)))))
+  (cond [(rational? x2)
+         (let ([x2  (+ x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2))))])
+           (cond [(rational? x2)
+                  (values x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2))))]
+                 [else
+                  (values x2 0.0)]))]
+        [else  (values x2 0.0)]))
 
 (: fl2->bigfloat (Flonum Flonum -> Bigfloat))
 (define (fl2->bigfloat x2 x1)
-  (bf+ (flonum->bigfloat x1) (flonum->bigfloat x2)))
+  (cond [(fl= x1 0.0)  (bf x2)]
+        [else  (bf+ (flonum->bigfloat x1) (flonum->bigfloat x2))]))
 
 (provide
  ;; Parameters
diff --git a/collects/math/private/flonum/expansion/expansion-base.rkt b/collects/math/private/flonum/expansion/expansion-base.rkt
index 7c0a244aab..b597b9785f 100644
--- a/collects/math/private/flonum/expansion/expansion-base.rkt
+++ b/collects/math/private/flonum/expansion/expansion-base.rkt
@@ -6,81 +6,124 @@ Arithmetic based on:
 Jonathan Richard Shewchuk
 Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates
 Discrete & Computational Geometry 18(3):305–363, October 1997
+
+Other parts shamelessly stolen from crlibm (which is LGPL)
 |#
 
 (require "../flonum-functions.rkt"
-         "../flonum-syntax.rkt")
+         "../flonum-bits.rkt"
+         "../flonum-error.rkt"
+         "../flonum-constants.rkt"
+         "../utils.rkt")
 
-(provide fl2 fl2->real
-         fl2+ fl2- fl2*split-fl fl2* fl2/
-         flsqrt/error fl2sqrt)
+(provide fl2? fl2 fl2->real
+         fl2ulp fl2ulp-error
+         +max-fl2-subnormal.0 -max-fl2-subnormal.0
+         fl2abs fl2+ fl2- fl2*split-fl fl2* fl2sqr fl2/
+         fl2sqrt flsqrt/error)
+
+(: floverlapping? (Flonum Flonum -> Boolean))
+(define (floverlapping? x2 x1)
+  (define-values (s2 e2) (flonum->sig+exp (flabs x2)))
+  (define-values (s1 e1) (flonum->sig+exp (flabs x1)))
+  (define-values (n1 n2)
+    (if (e2 . > . e1)
+        (values s1 (arithmetic-shift s2 (- e2 e1)))
+        (values (arithmetic-shift s1 (- e1 e2)) s2)))
+  (not (= (bitwise-ior n1 n2)
+          (bitwise-xor n1 n2))))
+
+(: fl2? (Flonum Flonum -> Boolean))
+(define (fl2? x2 x1)
+  (cond [(flrational? x2)
+         (cond [(flrational? x1)
+                (cond [((flabs x2) . < . (flabs x1))  #f]
+                      [else  (not (floverlapping? x2 x1))])]
+               [else  #f])]
+        [else  (fl= x1 0.0)]))
 
 ;; ===================================================================================================
 ;; Conversion
 
-(: fl2 (case-> (Real -> (Values Flonum Flonum))
-               (Flonum Flonum -> (Values Flonum Flonum))))
-(define fl2
-  (case-lambda
-    [(x)
-     (cond [(flonum? x)  (values x 0.0)]
-           [(single-flonum? x)  (values (fl x) 0.0)]
-           [else  (let ([x2  (fl x)])
-                    (values x2 (fl (- x (inexact->exact x2)))))])]
-    [(x y)
-     (fast-fl+/error x y)]))
+(: fl2 (Real -> (Values Flonum Flonum)))
+(define (fl2 x)
+  (cond [(flonum? x)  (values x 0.0)]
+        [(single-flonum? x)  (values (fl x) 0.0)]
+        [else  (let ([x2  (fl x)])
+                 (cond [(flrational? x2)
+                        (let ([x2  (fl+ x2 (fl (- x (inexact->exact x2))))])
+                          (cond [(flrational? x2)
+                                 (values x2 (fl (- x (inexact->exact x2))))]
+                                [else  (values x2 0.0)]))]
+                       [else  (values x2 0.0)]))]))
 
 (: fl2->real (Flonum Flonum -> Real))
 (define (fl2->real x2 x1)
-  (cond [(and (x1 . fl> . -inf.0) (x1 . fl< . +inf.0)
-              (x2 . fl> . -inf.0) (x2 . fl< . +inf.0))
+  (cond [(and (flrational? x2) (flrational? x1))
          (+ (inexact->exact x2) (inexact->exact x1))]
         [else  (fl+ x1 x2)]))
 
-(: fl3->fl2 (Flonum Flonum Flonum -> (Values Flonum Flonum)))
-(define (fl3->fl2 e3 e2 e1)
-  (values e3 (fl+ e2 e1)))
-
 (: fl4->fl2 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)))
 (define (fl4->fl2 e4 e3 e2 e1)
   (values e4 (fl+ e3 (fl+ e2 e1))))
 
+;; ===================================================================================================
+;; Error
+
+(: fl2ulp (Flonum Flonum -> Flonum))
+(define (fl2ulp x2 x1)
+  (cond [(fl= x2 0.0)  0.0]
+        [else  (flmax +min.0 (fl* (flulp x2) (flexpt 2.0 -52.0)))]))
+
+(: fl2ulp-error (Flonum Flonum Real -> Flonum))
+(define (fl2ulp-error x2 x1 r)
+  (define x (fl2->real x2 x1))
+  (define-values (r2 r1) (fl2 r))
+  (cond [(eqv? x r)  0.0]
+        [(and (fl= x2 0.0) (fl= r2 0.0))  0.0]
+        [(and (fl= x2 +inf.0) (fl= r2 +inf.0))  0.0]
+        [(and (fl= x2 -inf.0) (fl= r2 -inf.0))  0.0]
+        [(zero? r)  +inf.0]
+        [(and (rational? x) (flrational? r2))
+         (flabs (fl (/ (- (inexact->exact x) (inexact->exact r))
+                       (inexact->exact (flmax +min.0 (fl2ulp r2 r1))))))]
+        [else  +inf.0]))
+
+(define +max-fl2-subnormal.0 1.0020841800044863e-292)  ; (flprev (flexpt 2.0 -970.0))
+(define -max-fl2-subnormal.0 (- +max-fl2-subnormal.0))
+
+;; ===================================================================================================
+;; Absolute value
+
+(: fl2abs (case-> (Flonum -> (Values Flonum Flonum))
+                  (Flonum Flonum -> (Values Flonum Flonum))))
+(define (fl2abs x2 [x1 0.0])
+  (define x (fl+ x2 x1))
+  (cond [(not (flrational? x))  (values (flabs x2) 0.0)]
+        [(fl= x 0.0)  (values 0.0 0.0)]
+        [(fl> x 0.0)  (values x2 x1)]
+        [else  (values (- x2) (- x1))]))
+
 ;; ===================================================================================================
 ;; Addition and subtraction
 
-(: raw-fl2+fl (Flonum Flonum Flonum -> (Values Flonum Flonum Flonum)))
-(define (raw-fl2+fl e2 e1 b)
-  (let*-values ([(Q h1)  (fast-fl+/error b e1)]
-                [(h3 h2)  (fast-fl+/error Q e2)])
-    (values h3 h2 h1)))
-
-(: raw-fl2+ (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum Flonum Flonum)))
-(define (raw-fl2+ e2 e1 f2 f1)
-  (let*-values ([(h3 h2 h1)  (raw-fl2+fl e2 e1 f1)]
-                [(h4 h3 h2)  (raw-fl2+fl h3 h2 f2)])
-    (values h4 h3 h2 h1)))
-
 (: fl2+ (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
                 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
-(define fl2+
-  (case-lambda
-    [(e2 e1 b)
-     (let-values ([(h3 h2 h1)  (raw-fl2+fl e2 e1 b)])
-       (fl3->fl2 h3 h2 h1))]
-    [(x2 x1 y2 y1)
-     (let*-values ([(e4 e3 e2 e1)  (raw-fl2+ x2 x1 y2 y1)])
-       (fl4->fl2 e4 e3 e2 e1))]))
+(define (fl2+ x2 x1 y2 [y1 0.0])
+  (define r (fl+ x2 y2))
+  (cond [(not (flrational? r))  (values r 0.0)]
+        [(and (fl= x2 0.0) (fl= y2 0.0))  (values r 0.0)]
+        [else
+         (define s (if ((flabs x2) . fl> . (flabs y2))
+                       (fl+ (fl+ (fl+ (fl- x2 r) y2) y1) x1)
+                       (fl+ (fl+ (fl+ (fl- y2 r) x2) x1) y1)))
+         (define z2 (fl+ r s))
+         (values z2 (fl+ (fl- r z2) s))]))
 
 (: fl2- (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
                 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
-(define fl2-
-  (case-lambda
-    [(e2 e1 b)
-     (let-values ([(h3 h2 h1)  (raw-fl2+fl e2 e1 (- b))])
-       (fl3->fl2 h3 h2 h1))]
-    [(x2 x1 y2 y1)
-     (let*-values ([(e4 e3 e2 e1)  (raw-fl2+ x2 x1 (- y2) (- y1))])
-       (fl4->fl2 e4 e3 e2 e1))]))
+(define (fl2- x2 x1 y2 [y1 0.0])
+  (fl2+ x2 x1 (- y2) (- y1)))
 
 ;; ===================================================================================================
 ;; Multiplication and division
@@ -119,56 +162,104 @@ Discrete & Computational Geometry 18(3):305–363, October 1997
 
 (: fl2* (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
                 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
-(define fl2*
+(define (fl2* x2 x1 y2 [y1 0.0])
+  (define z (fl* x2 y2))
+  (cond [(and (flrational? z) (not (fl= z 0.0)) (not (flsubnormal? z)))
+         (define dx (near-pow2 x2))
+         (define dy (near-pow2 y2))
+         (define d (fl* dx dy))
+         (let* ([x2  (fl/ x2 dx)]
+                [x1  (fl/ x1 dx)]
+                [y2  (fl/ y2 dy)]
+                [y1  (fl/ y1 dy)]
+                [up  (fl* x2 (fl+ 1.0 (flexpt 2.0 27.0)))]
+                [vp  (fl* y2 (fl+ 1.0 (flexpt 2.0 27.0)))]
+                [u1  (fl+ (fl- x2 up) up)]
+                [v1  (fl+ (fl- y2 vp) vp)]
+                [u2  (fl- x2 u1)]
+                [v2  (fl- y2 v1)]
+                [m2  (fl* x2 y2)]
+                [m1  (fl+ (fl+ (fl+ (fl+ (fl+ (fl- (fl* u1 v1) m2)
+                                              (fl* u1 v2))
+                                         (fl* u2 v1))
+                                    (fl* u2 v2))
+                               (fl* x2 y1))
+                          (fl* x1 y2))]
+                [z2  (fl+ m2 m1)]
+                [z1  (fl+ (fl- m2 z2) m1)]
+                [z2  (fl* d z2)])
+           (values z2 (if (flrational? z2) (fl* d z1) 0.0)))]
+        [else
+         (values z 0.0)]))
+
+(: fl2sqr (case-> (Flonum -> (Values Flonum Flonum))
+                  (Flonum Flonum -> (Values Flonum Flonum))))
+;; Derived from fl2*
+(define fl2sqr
   (case-lambda
-    [(e2 e1 b)
-     (let-values ([(b-hi b-lo)  (flsplit b)])
-       (fl2*split-fl e2 e1 b-hi b-lo))]
-    [(x2 x1 y2 y1)
-     (let*-values ([(x2-hi x2-lo)  (flsplit x2)]
-                   [(x1-hi x1-lo)  (flsplit x1)]
-                   [(y2-hi y2-lo)  (flsplit y2)]
-                   [(y1-hi y1-lo)  (flsplit y1)]
-                   [(a2 a1)  (split-fl2*split-fl x2-hi x2-lo x1-hi x1-lo y1-hi y1-lo)]
-                   [(b2 b1)  (split-fl2*split-fl x2-hi x2-lo x1-hi x1-lo y2-hi y2-lo)])
-       (fl2+ a2 a1 b2 b1))]))
+    [(x)  (flsqr/error x)]
+    [(x2 x1)
+     (define z (fl* x2 x2))
+     (cond [(and (flrational? z) (not (fl= z 0.0)) (not (flsubnormal? z)))
+            (define d (near-pow2 x2))
+            (define d^2 (fl* d d))
+            (let* ([x2  (fl/ x2 d)]
+                   [x1  (fl/ x1 d)]
+                   [up  (fl* x2 (fl+ 1.0 (flexpt 2.0 27.0)))]
+                   [u1  (fl+ (fl- x2 up) up)]
+                   [u2  (fl- x2 u1)]
+                   [m2  (fl* x2 x2)]
+                   [m1  (fl+ (fl+ (fl+ (fl- (fl* u1 u1) m2)
+                                       (fl* 2.0 (fl* u1 u2)))
+                                  (fl* u2 u2))
+                             (fl* 2.0 (fl* x2 x1)))]
+                   [z2  (fl+ m2 m1)]
+                   [z1  (fl+ (fl- m2 z2) m1)]
+                   [z2  (fl* d^2 z2)])
+              (values z2 (if (flrational? z2) (fl* d^2 z1) 0.0)))]
+           [else
+            (values z 0.0)])]))
 
 (: fl2/ (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
                 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
-(define fl2/
-  (case-lambda
-    [(x2 x1 y)
-     (let*-values ([(a2 a1)  (fast-fl//error x1 y)]
-                   [(b2 b1)  (fast-fl//error x2 y)])
-       (fl2+ a2 a1 b2 b1))]
-    [(x2 x1 y2 y1)
-     ;; Compute three "digits" (flonums) of two-flonum long division; the third ensures the result is
-     ;; correctly rounded
-     (let*-values ([(z2)  (fl/ x2 y2)]
-                   [(w2 w1)  (fl2* y2 y1 z2)]
-                   [(x2 x1)  (fl2- x2 x1 w2 w1)]
-                   [(z1)  (fl/ x2 y2)]
-                   [(w2 w1)  (fl2* y2 y1 z1)]
-                   [(x2 x1)  (fl2- x2 x1 w2 w1)])
-       (fl3->fl2 z2 z1 (/ x2 y2)))]))
+(define (fl2/ x2 x1 y2 [y1 0.0])
+  (define z (fl/ x2 y2))
+  (cond [(and (flrational? z) (not (fl= z 0.0)) (flrational? y2))
+         (define d (near-pow2/div x2 y2))
+         (let*-values ([(x2 x1)  (values (fl/ x2 d) (fl/ x1 d))]
+                       [(y2 y1)  (values (fl/ y2 d) (fl/ y1 d))]
+                       [(c2)  (fl/ x2 y2)]
+                       [(u2 u1)  (fl*/error c2 y2)]
+                       [(c1)  (fl/ (fl- (fl+ (fl- (fl- x2 u2) u1) x1) (fl* c2 y1)) y2)]
+                       [(z2)  (fl+ c2 c1)])
+           (values z2 (if (flrational? z2) (fl+ (fl- c2 z2) c1) 0.0)))]
+        [else
+         (values z 0.0)]))
 
 ;; ===================================================================================================
 ;; Square roots
 
-(: flsqrt/error (Flonum -> (Values Flonum Flonum)))
+(: fl2sqrt (case-> (Flonum -> (Values Flonum Flonum))
+                   (Flonum Flonum -> (Values Flonum Flonum))))
 ;; One-flonum estimate followed by one Newton's method iteration
-;; This could be a little faster if `y' were split only once
-(define (flsqrt/error x)
-  (let*-values ([(y)  (flsqrt x)]
-                [(z2 z1)  (fast-flsqr/error y)]
-                [(dy2 dy1)  (fl2+ (- z2) (- z1) x)]
-                [(dy2 dy1)  (fl2/ dy2 dy1 y)])
-    (fl2+ (* 0.5 dy2) (* 0.5 dy1) y)))
+(define (fl2sqrt x2 [x1 0.0])
+  (cond [(and (flrational? x2) (not (fl= x2 0.0)))
+         (define-values (d^2 d)
+           (cond [(x2 . fl<= . +max-fl2-subnormal.0)  (values (flexpt 2.0 -104.0)
+                                                              (flexpt 2.0 -52.0))]
+                 [(x2 . fl> . 1e300)  (values (flexpt 2.0 104.0)
+                                              (flexpt 2.0 52.0))]
+                 [else  (values 1.0 1.0)]))
+         (let*-values ([(x2 x1)  (values (fl/ x2 d^2) (fl/ x1 d^2))]
+                       [(y)  (flsqrt (fl+ x2 x1))]
+                       [(z2 z1)  (fast-flsqr/error y)]
+                       [(dy2 dy1)  (fl2- x2 x1 z2 z1)]
+                       [(dy2 dy1)  (fl2/ dy2 dy1 y)]
+                       [(y2 y1)  (fl2+ (fl* 0.5 dy2) (fl* 0.5 dy1) y)]
+                       [(y2)  (fl* y2 d)])
+           (values y2 (if (flrational? y2) (fl* y1 d) 0.0)))]
+        [else
+         (values (flsqrt x2) 0.0)]))
 
-(: fl2sqrt (Flonum Flonum -> (Values Flonum Flonum)))
-(define (fl2sqrt x2 x1)
-  (let*-values ([(y)  (flsqrt (fl+ x1 x2))]
-                [(z2 z1)  (fast-flsqr/error y)]
-                [(dy2 dy1)  (fl2- x2 x1 z2 z1)]
-                [(dy2 dy1)  (fl2/ dy2 dy1 y)])
-    (fl2+ (* 0.5 dy2) (* 0.5 dy1) y)))
+(: flsqrt/error (Flonum -> (Values Flonum Flonum)))
+(define (flsqrt/error x) (fl2sqrt x 0.0))
diff --git a/collects/math/private/flonum/expansion/expansion-exp.rkt b/collects/math/private/flonum/expansion/expansion-exp.rkt
index ddf38a2309..a2df696b2b 100644
--- a/collects/math/private/flonum/expansion/expansion-exp.rkt
+++ b/collects/math/private/flonum/expansion/expansion-exp.rkt
@@ -17,7 +17,7 @@ Note to self: exp(x+y)-1 = (exp(x)-1) * (exp(y)-1) + (exp(x)-1) + (exp(y)-1)
          "../../../base.rkt"
          "../../../bigfloat.rkt"
          "../flonum-functions.rkt"
-         "../flonum-syntax.rkt"
+         "../flonum-error.rkt"
          "../flonum-exp.rkt"
          "../flonum-constants.rkt"
          "expansion-base.rkt"
diff --git a/collects/math/private/flonum/expansion/expansion-log.rkt b/collects/math/private/flonum/expansion/expansion-log.rkt
index fe41ce1504..19e697f850 100644
--- a/collects/math/private/flonum/expansion/expansion-log.rkt
+++ b/collects/math/private/flonum/expansion/expansion-log.rkt
@@ -8,7 +8,7 @@ perform one Newton step.
 |#
 
 (require "../flonum-functions.rkt"
-         "../flonum-syntax.rkt"
+         "../flonum-error.rkt"
          "../flonum-log.rkt"
          "expansion-base.rkt"
          "expansion-exp.rkt")
diff --git a/collects/math/private/flonum/flonum-bits.rkt b/collects/math/private/flonum/flonum-bits.rkt
index f4f6a74448..691ed4b7c1 100644
--- a/collects/math/private/flonum/flonum-bits.rkt
+++ b/collects/math/private/flonum/flonum-bits.rkt
@@ -4,6 +4,8 @@
          racket/performance-hint)
 
 (provide flonum->bit-field bit-field->flonum
+         flonum->fields fields->flonum
+         flonum->sig+exp sig+exp->flonum
          flonum->ordinal ordinal->flonum
          flstep flnext flprev flonums-between
          flulp)
@@ -20,6 +22,48 @@
         [else
          (raise-argument-error 'bit-field->flonum "Integer in [0 .. #xffffffffffffffff]" i)]))
 
+(define implicit-leading-one (arithmetic-shift 1 52))
+(define max-significand (- implicit-leading-one 1))
+(define max-exponent 2047)
+(define max-signed-exponent 1023)
+(define min-signed-exponent -1074)
+
+(: flonum->fields (Flonum -> (Values (U 0 1) Index Natural)))
+(define (flonum->fields x)
+  (define n (flonum->bit-field x))
+  (values (if (zero? (bitwise-bit-field n 63 64)) 0 1)
+          (assert (bitwise-bit-field n 52 63) index?)
+          (bitwise-bit-field n 0 52)))
+
+(: fields->flonum (Integer Integer Integer -> Flonum))
+(define (fields->flonum s e m)
+  (cond [(not (or (= s 0) (= s 1)))
+         (raise-argument-error 'fields->flonum "(U 0 1)" 0 s e m)]
+        [(or (e . < . 0) (e . > . max-exponent))
+         (raise-argument-error 'fields->flonum (format "Natural <= ~e" max-exponent) 1 s e m)]
+        [(or (m . < . 0) (m . > . max-significand))
+         (raise-argument-error 'fields->flonum (format "Natural <= ~a" max-significand) 2 s e m)]
+        [else
+         (bit-field->flonum (bitwise-ior (arithmetic-shift s 63)
+                                         (arithmetic-shift e 52)
+                                         m))]))
+
+(: flonum->sig+exp (Flonum -> (Values Integer Fixnum)))
+(define (flonum->sig+exp x)
+  (define-values (s e m) (flonum->fields x))
+  (let-values ([(sig exp)  (if (= e 0)
+                               (values m -1074)
+                               (values (bitwise-ior m implicit-leading-one)
+                                       (assert (- e 1075) fixnum?)))])
+    (values (if (zero? s) sig (- sig)) exp)))
+
+(: sig+exp->flonum (Integer Integer -> Flonum))
+(define (sig+exp->flonum sig exp)
+  (cond [(= sig 0)  0.0]
+        [(exp . > . max-signed-exponent)  (if (sig . < . 0) -inf.0 +inf.0)]
+        [(exp . < . min-signed-exponent)  (if (sig . < . 0) -0.0 0.0)]
+        [else  (real->double-flonum (* sig (expt 2 exp)))]))
+
 (: flonum->ordinal (Flonum -> Integer))
 (define (flonum->ordinal x)
   (cond [(x . fl< . 0.0)  (- (flonum->bit-field (fl- 0.0 x)))]
diff --git a/collects/math/private/flonum/flonum-error.rkt b/collects/math/private/flonum/flonum-error.rkt
new file mode 100644
index 0000000000..45b4cc810a
--- /dev/null
+++ b/collects/math/private/flonum/flonum-error.rkt
@@ -0,0 +1,184 @@
+#lang racket/base
+
+(provide flsplit
+         fast-mono-fl+/error
+         fast-mono-fl-/error
+         fast-fl+/error
+         fast-fl-/error
+         fast-fl*/error
+         fast-flsqr/error
+         fast-fl//error
+         fl+/error
+         fl-/error
+         fl*/error
+         flsqr/error
+         fl//error)
+
+(module untyped-defs racket/base
+  (require (for-syntax racket/base)
+           "flonum-functions.rkt")
+  
+  (provide (all-defined-out))
+  
+  ;(: flsplit (Flonum -> (Values Flonum Flonum)))
+  ;; Splits a flonum into a two flonums `hi' and `lo' with 26 bits precision each, such that
+  ;; |hi| >= |lo| and hi + lo = a. (The extra sign bit accounts for the missing bit.)
+  ;; This function returns (values +nan.0 +nan.0) for |a| >= 1.3393857490036326e+300.
+  (define-syntax-rule (flsplit a-expr)
+    (let ([a a-expr])
+      (let* ([c   (fl* a (fl+ 1.0 (flexpt 2.0 27.0)))]
+             [x2  (fl- c (fl- c a))])
+        (values x2 (fl- a x2)))))
+  
+  ;; =================================================================================================
+  ;; Fast monotone addition and subtraction
+  
+  ;(: fast-mono-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a+b and its rounding error
+  ;; Assumes |a| >= |b|
+  (define-syntax-rule (fast-mono-fl+/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let ([x2  (+ a b)])
+        (values x2 (- b (- x2 a))))))
+  
+  ;(: fast-mono-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a+b and its rounding error
+  ;; Assumes |a| >= |b|
+  (define-syntax-rule (fast-mono-fl-/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let ([x2  (- a b)])
+        (values x2 (- (- a x2) b)))))
+  
+  ;; =================================================================================================
+  ;; Fast arithmetic that returns rounding error
+  
+  ;(: fast-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a+b and its rounding error
+  (define-syntax-rule (fast-fl+/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let* ([x2  (fl+ a b)]
+             [v   (fl- x2 a)])
+        (values x2 (fl+ (fl- a (fl- x2 v)) (fl- b v))))))
+  
+  ;(: fast-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a-b and its rounding error
+  (define-syntax-rule (fast-fl-/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let* ([x2  (fl- a b)]
+             [v   (fl- x2 a)])
+        (values x2 (fl- (fl- a (fl- x2 v)) (fl+ b v))))))
+  
+  ;(: fast-fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a*b and its rounding error
+  (define-syntax-rule (fast-fl*/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let*-values ([(x2)  (fl* a b)]
+                    [(a2 a1)  (flsplit a)]
+                    [(b2 b1)  (flsplit b)])
+        (values x2 (- (fl- (fl- (fl- (fl- x2 (fl* a2 b2))
+                                     (fl* a1 b2))
+                                (fl* a2 b1))
+                           (fl* a1 b1)))))))
+  
+  ;(: fast-flfma/error (Flonum Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a*b+c and its rounding error
+  (define-syntax-rule (fast-flfma/error a-expr b-expr c-expr)
+    (let*-values ([(y2 y1)  (fast-fl*/error a-expr b-expr)]
+                  [(h0 h1)   (fast-fl+/error c-expr y1)]
+                  [(h3 h2)  (fast-fl+/error h0 y2)])
+      (values h3 (fl+ h2 h1))))
+  
+  #;; If we had hardware fused multiply-add:
+  (define-syntax-rule (fast-fl*/error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let ([x2  (fl* a b)])
+        (values x2 (flfma a b (- x2))))))
+
+  ;(: fast-flsqr/error (Flonum -> (Values Flonum Flonum)))
+  ;; Returns a*a and its rounding error
+  (define-syntax-rule (fast-flsqr/error a-expr)
+    (let ([a a-expr])
+      (let*-values ([(x2)  (fl* a a)]
+                    [(a2 a1)  (flsplit a)])
+        (values x2 (- (fl- (fl- (fl- x2 (fl* a2 a2))
+                                (fl* 2.0 (fl* a2 a1)))
+                           (fl* a1 a1)))))))
+  
+  ;(: fast-fl//error (Flonum Flonum -> (Values Flonum Flonum)))
+  ;; Returns a/b and its rounding error
+  (define-syntax-rule (fast-fl//error a-expr b-expr)
+    (let ([a a-expr] [b b-expr])
+      (let*-values ([(x2)  (fl/ a b)]
+                    [(w2 w1)  (fast-fl*/error x2 b)])
+        (fast-mono-fl+/error x2 (fl/ (fl- (fl- a w2) w1) b)))))
+  
+  )  ; module
+
+(require (submod "." untyped-defs))
+
+(module typed-defs typed/racket/base
+  (require racket/performance-hint
+           (submod ".." untyped-defs)
+           "flonum-functions.rkt"
+           "utils.rkt")
+  
+  (provide (all-defined-out))
+  
+  ;; =================================================================================================
+  ;; Function versions of the above that are well-defined for the largest domain, and return 0.0 as
+  ;; the second argument whenever the first isn't rational
+  
+  (begin-encourage-inline
+    
+    (: fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
+    (define (fl+/error a b)
+      (let-values ([(x2 x1)  (fast-fl+/error a b)])
+        (values x2 (if (flrational? x2) x1 0.0))))
+    
+    (: fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
+    (define (fl-/error a b)
+      (let-values ([(x2 x1)  (fast-fl-/error a b)])
+        (values x2 (if (flrational? x2) x1 0.0))))
+    
+    (: fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
+    (define (fl*/error a b)
+      (let ([x2  (fl* a b)])
+        (values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
+                       (let*-values ([(da db)  (values (near-pow2 a) (near-pow2 b))]
+                                     [(d)  (fl* da db)]
+                                     [(a2 a1)  (flsplit (fl/ a da))]
+                                     [(b2 b1)  (flsplit (fl/ b db))])
+                         (fl* d (- (fl- (fl- (fl- (fl- (fl/ x2 d) (fl* a2 b2))
+                                                  (fl* a1 b2))
+                                             (fl* a2 b1))
+                                        (fl* a1 b1)))))
+                       0.0))))
+    
+    (: flsqr/error (Flonum -> (Values Flonum Flonum)))
+    (define (flsqr/error a)
+      (let ([x2  (fl* a a)])
+        (values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
+                       (let*-values ([(d)  (near-pow2 a)]
+                                     [(d^2)  (fl* d d)]
+                                     [(a2 a1)  (flsplit (fl/ a d))])
+                         (fl* d^2 (- (fl- (fl- (fl- (fl/ x2 d^2) (fl* a2 a2))
+                                               (fl* 2.0 (fl* a1 a2)))
+                                          (fl* a1 a1)))))
+                       0.0))))
+    
+    (: fl//error (Flonum Flonum -> (Values Flonum Flonum)))
+    (define (fl//error a b)
+      (let ([x2  (fl/ a b)])
+        (values x2 (if (and (flrational? x2) (flrational? b))
+                       (let* ([d  (near-pow2/div a b)]
+                              [a  (fl/ a d)]
+                              [b  (fl/ b d)])
+                         (let-values ([(w2 w1)  (fl*/error x2 b)])
+                           (fl/ (fl- (fl- a w2) w1) b)))
+                       0.0))))
+    
+    )  ; begin-encourage-inline
+  
+  )  ; module
+
+(require (submod "." typed-defs))
diff --git a/collects/math/private/flonum/flonum-factorial.rkt b/collects/math/private/flonum/flonum-factorial.rkt
index 4f7d9fa9be..f33fbf5911 100644
--- a/collects/math/private/flonum/flonum-factorial.rkt
+++ b/collects/math/private/flonum/flonum-factorial.rkt
@@ -8,7 +8,7 @@
          "flonum-functions.rkt"
          "flonum-log.rkt"
          "flonum-more-functions.rkt"
-         "flonum-syntax.rkt")
+         "flonum-error.rkt")
 
 (provide flfactorial
          flbinomial
diff --git a/collects/math/private/flonum/flonum-functions.rkt b/collects/math/private/flonum/flonum-functions.rkt
index cb1361579a..0fafa850e4 100644
--- a/collects/math/private/flonum/flonum-functions.rkt
+++ b/collects/math/private/flonum/flonum-functions.rkt
@@ -83,13 +83,15 @@
 
 (: flulp-error (Flonum Real -> Flonum))
 (define (flulp-error x r)
+  (define r.0 (fl r))
   (cond [(eqv? x r)  0.0]
-        [(and (fl= x 0.0) (zero? r))  0.0]
+        [(and (fl= x 0.0) (fl= r.0 0.0))  0.0]
+        [(and (fl= x +inf.0) (fl= r.0 +inf.0))  0.0]
+        [(and (fl= x -inf.0) (fl= r.0 -inf.0))  0.0]
         [(zero? r)  +inf.0]
-        [(and (rational? x) (rational? r))
-         (flabs (real->double-flonum
-                 (/ (- (inexact->exact x) (inexact->exact r))
-                    (inexact->exact (flulp x)))))]
+        [(and (flrational? x) (flrational? r.0))
+         (flabs (fl (/ (- (inexact->exact x) (inexact->exact r))
+                       (inexact->exact (flmax +min.0 (flulp r.0))))))]
         [else  +inf.0]))
 
 ;; ===================================================================================================
diff --git a/collects/math/private/flonum/flonum-log.rkt b/collects/math/private/flonum/flonum-log.rkt
index 1ca75ccdc9..a3290b804d 100644
--- a/collects/math/private/flonum/flonum-log.rkt
+++ b/collects/math/private/flonum/flonum-log.rkt
@@ -5,7 +5,7 @@
          "flonum-functions.rkt"
          "flonum-constants.rkt"
          "flonum-exp.rkt"
-         "flonum-syntax.rkt"
+         "flonum-error.rkt"
          "flvector.rkt")
 
 (provide fllog1p fllog+
diff --git a/collects/math/private/flonum/flonum-more-functions.rkt b/collects/math/private/flonum/flonum-more-functions.rkt
index a8d21903e8..5cc3fdf91b 100644
--- a/collects/math/private/flonum/flonum-more-functions.rkt
+++ b/collects/math/private/flonum/flonum-more-functions.rkt
@@ -5,7 +5,7 @@
          "flonum-constants.rkt"
          "flonum-exp.rkt"
          "flonum-log.rkt"
-         "flonum-syntax.rkt"
+         "flonum-error.rkt"
          "flvector.rkt")
 
 (provide flsqrt1pm1
diff --git a/collects/math/private/flonum/flonum-syntax.rkt b/collects/math/private/flonum/flonum-syntax.rkt
deleted file mode 100644
index c3a2d67184..0000000000
--- a/collects/math/private/flonum/flonum-syntax.rkt
+++ /dev/null
@@ -1,105 +0,0 @@
-#lang racket/base
-
-(require (for-syntax racket/base)
-         racket/flonum)
-
-(provide flsplit
-         fast-mono-fl+/error fast-fl+/error fl+/error
-         fast-mono-fl-/error fast-fl-/error fl-/error
-         fast-fl*/error fast-flsqr/error fl*/error
-         fast-fl//error)
-
-;(: flsplit (Flonum -> (Values Flonum Flonum)))
-;; Splits a flonum into a two flonums `hi' and `lo' with 26 bits precision each, such that
-;; |hi| >= |lo| and hi + lo = a. (The extra sign bit accounts for the missing bit.)
-(define-syntax-rule (flsplit a-expr)
-  (let ([a a-expr])
-    (let* ([s  (if ((flabs a) . fl< . 1e300) 1.0 (flexpt 2.0 52.0))]
-           [a  (fl/ a s)]
-           [c  (fl* a (fl+ 1.0 (flexpt 2.0 27.0)))]
-           [a-hi  (fl- c (fl- c a))]
-           [a-lo  (fl- a a-hi)])
-      (values (fl* a-hi s)
-              (fl* a-lo s)))))
-
-;(: fast-mono-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a+b and its rounding error
-;; Assumes |a| >= |b|
-(define-syntax-rule (fast-mono-fl+/error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let ([x  (+ a b)])
-      (values x (- b (- x a))))))
-
-;(: fast-mono-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a+b and its rounding error
-;; Assumes |a| >= |b|
-(define-syntax-rule (fast-mono-fl-/error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let ([x  (- a b)])
-      (values x (- (- a x) b)))))
-
-;(: fast-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a+b and its rounding error
-(define-syntax-rule (fast-fl+/error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let* ([x  (fl+ a b)]
-           [v  (fl- x a)])
-      (values x (fl+ (fl- a (fl- x v)) (fl- b v))))))
-
-;(: fast-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a-b and its rounding error
-(define-syntax-rule (fast-fl-/error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let* ([x  (fl- a b)]
-           [v  (fl- x a)])
-      (values x (fl- (fl- a (fl- x v)) (fl+ b v))))))
-
-;(: fast-fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a*b and its rounding error
-(define-syntax-rule (fast-fl*/error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let*-values ([(x)  (fl* a b)]
-                  [(a-hi a-lo)  (flsplit a)]
-                  [(b-hi b-lo)  (flsplit b)])
-      (values x (- (fl- (fl- (fl- (fl- x (fl* a-hi b-hi))
-                                  (fl* a-lo b-hi))
-                             (fl* a-hi b-lo))
-                        (fl* a-lo b-lo)))))))
-
-(define-syntax-rule (fast-flsqr/error a-expr)
-  (let ([a a-expr])
-    (let*-values ([(x)  (fl* a a)]
-                  [(a-hi a-lo)  (flsplit a)])
-      (values x (- (fl- (fl- (fl- x (fl* a-hi a-hi))
-                             (fl* 2.0 (fl* a-hi a-lo)))
-                        (fl* a-lo a-lo)))))))
-
-;(: fast-fl//error (Flonum Flonum -> (Values Flonum Flonum)))
-;; Returns a/b and its rounding error
-(define-syntax-rule (fast-fl//error a-expr b-expr)
-  (let ([a a-expr] [b b-expr])
-    (let*-values ([(q-hi)  (fl/ a b)]
-                  [(q0 q1)  (fast-fl*/error q-hi b)])
-      (values q-hi (fl/ (fl+ (- q1) (fl- a q0)) b)))))
-
-#;; If we had a fused multiply-add:
-(define (fast-fl*/error a b)
-  (let ([p  (* a b)])
-    (values p (flfma a b (- p)))))
-
-(define-syntax-rule (define-slow-wrapper name fast-flop/error op)
-  (define-syntax-rule (name a-expr b-expr)
-    (let ([a a-expr] [b b-expr])
-      (let-values ([(x e)  (fast-flop/error a b)])
-        (if (and (x . > . -inf.0) (x . < . +inf.0)
-                 (e . > . -inf.0) (e . < . +inf.0))
-            (values x e)
-            (let* ([v  (op (inexact->exact a) (inexact->exact b))]
-                   [x  (real->double-flonum v)])
-              (if (and (x . > . -inf.0) (x . < . +inf.0))
-                  (values x (real->double-flonum (- v (inexact->exact x))))
-                  (values x 0.0))))))))
-
-(define-slow-wrapper fl+/error fast-fl+/error +)
-(define-slow-wrapper fl-/error fast-fl-/error -)
-(define-slow-wrapper fl*/error fast-fl*/error *)
diff --git a/collects/math/private/flonum/utils.rkt b/collects/math/private/flonum/utils.rkt
new file mode 100644
index 0000000000..f7e374af1d
--- /dev/null
+++ b/collects/math/private/flonum/utils.rkt
@@ -0,0 +1,22 @@
+#lang racket/base
+
+(provide near-pow2
+         near-pow2/div)
+
+(module untyped-defs racket/base
+  (require racket/flonum)
+  
+  (provide (all-defined-out))
+  
+  (define-syntax-rule (near-pow2 a)
+    (flexpt 2.0 (flmax -511.0 (flmin 511.0 (flround (fl/ (fllog (flabs a)) (fllog 2.0)))))))
+  
+  (define-syntax-rule (near-pow2/div a b)
+    ;; Clamping both values makes this work properly when a or b is infinite or zero
+    (let ([ea  (flmax -511.0 (flmin 511.0 (fl/ (fllog (flabs a)) (fllog 2.0))))]
+          [eb  (flmax -511.0 (flmin 511.0 (fl/ (fllog (flabs b)) (fllog 2.0))))])
+      (flexpt 2.0 (flround (fl* 0.5 (fl+ ea eb))))))
+  
+  )  ; module
+
+(require (submod "." untyped-defs))
diff --git a/collects/math/private/utils/flonum-tests.rkt b/collects/math/private/utils/flonum-tests.rkt
new file mode 100644
index 0000000000..3f41f3e355
--- /dev/null
+++ b/collects/math/private/utils/flonum-tests.rkt
@@ -0,0 +1,418 @@
+#lang typed/racket/base
+
+(require racket/math
+         racket/flonum
+         racket/list
+         typed/rackunit
+         "../base/base-random.rkt"
+         "../flonum/expansion/expansion-base.rkt"
+         "../flonum/flonum-functions.rkt"
+         "../flonum/flonum-constants.rkt"
+         "../flonum/flonum-bits.rkt"
+         "../flonum/flonum-error.rkt"
+         "../distributions/dist-struct.rkt"
+         "../distributions/geometric-dist.rkt"
+         "../../bigfloat.rkt")
+
+(provide print-test-progress?
+         test-fpu-arith
+         test-fpu-trig
+         test-fpu-non-trig
+         test-fpu-arith/error
+         test-fpu-arith/fl2
+         test-fpu-non-trig/fl2
+         test-fpu)
+
+;; Allowable error for different kinds of functions, in ulps
+(define flonum-fun-ulps 0.5)
+(define flonum/error-fun-ulps 0.5)
+(define unary-fl2-fun-ulps 1.0)
+(define binary-fl2-fun-ulps 8.0)
+
+;; ===================================================================================================
+;; Helpers
+
+(: different-zero? (Real Real -> Boolean))
+(define (different-zero? x y)
+  (or (and (eqv? x -0.0) (eqv? y 0.0))
+      (and (eqv? x 0.0) (eqv? y -0.0))))
+
+(: fl2->real* (Flonum Flonum -> Real))
+;; Like `fl2->real', but returns signed flonum zeros
+(define (fl2->real* x2 x1)
+  (define x.0 (fl+ x2 x1))
+  (cond [(zero? x.0)  x2]
+        [else  (fl2->real x2 x1)]))
+
+(: bigfloat->real* (Bigfloat -> Real))
+;; Like `bigfloat->real*', but returns a signed infinity or a signed flonum zero if conversion would
+;; overflow or underflow a flonum
+(define (bigfloat->real* x)
+  (define x.0 (bigfloat->flonum x))
+  (cond [(fl= x.0 0.0)  x.0]
+        [(flrational? x.0)  (bigfloat->real x)]
+        [else  x.0]))
+
+(: filter/ulp-error (All (A B) ((Listof (List A (U B Flonum)))
+                                Flonum -> (Listof (List A (U B Flonum))))))
+(define (filter/ulp-error xes ulps)
+  (filter (λ: ([xe : (List A (U B Flonum))])
+            (define e (second xe))
+            (or (not (flonum? e)) (e . fl> . ulps)))
+          xes))
+
+(: print-test-progress? (Parameterof Boolean))
+(define print-test-progress? (make-parameter #t))
+
+(define progress-chunk-size 200)
+(define progress-superchunk-chunks 5)
+
+(: maybe-print-progress (Symbol Integer Natural -> Void))
+(define (maybe-print-progress name i m)
+  (when (and (print-test-progress?) (i . > . 0) (i . <= . m))
+    (let* ([flush?  (cond [(= i 1)  (printf "~a: " name)]
+                          [else  #f])]
+           [flush?  (cond [(= 0 (modulo i progress-chunk-size))
+                           (cond [(= 0 (modulo i (* progress-superchunk-chunks
+                                                    progress-chunk-size)))
+                                  (printf "* ~a " i)]
+                                 [else  (printf "*")])]
+                          [else  flush?])]
+           [flush?  (cond [(= i m)  (printf "* ~a~n" m)]
+                          [else  flush?])])
+      (when flush? (flush-output)))))
+
+;; ===================================================================================================
+;; Test case generation
+
+;; Deteriministic test cases
+
+(define standard-xs
+  (list
+   ;; Test the sign of the return value of `flexpt'
+   -1001.0 -10.0 -0.1 +0.1 +10.0 +1001.0
+   ;; Test squaring
+   (- (flsqrt +min.0)) (flsqrt +min.0)
+   (- (flsqrt +max-subnormal.0)) (flsqrt +max-subnormal.0)
+   (- (flsqrt +max-fl2-subnormal.0)) (flsqrt +max-fl2-subnormal.0)
+   (- (flsqrt +max.0)) (flsqrt +max.0)
+   ;; Test exp limits
+   (fllog +min.0) (fllog +max-subnormal.0) (fllog +max-fl2-subnormal.0) (fllog +max.0)
+   ;; Standard special values
+   -inf.0 -max.0 -1.0 -max-fl2-subnormal.0 -max-subnormal.0 -min.0 -0.0
+   +inf.0 +max.0 +1.0 +max-fl2-subnormal.0 +max-subnormal.0 +min.0 +0.0
+   +nan.0))
+
+(define standard-rs
+  (append standard-xs
+          (list +10 +1 +1/7 +1/10 +1/13
+                -10 -1 -1/7 -1/10 -1/13
+                0)))
+
+(: product (All (A B) ((Listof A) (Listof B) -> (Values (Listof A) (Listof B)))))
+(define (product as bs)
+  (define abs
+    (append*
+     (for/list: : (Listof (Listof (Pair A B))) ([a  (in-list as)])
+       (for/list: : (Listof (Pair A B)) ([b  (in-list bs)])
+         (cons a b)))))
+  (values (map (inst car A B) abs) (map (inst cdr A B) abs)))
+
+;; Random test cases
+
+(define min-subnormal-ord (flonum->ordinal -max-subnormal.0))
+(define max-subnormal-ord (+ 1 (flonum->ordinal +max-subnormal.0)))
+
+(define min-fl2-subnormal-ord (flonum->ordinal -max-fl2-subnormal.0))
+(define max-fl2-subnormal-ord (+ 1 (flonum->ordinal +max-fl2-subnormal.0)))
+
+(: sample-flonum (case-> (Integer -> (Listof Flonum))
+                         (Integer Flonum Flonum -> (Listof Flonum))))
+(define (sample-flonum n [mn -inf.0] [mx +inf.0])
+  (define min-ord (flonum->ordinal mn))
+  (define max-ord (+ 1 (flonum->ordinal mx)))
+  (let ([min-subnormal-ord  (max min-ord min-subnormal-ord)]
+        [max-subnormal-ord  (min max-ord max-subnormal-ord)]
+        [min-fl2-subnormal-ord  (max min-ord min-fl2-subnormal-ord)]
+        [max-fl2-subnormal-ord  (min max-ord max-fl2-subnormal-ord)])
+    (build-list
+     n (λ (_)
+         (define r (random))
+         (ordinal->flonum
+          (cond [(and (min-subnormal-ord . < . max-subnormal-ord) (r . < . 0.1))
+                 (random-integer min-subnormal-ord max-subnormal-ord)]
+                [(and (min-fl2-subnormal-ord . < . max-fl2-subnormal-ord) (r . < . 0.2))
+                 (random-integer min-subnormal-ord max-subnormal-ord)]
+                [else
+                 (random-integer min-ord max-ord)]))))))
+
+(define denom-dist (geometric-dist 1e-32))
+
+(: sample-rational (Integer -> (Listof Exact-Rational)))
+(define (sample-rational n)
+  (map (λ: ([f1 : Flonum] [d : Integer])
+         (+ (inexact->exact f1)
+            (* (if ((random) . > . 0.5) -1 1)
+               (/ (random-natural (+ 1 d)) d)
+               (expt 2 (- (exact-round (/ (fllog (flabs f1)) (fllog 2.0))) 52)))))
+       (sample-flonum n)
+       (map (λ: ([x : Flonum]) (+ 1 (exact-floor x))) (sample denom-dist n))))
+
+;; ===================================================================================================
+;; Flonum functions
+
+(: flonum-error (Flonum Bigfloat -> Any))
+(define (flonum-error z z0.bf)
+  (define z0 (bigfloat->real* z0.bf))
+  (cond [(different-zero? z z0)  (list 'different-zero? z z0)]
+        [else  (flulp-error z z0)]))
+
+(: unary-flonum-fun-error ((Flonum -> Flonum) (Bigfloat -> Bigfloat) Flonum -> Any))
+(define (unary-flonum-fun-error f g x)
+  (flonum-error (f x) (parameterize ([bf-precision 53])
+                        (g (bf x)))))
+
+(: test-unary-flonum-fun
+   (Symbol (Flonum -> Flonum) (Bigfloat -> Bigfloat) Integer Flonum Flonum
+           -> (Listof (List (List Symbol Flonum) Any))))
+(define (test-unary-flonum-fun name f g n mn mx)
+  (define xs (append standard-xs (sample-flonum n mn mx)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum) Any)) ([x  (in-list xs)]
+                                                          [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (list (list name x) (unary-flonum-fun-error f g x)))
+   flonum-fun-ulps))
+
+(: binary-flonum-fun-error
+   ((Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum -> Any))
+(define (binary-flonum-fun-error f g x y)
+  (flonum-error (f x y) (parameterize ([bf-precision 53])
+                          (g (bf x) (bf y)))))
+
+(: test-binary-flonum-fun
+   (Symbol (Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Integer
+           -> (Listof (List (List Symbol Flonum Flonum) Any))))
+(define (test-binary-flonum-fun name f g n)
+  (define-values (pre-xs pre-ys) (product standard-xs standard-xs))
+  (define xs (append pre-xs (sample-flonum n)))
+  (define ys (append pre-ys (sample-flonum n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x  (in-list xs)]
+                                                                 [y  (in-list ys)]
+                                                                 [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (list (list name x y) (binary-flonum-fun-error f g x y)))
+   flonum-fun-ulps))
+
+;; ===================================================================================================
+;; fl2 conversion
+
+(: fl2-error (Flonum Flonum Real -> Any))
+(define (fl2-error x2 x1 x)
+  (cond [(not (fl2? x2 x1))  (list 'not-fl2? x2 x1)]
+        [(different-zero? x2 x)  (list 'different-zero? x2 x)]
+        [else  (fl2ulp-error x2 x1 x)]))
+
+(: fl2-conversion-error (Real -> Any))
+(define (fl2-conversion-error x)
+  (define-values (x2 x1) (fl2 x))
+  (fl2-error x2 x1 x))
+
+(: test-fl2-conversion (Integer -> (Listof (List (List 'fl2 Real) Any))))
+(define (test-fl2-conversion n)
+  (define xs (append standard-rs (sample-rational n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List 'fl2 Real) Any)) ([x  (in-list xs)]
+                                                      [i  (in-naturals 1)])
+     (maybe-print-progress 'fl2 i m)
+     (list (list 'fl2 x) (fl2-conversion-error x)))
+   flonum/error-fun-ulps))
+
+;; ===================================================================================================
+;; Flonum arithmetic with error
+
+(: unary-flonum/error-fun-error ((Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Flonum
+                                                                    -> Any))
+(define (unary-flonum/error-fun-error f g x)
+  (define-values (z2 z1) (f x))
+  (fl2-error z2 z1 (parameterize ([bf-precision 256])
+                     (bigfloat->real* (g (bf x))))))
+
+(: binary-flonum/error-fun-error ((Flonum Flonum -> (Values Flonum Flonum))
+                                  (Bigfloat Bigfloat -> Bigfloat)
+                                  Flonum Flonum
+                                  -> Any))
+(define (binary-flonum/error-fun-error f g x y)
+  (define-values (z2 z1) (f x y))
+  (fl2-error z2 z1 (parameterize ([bf-precision 256])
+                     (bigfloat->real* (g (bf x) (bf y))))))
+
+(: test-unary-flonum/error-fun
+   (Symbol (Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer
+           -> (Listof (List (List Symbol Flonum) Any))))
+(define (test-unary-flonum/error-fun name f g n)
+  (define xs (append standard-xs (sample-flonum n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum) Any)) ([x  (in-list xs)]
+                                                          [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (list (list name x) (unary-flonum/error-fun-error f g x)))
+   flonum/error-fun-ulps))
+
+(: test-binary-flonum/error-fun
+   (Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Integer
+           -> (Listof (List (List Symbol Flonum Flonum) Any))))
+(define (test-binary-flonum/error-fun name f g n)
+  (define-values (pre-xs pre-ys) (product standard-xs standard-xs))
+  (define xs (append pre-xs (sample-flonum n)))
+  (define ys (append pre-ys (sample-flonum n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x  (in-list xs)]
+                                                                 [y  (in-list ys)]
+                                                                 [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (list (list name x y) (binary-flonum/error-fun-error f g x y)))
+   flonum/error-fun-ulps))
+
+;; ===================================================================================================
+;; Flonum expansions
+
+(: unary-fl2-fun-error ((Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat)
+                                                                  Flonum Flonum -> Any))
+(define (unary-fl2-fun-error f g x2 x1)
+  (define-values (z2 z1) (f x2 x1))
+  (fl2-error z2 z1 (parameterize ([bf-precision 256])
+                     (bigfloat->real* (g (bf (fl2->real* x2 x1)))))))
+
+(: test-unary-fl2-fun
+   (Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer
+           -> (Listof (List (List Symbol Flonum Flonum) Any))))
+(define (test-unary-fl2-fun name f g n)
+  (define xs (append standard-rs (sample-rational n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x  (in-list xs)]
+                                                                 [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (define-values (x2 x1) (fl2 x))
+     (list (list name x2 x1) (unary-fl2-fun-error f g x2 x1)))
+   unary-fl2-fun-ulps))
+
+(: binary-fl2-fun-error ((Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))
+                         (Bigfloat Bigfloat -> Bigfloat)
+                         Flonum Flonum Flonum Flonum
+                         -> Any))
+(define (binary-fl2-fun-error f g x2 x1 y2 y1)
+  (define-values (z2 z1) (f x2 x1 y2 y1))
+  (fl2-error z2 z1 (parameterize ([bf-precision 256])
+                     (bigfloat->real* (g (bf (fl2->real* x2 x1)) (bf (fl2->real* y2 y1)))))))
+
+(: test-binary-fl2-fun
+   (Symbol (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat)
+           Integer -> (Listof (List (List Symbol Flonum Flonum Flonum Flonum) Any))))
+(define (test-binary-fl2-fun name f g n)
+  (define-values (pre-xs pre-ys) (product standard-rs standard-rs))
+  (define xs (append pre-xs (sample-rational n)))
+  (define ys (append pre-ys (sample-rational n)))
+  (define m (length xs))
+  (filter/ulp-error
+   (for/list: : (Listof (List (List Symbol Flonum Flonum Flonum Flonum) Any)
+                        ) ([x  (in-list xs)]
+                           [y  (in-list ys)]
+                           [i  (in-naturals 1)])
+     (maybe-print-progress name i m)
+     (define-values (x2 x1) (fl2 x))
+     (define-values (y2 y1) (fl2 y))
+     (list (list name x2 x1 y2 y1) (binary-fl2-fun-error f g x2 x1 y2 y1)))
+   binary-fl2-fun-ulps))
+
+;; ===================================================================================================
+
+(: test-fpu-arith (Natural -> Any))
+(define (test-fpu-arith n)
+  (check-equal? (test-unary-flonum-fun 'flabs flabs bfabs n -inf.0 +inf.0)
+                '())
+  (check-equal? (test-binary-flonum-fun 'fl+ fl+ bf+ n)
+                '())
+  (check-equal? (test-binary-flonum-fun 'fl- fl- bf- n)
+                '())
+  (check-equal? (test-binary-flonum-fun 'fl* fl* bf* n)
+                '())
+  (check-equal? (test-binary-flonum-fun 'fl/ fl/ bf/ n)
+                '()))
+
+(: test-fpu-trig (Natural -> Any))
+(define (test-fpu-trig n)
+  (check-equal? (test-unary-flonum-fun 'flsin flsin bfsin n -inf.0 +inf.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'flcos flcos bfcos n -inf.0 +inf.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'fltan fltan bftan n -inf.0 +inf.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'flasin flasin bfasin n -1.0 1.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'flacos flacos bfacos n -1.0 1.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'flatan flatan bfatan n -inf.0 +inf.0)
+                '()))
+
+(: test-fpu-non-trig (Natural -> Any))
+(define (test-fpu-non-trig n)
+  (check-equal? (test-unary-flonum-fun 'flsqrt flsqrt bfsqrt n 0.0 +inf.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'fllog fllog bflog n 0.0 +inf.0)
+                '())
+  (check-equal? (test-unary-flonum-fun 'flexp flexp bfexp n -746.0 710.0)
+                '())
+  (check-equal? (test-binary-flonum-fun 'flexpt flexpt bfexpt n)
+                '()))
+
+(: test-fpu-arith/error (Natural -> Any))
+(define (test-fpu-arith/error n)
+  (check-equal? (test-binary-flonum/error-fun 'fl+/error fl+/error bf+ n)
+                '())
+  (check-equal? (test-binary-flonum/error-fun 'fl-/error fl-/error bf- n)
+                '())
+  (check-equal? (test-binary-flonum/error-fun 'fl*/error fl*/error bf* n)
+                '())
+  (check-equal? (test-unary-flonum/error-fun 'flsqr/error flsqr/error bfsqr n)
+                '())
+  (check-equal? (test-binary-flonum/error-fun 'fl//error fl//error bf/ n)
+                '()))
+
+(: test-fpu-arith/fl2 (Natural -> Any))
+(define (test-fpu-arith/fl2 n)
+  (check-equal? (test-fl2-conversion n)
+                '())
+  (check-equal? (test-unary-fl2-fun 'fl2abs fl2abs bfabs n)
+                '())
+  (check-equal? (test-binary-fl2-fun 'fl2+ fl2+ bf+ n)
+                '())
+  (check-equal? (test-binary-fl2-fun 'fl2- fl2- bf- n)
+                '())
+  (check-equal? (test-binary-fl2-fun 'fl2* fl2* bf* n)
+                '())
+  (check-equal? (test-unary-fl2-fun 'fl2sqr fl2sqr bfsqr n)
+                '())
+  (check-equal? (test-binary-fl2-fun 'fl2/ fl2/ bf/ n)
+                '()))
+
+(: test-fpu-non-trig/fl2 (Natural -> Any))
+(define (test-fpu-non-trig/fl2 n)
+  (check-equal? (test-unary-fl2-fun 'fl2sqrt fl2sqrt bfsqrt n)
+                '()))
+
+(: test-fpu (Natural -> Any))
+(define (test-fpu n)
+  (test-fpu-arith n)
+  (test-fpu-trig n)
+  (test-fpu-non-trig n)
+  (test-fpu-arith/error n)
+  (test-fpu-arith/fl2 n)
+  (test-fpu-non-trig/fl2 n))
diff --git a/collects/math/tests/flonum-tests.rkt b/collects/math/tests/flonum-tests.rkt
index 0740099935..d3ed5bf3a6 100644
--- a/collects/math/tests/flonum-tests.rkt
+++ b/collects/math/tests/flonum-tests.rkt
@@ -1,20 +1,151 @@
 #lang racket
 
 (require math/flonum
+         math/base
+         math/utils
          rackunit)
 
-(for* ([x  '(+inf.0 +nan.0 -inf.0)]
-       [y  '(+inf.0 +nan.0 -inf.0)])
-  (cond [(eqv? x y)
-         (check-eqv? (flulp-error x y) 0.0 (format "(flulp-error ~v ~v)" x y))]
-        [else
-         (check-eqv? (flulp-error x y) +inf.0 (format "(flulp-error ~v ~v)" x y))]))
+;; ===================================================================================================
+;; Test `flulp-error' heavily; it MUST be correct, or all the FPU tests are suspect
 
-(check-equal? (flulp-error 55123.135123 55123.135123)
-              0.0)
+;; Both arguments eqv?
+(check-equal? (flulp-error -inf.0 -inf.0) 0.0)
+(check-equal? (flulp-error -max.0 -max.0) 0.0)
+(check-equal? (flulp-error -1.0 -1.0) 0.0)
+(check-equal? (flulp-error -min.0 -min.0) 0.0)
+(check-equal? (flulp-error -0.0 -0.0) 0.0)
+(check-equal? (flulp-error +0.0 +0.0) 0.0)
+(check-equal? (flulp-error +min.0 +min.0) 0.0)
+(check-equal? (flulp-error +1.0 +1.0) 0.0)
+(check-equal? (flulp-error +max.0 +max.0) 0.0)
+(check-equal? (flulp-error +inf.0 +inf.0) 0.0)
 
-(check-equal? (flulp-error 1.0 (flnext 1.0))
-              1.0)
+;; Both arguments zero
+(check-equal? (flulp-error -0.0 +0.0) 0.0)
+(check-equal? (flulp-error +0.0 -0.0) 0.0)
 
-(check-equal? (flulp-error +max.0 (flprev +max.0))
-              1.0)
+;; LHS argument +inf.0
+(check-equal? (flulp-error +inf.0 -inf.0) +inf.0)
+(check-equal? (flulp-error +inf.0 -max.0) +inf.0)
+(check-equal? (flulp-error +inf.0 -1.0) +inf.0)
+(check-equal? (flulp-error +inf.0 -min.0) +inf.0)
+(check-equal? (flulp-error +inf.0 -0.0) +inf.0)
+(check-equal? (flulp-error +inf.0 +0.0) +inf.0)
+(check-equal? (flulp-error +inf.0 +min.0) +inf.0)
+(check-equal? (flulp-error +inf.0 +1.0) +inf.0)
+(check-equal? (flulp-error +inf.0 +max.0) +inf.0)
+
+;; LHS argument -inf.0
+(check-equal? (flulp-error -inf.0 -max.0) +inf.0)
+(check-equal? (flulp-error -inf.0 -1.0) +inf.0)
+(check-equal? (flulp-error -inf.0 -min.0) +inf.0)
+(check-equal? (flulp-error -inf.0 -0.0) +inf.0)
+(check-equal? (flulp-error -inf.0 +0.0) +inf.0)
+(check-equal? (flulp-error -inf.0 +min.0) +inf.0)
+(check-equal? (flulp-error -inf.0 +1.0) +inf.0)
+(check-equal? (flulp-error -inf.0 +max.0) +inf.0)
+(check-equal? (flulp-error -inf.0 +inf.0) +inf.0)
+
+;; RHS argument +inf.0
+(check-equal? (flulp-error -max.0 +inf.0) +inf.0)
+(check-equal? (flulp-error -1.0 +inf.0) +inf.0)
+(check-equal? (flulp-error -min.0 +inf.0) +inf.0)
+(check-equal? (flulp-error -0.0 +inf.0) +inf.0)
+(check-equal? (flulp-error +0.0 +inf.0) +inf.0)
+(check-equal? (flulp-error +min.0 +inf.0) +inf.0)
+(check-equal? (flulp-error +1.0 +inf.0) +inf.0)
+(check-equal? (flulp-error +max.0 +inf.0) +inf.0)
+
+;; RHS argument -inf.0
+(check-equal? (flulp-error -max.0 -inf.0) +inf.0)
+(check-equal? (flulp-error -1.0 -inf.0) +inf.0)
+(check-equal? (flulp-error -min.0 -inf.0) +inf.0)
+(check-equal? (flulp-error -0.0 -inf.0) +inf.0)
+(check-equal? (flulp-error +0.0 -inf.0) +inf.0)
+(check-equal? (flulp-error +min.0 -inf.0) +inf.0)
+(check-equal? (flulp-error +1.0 -inf.0) +inf.0)
+(check-equal? (flulp-error +max.0 -inf.0) +inf.0)
+
+;; RHS argument 0.0
+(check-equal? (flulp-error -max.0 +0.0) +inf.0)
+(check-equal? (flulp-error -1.0 +0.0) +inf.0)
+(check-equal? (flulp-error -min.0 +0.0) +inf.0)
+(check-equal? (flulp-error +min.0 +0.0) +inf.0)
+(check-equal? (flulp-error +1.0 +0.0) +inf.0)
+(check-equal? (flulp-error +max.0 +0.0) +inf.0)
+
+;; RHS argument -0.0
+(check-equal? (flulp-error -max.0 -0.0) +inf.0)
+(check-equal? (flulp-error -1.0 -0.0) +inf.0)
+(check-equal? (flulp-error -min.0 -0.0) +inf.0)
+(check-equal? (flulp-error +min.0 -0.0) +inf.0)
+(check-equal? (flulp-error +1.0 -0.0) +inf.0)
+(check-equal? (flulp-error +max.0 -0.0) +inf.0)
+
+;; Small errors
+(check-equal? (flulp-error +0.0 -min.0) 1.0)
+(check-equal? (flulp-error +0.0 +min.0) 1.0)
+(check-equal? (flulp-error -0.0 -min.0) 1.0)
+(check-equal? (flulp-error -0.0 +min.0) 1.0)
+(check-equal? (flulp-error -min.0 +min.0) 2.0)
+(check-equal? (flulp-error +min.0 -min.0) 2.0)
+
+(define large-flulp-error-xys
+  (list (list -1.0 -max.0)
+        (list -1.0 +max.0)
+        (list -min.0 -max.0)
+        (list -min.0 +max.0)
+        (list -1.0 +1.0)
+        (list -max.0 +max.0)
+        (list -max.0 -1.0)
+        (list -max.0 -min.0)
+        (list -max.0 +min.0)
+        (list -max.0 +1.0)
+        (list -1.0 -min.0)
+        (list -1.0 +min.0)
+        (list -min.0 -1.0)
+        (list -min.0 +1.0)
+        (list -0.0 -max.0)
+        (list -0.0 -1.0)
+        (list -0.0 +1.0)
+        (list -0.0 +max.0)
+        (list +0.0 -max.0)
+        (list +0.0 -1.0)
+        (list +0.0 +1.0)
+        (list +0.0 +max.0)
+        (list +min.0 -max.0)
+        (list +min.0 -1.0)
+        (list +min.0 +1.0)
+        (list +min.0 +max.0)
+        (list +1.0 -max.0)
+        (list +1.0 -1.0)
+        (list +1.0 -min.0)
+        (list +1.0 +min.0)
+        (list +1.0 +max.0)
+        (list +max.0 -max.0)
+        (list +max.0 -1.0)
+        (list +max.0 -min.0)
+        (list +max.0 +min.0)
+        (list +max.0 +1.0)))
+
+;; Large errors
+(for ([xy  (in-list large-flulp-error-xys)])
+  (match-define (list x y) xy)
+  (check-true ((flulp-error x y) . >= . (expt 2 52))
+              (format "x = ~a  y = ~a" x y)))
+
+(check-equal? (flulp-error 1.0 (flnext 1.0)) 1.0)
+(check-equal? (flulp-error +max.0 (flprev +max.0)) 1.0)
+
+(for ([_  (in-range 1000)])
+  (define s (random))
+  (define e (fl (random-integer -1074 1024)))
+  (define x (* s (flexpt 2.0 e)))
+  (check-equal? (flulp-error x (flnext x)) 1.0
+                (format "x = ~a" x)))
+
+;; ===================================================================================================
+;; FPU testing
+
+(parameterize ([print-test-progress? #f])
+  (test-fpu 1000))
diff --git a/collects/math/utils.rkt b/collects/math/utils.rkt
index d7e563c155..9581168c70 100644
--- a/collects/math/utils.rkt
+++ b/collects/math/utils.rkt
@@ -1,6 +1,8 @@
 #lang racket
 
-(require "private/parameters.rkt")
+(require "private/parameters.rkt"
+         "private/utils/flonum-tests.rkt")
 
 (provide (all-from-out
-          "private/parameters.rkt"))
+          "private/parameters.rkt"
+          "private/utils/flonum-tests.rkt"))