diff --git a/collects/math/private/bigfloat/bigfloat-mpfr.rkt b/collects/math/private/bigfloat/bigfloat-mpfr.rkt
index 0ab57d4c48..b8e386795b 100644
--- a/collects/math/private/bigfloat/bigfloat-mpfr.rkt
+++ b/collects/math/private/bigfloat/bigfloat-mpfr.rkt
@@ -8,8 +8,6 @@
 
 (require/typed
  "mpfr.rkt"
- ;; Library stuffs
- [mpfr-available?  (-> Boolean)]
  ;; Parameters
  [bf-rounding-mode  (Parameterof Rounding-Mode)]
  [bf-min-precision  Exact-Positive-Integer]
@@ -48,7 +46,7 @@
  [bfsum  ((Listof Bigfloat) -> Bigfloat)]
  [bfmax2  (Bigfloat Bigfloat -> Bigfloat)]
  [bfmin2  (Bigfloat Bigfloat -> Bigfloat)]
- [bfeqv?  (Bigfloat Bigfloat -> Boolean)]
+ [bf=?    (Bigfloat Bigfloat -> Boolean)]
  [bflt?   (Bigfloat Bigfloat -> Boolean)]
  [bflte?  (Bigfloat Bigfloat -> Boolean)]
  [bfgt?   (Bigfloat Bigfloat -> Boolean)]
@@ -126,7 +124,7 @@
     (: bfpred? (Bigfloat Bigfloat * -> Boolean))
     (define (bfpred? x . xs) (fold-binary-pred bfpred2? x xs))))
 
-(define-nary-pred bf=  bfeqv?)
+(define-nary-pred bf=  bf=?)
 (define-nary-pred bf<  bflt?)
 (define-nary-pred bf<= bflte?)
 (define-nary-pred bf>  bfgt?)
@@ -147,8 +145,6 @@
   (bf+ (flonum->bigfloat x1) (flonum->bigfloat x2)))
 
 (provide
- ;; Library stuffs
- mpfr-available?
  ;; Parameters
  bf-rounding-mode
  bf-min-precision
diff --git a/collects/math/private/bigfloat/gmp.rkt b/collects/math/private/bigfloat/gmp.rkt
new file mode 100644
index 0000000000..1bdf823385
--- /dev/null
+++ b/collects/math/private/bigfloat/gmp.rkt
@@ -0,0 +1,166 @@
+#lang racket/base
+
+(require (for-syntax racket/base)
+         ffi/unsafe
+         ffi/unsafe/define
+         racket/runtime-path)
+
+(provide
+ _mp_size_t
+ _limb_t
+ sizeof-limb_t
+ gmp-limb-bits
+ _mpz-pointer
+ new-mpz
+ integer->mpz
+ mpz->integer)
+
+;; ===================================================================================================
+;; Setup
+
+(define-runtime-path libgmp-so 
+  (case (system-type)
+    [(macosx) '(so "libgmp.10.dylib")]
+    [(windows) '(so "libgmp-10.dll")]
+    [else '(so "libgmp")]))
+
+(define gmp-lib (ffi-lib libgmp-so '("10" "3" "") #:fail (λ () #f)))
+
+(define-syntax-rule (get-gmp-fun name type)
+  (get-ffi-obj name gmp-lib type (make-not-available name)))
+
+;; ===================================================================================================
+;; Types
+
+;; Size header, precision, sign, exponent, limb (a "digit" in an _mpz)
+(define _mp_size_t _long)
+(define _limb_t _ulong)
+(define _mp_bitcnt_t _ulong)
+
+(define sizeof-limb_t (ctype-sizeof _limb_t))
+
+;; Number of bits in a limb
+(define gmp-limb-bits (* 8 sizeof-limb_t))
+;; We don't support "nail" builds, which haven't worked since before GMP 4.3 anyway (a "nail" takes
+;; two bits from each limb, and is supposed to speed up carries between limbs on certain systems)
+
+;; GMP's big integers
+(define-cstruct _mpz ([alloc _int] [size _int] [limbs _pointer]))
+
+;; ALWAYS let GMP manage the memory pointed at by the limbs field in an _mpz. It likes to reallocate
+;; limbs, particularly when the _mpz is an output argument.
+
+;; This means that _mpz instances take up memory that Racket can't account for, so they should only
+;; be used as temporary values. Possibly permanent values should be either instances of _mpfr or
+;; Racket's own exact integers.
+
+(define mpz-init (get-gmp-fun '__gmpz_init (_fun _mpz-pointer -> _void)))
+(define mpz-init2 (get-gmp-fun '__gmpz_init2 (_fun _mpz-pointer _mp_bitcnt_t -> _void)))
+(define mpz-clear (get-gmp-fun '__gmpz_clear (_fun _mpz-pointer -> _void)))
+
+;; new-mpz : [_mp_bitcnt_t] -> _mpz-pointer
+;; Creates an _mpz that is managed by the garbage collector, but whose limbs are not, and which is
+;; not relocatable. These are always safe to pass to mpz_* functions.
+(define (new-mpz [bits 0])
+  (define z (make-mpz 0 0 #f))
+  ;; Same as above, but don't let it move:
+  ;(define z (ptr-ref (malloc _mpz 'atomic-interior) _mpz))
+  (if (= bits 0) (mpz-init z) (mpz-init2 z bits))
+  (register-finalizer z (λ (z) (mpz-clear z)))
+  z)
+
+;; ===================================================================================================
+;; Conversions
+
+(define mpz-get-si (get-gmp-fun '__gmpz_get_si (_fun _mpz-pointer -> _long)))
+(define mpz-fits-long? (get-gmp-fun '__gmpz_fits_slong_p (_fun _mpz-pointer -> _int)))
+
+;; integer->mpz : integer -> _mpz
+;; Converts a Racket integer to an _mpz.
+(define (integer->mpz n)
+  (define abs-n (abs n))
+  (define z (new-mpz (integer-length abs-n)))
+  (define len (mpz-alloc z))
+  (define limbs (mpz-limbs z))
+  (let loop ([i 0])
+    (when (i . < . len)
+      (define bit (* i gmp-limb-bits))
+      (ptr-set! limbs _limb_t i (bitwise-bit-field abs-n bit (+ bit gmp-limb-bits)))
+      (loop (+ i 1))))
+  (define size (ceiling (/ (integer-length abs-n) gmp-limb-bits)))
+  (set-mpz-size! z (if (< n 0) (- size) size))
+  z)
+
+;; size+limbs->integer : integer (listof integer) -> integer
+;; Converts a size (which may be negative) and a limb list into an integer.
+(define (size+limbs->integer size limbs)
+  (define len (abs size))
+  (define num
+    (let loop ([i 0] [res  0])
+      (cond [(i . < . len)
+             (define v (ptr-ref limbs _limb_t i))
+             (loop (+ i 1) (bitwise-ior res (arithmetic-shift v (* i gmp-limb-bits))))]
+            [else  res])))
+  (if (negative? size) (- num) num))
+
+;; mpz->integer : _mpz -> integer
+;; Converts an mpz_t to a Racket integer.
+(define (mpz->integer z)
+  (if (zero? (mpz-fits-long? z))
+      (size+limbs->integer (mpz-size z) (mpz-limbs z))
+      (mpz-get-si z)))
+
+;; ===================================================================================================
+;; Number Theoretic Functions
+;; http://gmplib.org/manual/Number-Theoretic-Functions.html#Number-Theoretic-Functions
+
+(define mpz-probab-prime-p 
+  (get-gmp-fun '__gmpz_probab_prime_p  (_fun _mpz-pointer _int -> _int)))
+
+(define (probably-prime n [repetitions 10])
+  (case (mpz-probab-prime-p (integer->mpz n) repetitions)
+    [(0) 'composite]
+    [(1) 'probably-prime]
+    [(2) 'prime]))
+
+(define (prime? n)
+  (case (probably-prime n)
+    [(probably-prime prime) #t]
+    [else #f]))
+
+(define mpz-nextprime
+  (get-gmp-fun '__gmpz_nextprime  (_fun _mpz-pointer _mpz-pointer -> _void)))
+
+(define (next-prime op)
+  (define result (new-mpz))
+  (mpz-nextprime result (integer->mpz op))
+  (mpz->integer result))
+
+(define mpz-invert
+  (get-gmp-fun '__gmpz_invert  (_fun _mpz-pointer _mpz-pointer _mpz-pointer -> _void)))
+
+(define (mod-inverse n modulus)
+  (define result (new-mpz))
+  (mpz-invert result (integer->mpz n) (integer->mpz modulus))
+  (mpz->integer result))
+
+(define mpz-jacobi
+  (get-gmp-fun '__gmpz_jacobi  (_fun _mpz-pointer _mpz-pointer -> _int)))
+
+(define (jacobi a b)
+  (mpz-jacobi (integer->mpz a) (integer->mpz b)))
+
+(define mpz-legendre
+  (get-gmp-fun '__gmpz_legendre  (_fun _mpz-pointer _mpz-pointer -> _int)))
+
+(define (legendre a b)
+  (mpz-legendre (integer->mpz a) (integer->mpz b)))
+
+(define mpz-remove
+  (get-gmp-fun '__gmpz_remove  (_fun _mpz-pointer _mpz-pointer _mpz-pointer -> _int)))
+
+(define (remove-factor z f)
+  (define result (new-mpz))
+  (define n (mpz-remove result (integer->mpz z) (integer->mpz f)))
+  (define z/f (mpz->integer result))
+  (values z/f n))
diff --git a/collects/math/private/bigfloat/mpfr.rkt b/collects/math/private/bigfloat/mpfr.rkt
index a7e158cf95..9ad0c9f51a 100644
--- a/collects/math/private/bigfloat/mpfr.rkt
+++ b/collects/math/private/bigfloat/mpfr.rkt
@@ -1,21 +1,20 @@
 #lang racket/base
 
-(require ffi/unsafe
+(require (for-syntax racket/base)
+         ffi/unsafe
          ffi/unsafe/cvector
          ffi/unsafe/custodian
          ffi/unsafe/define
+         racket/math
          racket/runtime-path
-         racket/list
          racket/promise
          racket/serialize
-         (for-syntax racket/base))
-
-(require (only-in rnrs/arithmetic/bitwise-6
-                  bitwise-first-bit-set))
+         (only-in rnrs/arithmetic/bitwise-6
+                  bitwise-first-bit-set)
+         "gmp.rkt"
+         "utils.rkt")
 
 (provide
- ;; Library stuffs
- mpfr-available?
  ;; Parameters
  bf-rounding-mode
  bf-min-precision
@@ -51,42 +50,25 @@
 ;; ===================================================================================================
 ;; Setup/takedown
 
-;; All MPFR functions and constants are delayed so that libmpfr and libgmp are loaded on first use.
-
-;; This allows the `math' collection to export `math/bigfloat', and for `math/special-functions' to
-;; use MPFR for functions that don't have a Typed Racket implementation yet. On systems without MPFR,
-;; no exceptions will be raised unless a user tries to use those functions.
-
-(define-runtime-path libgmp-so 
-  (case (system-type)
-    [(macosx) '(so "libgmp.10.dylib")]
-    [(windows) '(so "libgmp-10.dll")]
-    [else '(so "libgmp")]))
 (define-runtime-path libmpfr-so
   (case (system-type)
     [(macosx) '(so "libmpfr.4.dylib")]
     [(windows) '(so "libmpfr-4.dll")]
     [else '(so "libmpfr")]))
 
-(define gmp-lib (ffi-lib libgmp-so '("10" "3" "") #:fail (λ () #f)))
 (define mpfr-lib (ffi-lib libmpfr-so '("4" "1" "") #:fail (λ () #f)))
 
-(define-syntax-rule (get-gmp-fun name type)
-  (get-ffi-obj name gmp-lib type (make-not-available name)))
-
 (define-syntax get-mpfr-fun
   (syntax-rules ()
     [(_ name type) (get-mpfr-fun name type (make-not-available name))]
     [(_ name type fail-thunk) (get-ffi-obj name mpfr-lib type fail-thunk)]))
                   
 (define mpfr-free-cache (get-mpfr-fun 'mpfr_free_cache (_fun -> _void)))
-#;; This may be crashing Racket
+
+;; This may be crashing Racket
 (define mpfr-shutdown (register-custodian-shutdown
                        mpfr-free-cache (λ (free) (free))))
 
-(define (mpfr-available?)
-  (and gmp-lib mpfr-lib))
-
 ;; ===================================================================================================
 ;; Parameters: rounding mode, bit precision, printing
 ;; Exponent min and max are not included; they can't be made into parameters, and if we tried they
@@ -97,8 +79,8 @@
 
 ;; minimum precision (1 bit can't be rounded correctly)
 (define bf-min-precision 2)
-;; maximum precision (the number on 64-bit platforms is ridiculously large)
-(define bf-max-precision (- (expt 2 (- (* (ctype-sizeof _long) 8) 1)) 1))
+;; maximum precision (the number when longs are 64 bits is ridiculously large)
+(define bf-max-precision _long-max)
 
 (define bf-precision
   (make-parameter 128 (λ (p) (cond [(p . < . bf-min-precision)  bf-min-precision]
@@ -111,20 +93,9 @@
 ;; Rounding modes (not all of them, just the useful/well-supported ones)
 (define _rnd_t (_enum '(nearest zero up down)))
 
-;; Size header, precision, sign, exponent, limb (part of a bigint)
-(define _mp_size_t _long)
 (define _prec_t _long)
 (define _sign_t _int)
 (define _exp_t _long)
-(define _limb_t _ulong)
-
-;; Also the size of a union of _mp_size_t and _limb_t (which MPFR uses internally)
-(define sizeof-limb_t (ctype-sizeof _limb_t))
-
-;; Number of bits in a limb
-(define gmp-limb-bits (* 8 sizeof-limb_t))
-;; We don't support "nail" builds, which haven't worked since before GMP 4.3 anyway (a "nail" takes
-;; two bits from each limb, and is supposed to speed up carries between limbs on certain systems)
 
 ;; The "limbs" of a bigint are an array of _limb_t, where the element 0 is the length of the array
 ;; in bytes. Entirely reasonable... except that MPFR's bigfloats point at *element 1* for the
@@ -147,9 +118,8 @@
    ))
 
 (define (bigfloat-equal? x1 x2 _)
-  (and (= (bigfloat-sign x1) (bigfloat-sign x2))
-       (or (bfeqv? x1 x2) 
-           (and (bfnan? x1) (bfnan? x2)))))
+  (or (and (bfnan? x1) (bfnan? x2))
+      (bf=? x1 x2)))
 
 (define (canonicalize-sig+exp sig exp)
   (cond [(zero? sig)  (values 0 0)]
@@ -164,16 +134,16 @@
 (define (bfcanonicalize x)
   (cond [(bfzero? x)  (if (zero? (bigfloat-sign x)) (force 0.bf) (force -0.bf))]
         [(bfnan? x)  (force +nan.bf)]
-        [(bfinfinite? x)  (if (bfnegative? x) (force -inf.bf) (force +inf.bf))]
+        [(bfinfinite? x)  (if (zero? (bigfloat-sign x)) (force +inf.bf) (force -inf.bf))]
         [else
          (let*-values ([(sig exp)  (bigfloat-sig+exp x)]
                        [(sig exp)  (canonicalize-sig+exp sig exp)])
            (parameterize ([bf-precision (integer-length sig)])
-             (bf sig exp)))]))
+             (sig+exp->bigfloat sig exp)))]))
 
 (define (bigfloat-hash x recur-hash)
-  (let*-values ([(sig exp)  (bigfloat-sig+exp x)]
-                [(sig exp)  (canonicalize-sig+exp sig exp)])
+  (let*-values ([(x)  (bfcanonicalize x)]
+                [(sig exp)  (bigfloat-sig+exp x)])
     (recur-hash (vector (bigfloat-sign x) sig exp))))
 
 (define bigfloat-deserialize
@@ -193,7 +163,7 @@
      (unless (exact-integer? exp)
        (raise-argument-error 'bigfloat-deserialize "Integer" 2 p sig exp))
      (parameterize ([bf-precision p])
-       (bf sig exp))]))
+       (sig+exp->bigfloat sig exp))]))
 
 (define bigfloat-deserialize-info
   (make-deserialize-info
@@ -216,6 +186,7 @@
 
 ;; mpfr_t: a multi-precision float with rounding (the main data type)
 (define-cstruct _mpfr ([prec _prec_t] [sign _sign_t] [exp _exp_t] [d (_gcable _mpfr_limbs)])
+  #:property prop:custom-print-quotable 'never
   #:property prop:custom-write (λ (b port mode) (bigfloat-custom-write b port mode))
   #:property prop:equal+hash (list bigfloat-equal? bigfloat-hash bigfloat-hash)
   #:property prop:serializable bigfloat-serialize-info)
@@ -231,39 +202,25 @@
 (define (new-mpfr prec)
   (define n (add1 (quotient (- prec 1) gmp-limb-bits)))
   (define size (* sizeof-limb_t (+ n 1)))
-  ;; Allocate d so that it won't be traced (atomic) or moved (interior)
+  ;; Allocate d and x so they won't be traced (atomic) or moved (interior)
   (define d (make-cvector* (malloc size 'atomic-interior) _limb_t (+ n 1)))
   (define x (make-mpfr prec 1 0 d))
-  ;; Use a finalizer to keep a reference to d as long as x is alive (x's memory isn't traced because
-  ;; it's allocated using make-mpfr; this is equivalent to tracing through d only)
+  ;; The next five lines are like the one above, but malloc uses 'atomic-interior instead of 'atomic
+  ;(define x (ptr-ref (malloc (ctype-sizeof _mpfr) 'atomic-interior) _mpfr))
+  ;(set-mpfr-prec! x prec)
+  ;(set-mpfr-sign! x 1)
+  ;(set-mpfr-exp! x 0)
+  ;(set-mpfr-d! x d)
+  ;; Use a finalizer to keep a reference to d as long as x is alive (this is equivalent to tracing
+  ;; through d)
   (register-finalizer x (λ (x) d))
   (mpfr-set-nan x)
   x)
 
-;; We always create mpfr_ts using new-mpfr. In doing so, we assume that no mpfr_* function will ever
-;; try to reallocate limbs. This is a good assumption because an mpfr_t's precision is fixed from
-;; when it's allocated to when it's deallocated. (There's no reason to allocate new limbs for an
-;; mpfr_t without changing its precision.)
-
-;; Big integers, big rationals
-(define-cstruct _mpz ([alloc _int] [size _int] [limbs _pointer]))
-(define-cstruct _mpq ([num _mpz] [den _mpz]))
-
-;; BE CAREFUL WITH THESE. If you make one with make-mpz or make-mpq, DO NOT send it to a function
-;; that will reallocate its limbs. In particular, NEVER use it as an output argument. However, you
-;; can generally use it as an input argument.
-
-;; MPFR memory management for mpz_t
-(define mpz-init (get-gmp-fun '__gmpz_init (_fun _mpz-pointer -> _void)))
-(define mpz-clear (get-gmp-fun '__gmpz_clear (_fun _mpz-pointer -> _void)))
-
-;; raw-mpz : -> _mpz-pointer
-;; Creates an mpz_t that is managed by the garbage collector, but whose limbs are not. These are
-;; always safe to pass to mpz_* functions. We use them for output parameters.
-(define (raw-mpz)
-  (define x (ptr-ref (malloc _mpz 'atomic-interior) _mpz))
-  (mpz-init x)
-  x)
+;; We always create _mpfr instances using new-mpfr. In doing so, we assume that no mpfr_* function
+;; will ever try to reallocate limbs. This is a good assumption because an _mpfr's precision is
+;; fixed from when it's allocated to when it's deallocated. (There's no reason to allocate new limbs
+;; for an _mpfr without changing its precision.)
 
 ;; ===================================================================================================
 ;; Accessors
@@ -293,10 +250,9 @@
 ;; bigfloat-sig+exp : bigfloat -> integer integer
 ;; Returns the signed significand and exponent of a bigfloat.
 (define (bigfloat-sig+exp x)
-  (define z (raw-mpz))
+  (define z (new-mpz))
   (define exp (mpfr-get-z-2exp z x))
   (define sig (mpz->integer z))
-  (mpz-clear z)
   (values sig exp))
 
 ;; bigfloat-significand : bigfloat -> integer
@@ -311,35 +267,14 @@
 (define mpfr-set-d  (get-mpfr-fun 'mpfr_set_d  (_fun _mpfr-pointer _double _rnd_t -> _void)))
 (define mpfr-set-si (get-mpfr-fun 'mpfr_set_si (_fun _mpfr-pointer _long _rnd_t -> _void)))
 (define mpfr-set-z  (get-mpfr-fun 'mpfr_set_z  (_fun _mpfr-pointer _mpz-pointer _rnd_t -> _void)))
-(define mpfr-set-q  (get-mpfr-fun 'mpfr_set_q  (_fun _mpfr-pointer _mpq-pointer _rnd_t -> _void)))
+(define mpfr-set-z-2exp
+  (get-mpfr-fun 'mpfr_set_z_2exp (_fun _mpfr-pointer _mpz-pointer _exp_t _rnd_t -> _int)))
 
-;; integer->size+limbs : integer -> (values integer (listof integer))
-;; Returns a cvector of limbs and the size of the limbs. The size is negated when n is negative.
-(define (integer->size+limbs n)
-  ;; +1 because GMP expects the last limb to be 0
-  (define len (+ (ceiling (/ (integer-length (abs n)) gmp-limb-bits)) 1))
-  (define limbs (make-cvector _limb_t len))
-  (define an (abs n))
-  (let loop ([i 0])
-    (when (i . < . len)
-      (define bit (* i gmp-limb-bits))
-      (cvector-set! limbs i (bitwise-bit-field an bit (+ bit gmp-limb-bits)))
-      (loop (+ i 1))))
-  (define size (- len 1))
-  (values (if (< n 0) (- size) size)
-          (cvector-ptr limbs)))
-
-;; integer->mpz : integer -> _mpz
-;; Converts an integer to an _mpz. DO NOT send the result of this as an output argument!
-(define (integer->mpz n)
-  (let-values ([(size limbs)  (integer->size+limbs n)])
-    (make-mpz (abs size) size limbs)))
-
-;; rational->mpq : rational -> _mpz
-;; Converts a rational to an _mpq. DO NOT send the result of this as an output argument!
-(define (rational->mpq r)
-  (make-mpq (integer->mpz (numerator r))
-            (integer->mpz (denominator r))))
+;; sig+exp->bigfloat integer integer -> bigfloat
+(define (sig+exp->bigfloat n e)
+  (define y (new-mpfr (bf-precision)))
+  (mpfr-set-z-2exp y (integer->mpz n) e (bf-rounding-mode))
+  y)
 
 ;; flonum->bigfloat : float -> bigfloat
 ;; Converts a Racket inexact real to a bigfloat; rounds if bf-precision < 53.
@@ -352,51 +287,66 @@
 ;; Converts a Racket integer to a bigfloat; rounds if necessary.
 (define (integer->bigfloat value)
   (define x (new-mpfr (bf-precision)))
-  (if (fixnum? value)
+  (if (_long? value)
       (mpfr-set-si x value (bf-rounding-mode))
       (mpfr-set-z x (integer->mpz value) (bf-rounding-mode)))
   x)
 
+(define (round/mode q)
+  (case (bf-rounding-mode)
+    [(up)    (ceiling q)]
+    [(down)  (floor q)]
+    [(zero)  (truncate q)]
+    [else    (round q)]))
+
+(define (floor-log2 n) (max 0 (sub1 (integer-length n))))
+(define (ceiling-log2 n) (max 0 (integer-length (sub1 n))))
+
+(define (log2-lower-bound q)
+  (- (floor-log2 (numerator q))
+     (ceiling-log2 (denominator q))))
+
+(define (log2-upper-bound q)
+  (- (ceiling-log2 (numerator q))
+     (floor-log2 (denominator q))))
+
 ;; rational->bigfloat : rational -> bigfloat
 ;; Converts a Racket rational to a bigfloat; rounds if necessary.
-(define (rational->bigfloat value)
-  (define x (new-mpfr (bf-precision)))
-  (mpfr-set-q x (rational->mpq value) (bf-rounding-mode))
-  x)
+(define (rational->bigfloat q)
+  (define bits (bf-precision))
+  (define sgn-q (sgn q))
+  (define abs-q (abs q))
+  (define ipart (floor abs-q))
+  (cond [(zero? sgn-q)  (sig+exp->bigfloat 0 0)]
+        [(zero? ipart)
+         (define e-ub (- (log2-upper-bound abs-q) bits))
+         (define e-lb (- (log2-lower-bound abs-q) bits))
+         (let loop ([e  e-ub])
+           (define 2^-e (arithmetic-shift 1 (- e)))
+           (define sig (round/mode (* sgn-q abs-q 2^-e)))
+           (cond [(or (= e e-lb) (= bits (integer-length (abs sig))))
+                  (sig+exp->bigfloat sig e)]
+                 [else
+                  (loop (- e 1))]))]
+        [else
+         (define e (- (integer-length ipart) bits))
+         (define 2^-e (cond [(e . < . 0)  (arithmetic-shift 1 (- e))]
+                            [else  (/ 1 (arithmetic-shift 1 e))]))
+         (define sig (round/mode (* sgn-q abs-q 2^-e)))
+         (sig+exp->bigfloat sig e)]))
 
 ;; real->bigfloat : real -> bigfloat
 ;; Converts any real Racket value to a bigfloat; rounds if necessary.
 (define (real->bigfloat value)
   (cond [(inexact? value)  (flonum->bigfloat value)]
         [(integer? value)  (integer->bigfloat value)]
-        [(rational? value)  (rational->bigfloat value)]))
+        [else  (rational->bigfloat value)]))
 
 ;; ===================================================================================================
 ;; Conversion from mpfr_t to Racket data types
 
 (define mpfr-get-d (get-mpfr-fun 'mpfr_get_d (_fun _mpfr-pointer _rnd_t -> _double)))
 (define mpfr-get-z (get-mpfr-fun 'mpfr_get_z (_fun _mpz-pointer _mpfr-pointer _rnd_t -> _int)))
-(define mpz-get-si (get-mpfr-fun '__gmpz_get_si (_fun _mpz-pointer -> _long)))
-(define mpz-fits-long? (get-gmp-fun '__gmpz_fits_slong_p (_fun _mpz-pointer -> _int)))
-
-;; size+limbs->integer : integer (listof integer) -> integer
-;; Converts a size (which may be negative) and a limb list into an integer.
-(define (size+limbs->integer size limbs)
-  (define len (abs size))
-  (define num
-    (let loop ([i 0] [res  0])
-      (cond [(i . < . len)
-             (define v (ptr-ref limbs _limb_t i))
-             (loop (+ i 1) (bitwise-ior res (arithmetic-shift v (* i gmp-limb-bits))))]
-            [else  res])))
-  (if (negative? size) (- num) num))
-
-;; mpz->integer : _mpz -> integer
-;; Converts an mpz_t to an integer.
-(define (mpz->integer z)
-  (if (zero? (mpz-fits-long? z))
-      (size+limbs->integer (mpz-size z) (mpz-limbs z))
-      (mpz-get-si z)))
 
 ;; bigfloat->flonum : bigfloat -> float
 ;; Converts a bigfloat to a Racket float; rounds if necessary.
@@ -407,10 +357,9 @@
 ;; Converts a bigfloat to a Racket integer; rounds if necessary.
 (define (bigfloat->integer x)
   (unless (bfinteger? x) (raise-argument-error 'bigfloat->integer "bfinteger?" x))
-  (define z (raw-mpz))
+  (define z (new-mpz))
   (mpfr-get-z z x (bf-rounding-mode))
   (define res (mpz->integer z))
-  (mpz-clear z)
   res)
 
 ;; bigfloat->rational : bigfloat -> rational
@@ -522,33 +471,23 @@
 
 (define (bigfloat-custom-write x port mode)
   (write-string
-   (if (mpfr-available?)
-       (cond [(bfzero? x)  (if (= 0 (bigfloat-sign x)) "0.bf" "-0.bf")]
-             [(bfrational? x)
-              (define str (bigfloat->string x))
-              (cond [(regexp-match #rx"\\.|e" str)
-                     (define exp (bigfloat-exponent x))
-                     (define prec (bigfloat-precision x))
-                     (if ((abs exp) . > . (* prec 2))
-                         (format "(bf \"~a\")" str)
-                         (format "(bf #e~a)" str))]
-                    [else  (format "(bf ~a)" str)])]
-             [(bfinfinite? x)  (if (= 0 (bigfloat-sign x)) "+inf.bf" "-inf.bf")]
-             [else  "+nan.bf"])
-       "#<bigfloat>")
+   (cond [(bfzero? x)  (if (= 0 (bigfloat-sign x)) "0.bf" "-0.bf")]
+         [(bfrational? x)
+          (define str (bigfloat->string x))
+          (cond [(regexp-match #rx"\\.|e" str)
+                 (define exp (bigfloat-exponent x))
+                 (define prec (bigfloat-precision x))
+                 (if ((abs exp) . > . (* prec 2))
+                     (format "(bf \"~a\")" str)
+                     (format "(bf #e~a)" str))]
+                [else  (format "(bf ~a)" str)])]
+         [(bfinfinite? x)  (if (= 0 (bigfloat-sign x)) "+inf.bf" "-inf.bf")]
+         [else  "+nan.bf"])
    port))
 
 ;; ===================================================================================================
 ;; Main bigfloat constructor
 
-(define mpfr-set-z-2exp
-  (get-mpfr-fun 'mpfr_set_z_2exp (_fun _mpfr-pointer _mpz-pointer _exp_t _rnd_t -> _int)))
-
-(define (sig+exp->bigfloat n e)
-  (define y (new-mpfr (bf-precision)))
-  (mpfr-set-z-2exp y (integer->mpz n) e (bf-rounding-mode))
-  y)
-
 ;; bf : (or real string) -> bigfloat
 ;;    : integer integer -> bigfloat
 (define bf
@@ -804,13 +743,13 @@
 (define mpfr-yn (get-mpfr-fun 'mpfr_yn (_fun _mpfr-pointer _long _mpfr-pointer _rnd_t -> _int)))
 
 (define (bfbesj n x)
-  (unless (fixnum? n) (raise-argument-error 'bfbesj "Fixnum" 0 n x))
+  (unless (_long? n) (raise-argument-error 'bfbesj "_long?" 0 n x))
   (define y (new-mpfr (bf-precision)))
   (mpfr-jn y n x (bf-rounding-mode))
   y)
 
 (define (bfbesy n x)
-  (unless (fixnum? n) (raise-argument-error 'bfbesy "Fixnum" 0 n x))
+  (unless (_long? n) (raise-argument-error 'bfbesy "_long?" 0 n x))
   (define y (new-mpfr (bf-precision)))
   (mpfr-yn y n x (bf-rounding-mode))
   y)
@@ -818,7 +757,7 @@
 (define mpfr-root (get-mpfr-fun 'mpfr_root (_fun _mpfr-pointer _mpfr-pointer _ulong _rnd_t -> _int)))
 
 (define (bfroot x n)
-  (unless (and (fixnum? n) (n . >= . 0)) (raise-argument-error 'bfroot "Nonnegative-Fixnum" 1 x n))
+  (unless (and (_ulong? n) (n . >= . 0)) (raise-argument-error 'bfroot "_ulong?" 1 x n))
   (define y (new-mpfr (bf-precision)))
   (mpfr-root y x n (bf-rounding-mode))
   y)
@@ -908,7 +847,7 @@
   (begin (provide-2ary-pred name c-name) ...))
 
 (provide-2ary-preds
- [bfeqv? 'mpfr_equal_p]
+ [bf=? 'mpfr_equal_p]
  [bflt?  'mpfr_less_p]
  [bflte? 'mpfr_lessequal_p]
  [bfgt?  'mpfr_greater_p]
@@ -989,7 +928,7 @@
 
 (define (epsilon.bf)
   (define p (bf-precision))
-  (hash-ref! constant-hash (cons 'epsilon.bf p) (λ () (bfexpt (force 2.bf) (bf (- p))))))
+  (hash-ref! constant-hash (cons 'epsilon.bf p) (λ () (bfexpt (force 2.bf) (bf (- 1 p))))))
 
 (define (-max.bf) (bfnext (bf -inf.0)))
 (define (-min.bf) (bfprev (bf -0.0)))
@@ -1000,64 +939,3 @@
 (begin-for-syntax
   (set! 0ary-funs (list* #'phi.bf #'epsilon.bf #'-max.bf #'-min.bf #'+min.bf #'+max.bf
                          0ary-funs)))
-
-;; ===================================================================================================
-;; Number Theoretic Functions
-;; http://gmplib.org/manual/Number-Theoretic-Functions.html#Number-Theoretic-Functions
-
-(define mpz-probab-prime-p 
-  (get-gmp-fun '__gmpz_probab_prime_p  (_fun _mpz-pointer _int -> _int)))
-
-(define (probably-prime n [repetitions 10])
-  (case (mpz-probab-prime-p (integer->mpz n) repetitions)
-    [(0) 'composite]
-    [(1) 'probably-prime]
-    [(2) 'prime]))
-
-(define (prime? n)
-  (case (probably-prime n)
-    [(probably-prime prime) #t]
-    [else #f]))
-
-(define mpz-nextprime
-  (get-gmp-fun '__gmpz_nextprime  (_fun _mpz-pointer _mpz-pointer -> _void)))
-
-(define (next-prime op)
-  (define result (raw-mpz))
-  (mpz-nextprime result (integer->mpz op))
-  (begin0
-    (mpz->integer result)
-    (mpz-clear result)))
-
-(define mpz-invert
-  (get-gmp-fun '__gmpz_invert  (_fun _mpz-pointer _mpz-pointer _mpz-pointer -> _void)))
-
-(define (mod-inverse n modulus)
-  (define result (raw-mpz))
-  (mpz-invert result (integer->mpz n) (integer->mpz modulus))
-  (begin0
-    (mpz->integer result)
-    (mpz-clear result)))
-
-(define mpz-jacobi
-  (get-gmp-fun '__gmpz_jacobi  (_fun _mpz-pointer _mpz-pointer -> _int)))
-
-(define (jacobi a b)
-  (mpz-jacobi (integer->mpz a) (integer->mpz b)))
-
-(define mpz-legendre
-  (get-gmp-fun '__gmpz_legendre  (_fun _mpz-pointer _mpz-pointer -> _int)))
-
-(define (legendre a b)
-  (mpz-legendre (integer->mpz a) (integer->mpz b)))
-
-(define mpz-remove
-  (get-gmp-fun '__gmpz_remove  (_fun _mpz-pointer _mpz-pointer _mpz-pointer -> _int)))
-
-(define (remove-factor z f)
-  (define result (raw-mpz))
-  (define n (mpz-remove result (integer->mpz z) (integer->mpz f)))
-  (define z/f (mpz->integer result))
-  (mpz-clear result)
-  (values z/f n))
-
diff --git a/collects/math/private/bigfloat/utils.rkt b/collects/math/private/bigfloat/utils.rkt
index d9dcb2961b..88d043ffa7 100644
--- a/collects/math/private/bigfloat/utils.rkt
+++ b/collects/math/private/bigfloat/utils.rkt
@@ -1,8 +1,24 @@
-#lang typed/racket/base
+#lang racket/base
 
-(require (for-syntax racket/base racket/syntax syntax/strip-context))
+(require (for-syntax racket/base racket/syntax syntax/strip-context)
+         typed/racket/base
+         (only-in ffi/unsafe
+                  ctype-sizeof
+                  _long
+                  _ulong))
 
-(provide req/prov-uniform-collection)
+(provide (all-defined-out))
+
+(define (unsigned-max type) (- (expt 2 (* 8 (ctype-sizeof type))) 1))
+(define (signed-min type) (- (expt 2 (- (* 8 (ctype-sizeof type)) 1))))
+(define (signed-max type) (- (expt 2 (- (* 8 (ctype-sizeof type)) 1)) 1))
+
+(define _long-min (signed-min _long))
+(define _long-max (signed-max _long))
+(define _ulong-max (unsigned-max _ulong))
+
+(define (_ulong? n) (and (exact-integer? n) (<= 0 n _ulong-max)))
+(define (_long? n) (and (exact-integer? n) (<= _long-min n _long-max)))
 
 (define-syntax (req/prov-uniform-collection stx)
   (syntax-case stx ()
diff --git a/collects/math/scribblings/math-bigfloat.scrbl b/collects/math/scribblings/math-bigfloat.scrbl
index 79b27ee7c5..86b0fed1eb 100644
--- a/collects/math/scribblings/math-bigfloat.scrbl
+++ b/collects/math/scribblings/math-bigfloat.scrbl
@@ -28,8 +28,7 @@ represented by the type @racket[Bigfloat].
 
 MPFR is free and license-compatible with commercial software. It is distributed with Racket
 for Windows and Mac OS X, is installed on most Linux systems, and is
-@hyperlink["http://www.mpfr.org/ports.html"]{easy to install} on all major platforms.
-See @secref{loading} for details.
+@hyperlink["http://www.mpfr.org/ports.html"]{easy to install} on major Unix-like platforms.
 
 @local-table-of-contents[]
 
@@ -709,21 +708,4 @@ Canonicalizing bigfloats won't change answers computed from them.
                             (eval:result @racketresultfont{2}))]
 }
 
-@section[#:tag "loading"]{Loading the MPFR Library}
-
-The @racketmodname[math/bigfloat] library depends directly on external one library:
-MPFR, which on the filesystem is named something like @filepath{libmpfr.so}, but with a
-platform-specific suffix. MPFR itself depends on @hyperlink["http://gmplib.org/"]{GMP}, the
-GNU Multiple Precision Arithmetic Library, which on the filesystem is named something like
-@filepath{libgmp.so}.
-
-@defproc[(mpfr-available?) Boolean]{
-If @racket[(mpfr-available?)] is @racket[#t], both MPFR and GMP have been loaded
-by the Racket runtime, and everything in @racketmodname[math/bigfloat] should work.
-
-If @racket[(mpfr-available?)] is @racket[#t] and something doesn't work, it's likely an error
-in the foreign function interface. In that case, please email the author of this library or
-submit a bug report.
-}
-
 @(close-eval untyped-eval)
diff --git a/collects/math/tests/bigfloat-stress-test.rkt b/collects/math/tests/bigfloat-stress-test.rkt
new file mode 100644
index 0000000000..701f7f4f12
--- /dev/null
+++ b/collects/math/tests/bigfloat-stress-test.rkt
@@ -0,0 +1,35 @@
+#lang racket
+
+(require math/bigfloat
+         math/flonum
+         math/special-functions
+         math/private/bigfloat/bigfloat-hurwitz-zeta)
+
+(collect-garbage)
+(collect-garbage)
+(collect-garbage)
+(define start-mem (current-memory-use))
+
+(define num 60)
+
+(for ([x  (in-range 1 num)])
+  (let ([x  (* 0.25 x)])
+    (printf "x = ~v~n" x)
+    (for ([y  (in-range 1 num)])
+      (let ([y  (* 0.25 y)])
+        (define z (bigfloat->flonum (bfhurwitz-zeta (bf x) (bf y))))
+        (define err (flulp-error (flhurwitz-zeta x y) z))
+        (unless (err . <= . 5.0)
+          (printf "hurwitz-zeta ~v ~v; error = ~v ulps~n" x y err))))))
+
+(collect-garbage)
+(collect-garbage)
+(collect-garbage)
+(define end-mem (current-memory-use))
+
+(define diff-mem (- end-mem start-mem))
+(when (diff-mem . > . (* 2 1024 1024))
+  (printf "Significant higher memory use (more than 2MB):~n")
+  (printf "   start-mem = ~v~n" start-mem)
+  (printf "     end-mem = ~v~n" end-mem)
+  (printf "  difference = ~v~n" diff-mem))
diff --git a/collects/math/tests/bigfloat-tests.rkt b/collects/math/tests/bigfloat-tests.rkt
index 58fd5cb230..43bfa5ac7f 100644
--- a/collects/math/tests/bigfloat-tests.rkt
+++ b/collects/math/tests/bigfloat-tests.rkt
@@ -1,11 +1,11 @@
-#lang typed/racket
+#lang racket
 
 (require math/bigfloat
          math/flonum
-         typed/rackunit)
+         math/base
+         rackunit)
 
-;; ---------------------------------------------------------------------------------------------------
-;; Exact
+;; Exact tests
 
 (check-= (bigfloat->integer (bf (expt 2 64)))
          (expt 2 64)  0)
@@ -51,12 +51,11 @@
      "987931636889230098793127736178215424999229576351482208269895"
      "1936680331825288693984964651058209392398294887933203625094431"))
   
-  ;; truncates to 600 digits (many other numbers of digits, correct rounding makes it look wrong)
-  (check-equal? (substring (bigfloat->string (bfexp (bf 1))) 0 600)
+  ;; truncates to 600 digits (with other numbers of digits, correct rounding makes it look wrong)
+  (check-equal? (substring (bigfloat->string (bfexp 1.bf)) 0 600)
                 e-truncated600))
 
-;; ---------------------------------------------------------------------------------------------------
-;; Inexact
+;; Inexact tests
 
 (check-= (bigfloat->rational (bfsqr (bf 22/7)))
          (/ (sqr 22) (sqr 7))
@@ -73,3 +72,401 @@
 (check-= (bigfloat->flonum (bfcbrt (bf 2)))
          (flexpt 2.0 #i1/3)
          epsilon.0)
+
+;; Constants
+
+(check-equal? pi.bf (bf #e3.141592653589793238462643383279502884195))
+(check-equal? gamma.bf (bf #e0.5772156649015328606065120900824024310432))
+(check-equal? catalan.bf (bf #e0.9159655941772190150546035149323841107734))
+(check-equal? log2.bf (bf #e0.6931471805599453094172321214581765680748))
+(check-equal? phi.bf (bf #e1.61803398874989484820458683436563811772))
+
+;; Special values
+
+(check-eqv? (bigfloat->flonum +nan.bf) +nan.0)
+(check-eqv? (bigfloat->flonum +inf.bf) +inf.0)
+(check-eqv? (bigfloat->flonum -inf.bf) -inf.0)
+(check-eqv? (bigfloat->flonum -0.bf) -0.0)
+(check-eqv? (bigfloat->flonum 0.bf) 0.0)
+(check-equal? (bf- +inf.bf) -inf.bf)
+(check-equal? (bf- -inf.bf) +inf.bf)
+(check-equal? (bf- +nan.bf) +nan.bf)
+(check-equal? (bf- 0.bf) -0.bf)
+(check-equal? (bf- -0.bf) 0.bf)
+
+(for ([bits  '(2 4 8 16 32 64 128 3 5 7 11 13 129)])
+  (parameterize ([bf-precision bits])
+    ;; Epsilon
+    (check-not-equal? (bf+ 1.bf epsilon.bf) 1.bf)
+    (check-equal? (bf+ 1.bf (bf* (bf 0.5) epsilon.bf)) 1.bf)
+    ;; +max.bf/-max.bf
+    (check-equal? (bf- +max.bf) -max.bf)
+    (check-equal? (bf+ +max.bf (bf* epsilon.bf +max.bf)) +inf.bf)
+    (check-equal? (bf- -max.bf (bf* epsilon.bf +max.bf)) -inf.bf)
+    (check-equal? (- (bigfloat->ordinal +inf.bf) 1) (bigfloat->ordinal +max.bf))
+    (check-equal? (+ (bigfloat->ordinal -inf.bf) 1) (bigfloat->ordinal -max.bf))
+    (check-equal? (bfnext +max.bf) +inf.bf)
+    (check-equal? (bfprev +inf.bf) +max.bf)
+    (check-equal? (bfprev -max.bf) -inf.bf)
+    (check-equal? (bfnext -inf.bf) -max.bf)
+    ;; +min.bf/-min.bf
+    (check-equal? (bf- +min.bf) -min.bf)
+    (check-equal? (bf* (bf 0.5) +min.bf)  0.bf)
+    (check-equal? (bf* (bf 0.5) -min.bf) -0.bf)
+    (check-equal? (bigfloat->ordinal +min.bf)  1)
+    (check-equal? (bigfloat->ordinal -min.bf) -1)
+    (check-equal? (bfprev +min.bf) 0.bf)
+    (check-equal? (bfnext 0.bf) +min.bf)
+    (check-equal? (bfnext -min.bf) -0.bf)
+    (check-equal? (bfprev 0.bf) -min.bf)
+    ))
+
+;; Integer conversion
+
+(check-equal? (bigfloat->rational (integer->bigfloat 0)) 0)
+
+(for ([mode  (in-list '(nearest up down zero))])
+  (define eps (if (eq? mode 'nearest) (bf* (bf 0.5) epsilon.bf) epsilon.bf))
+  (parameterize ([bf-rounding-mode mode])
+    (for* ([s  (in-list '(-1 1))]
+           [i  (in-range 1000)])
+      (define n (* s (expt 10 i)))
+      (define n0 (bigfloat->integer (integer->bigfloat n)))
+      (define err (bf (relative-error n0 n)))
+      (unless (err . bf<= . eps)
+        (printf "bad integer conversion: mode = ~v  n = ~v~n" mode n))
+      (check-true (err . bf<= . eps)))))
+
+;; Rational conversion
+
+(for ([mode  (in-list '(nearest up down zero))])
+  (parameterize ([bf-rounding-mode mode])
+    (for* ([i  (in-range -37 37)]
+           [j  (in-range 1 37)])
+      (define q0 (bigfloat->rational (rational->bigfloat (/ i j))))
+      (define q1 (bigfloat->rational (bf/ (bf i) (bf j))))
+      (unless (equal? q0 q1)
+        (printf "bad rational conversion: mode = ~v  q = ~v/~v~n" mode i j))
+      (check-equal? q0 q1))))
+
+;; Flonum conversion
+
+(for ([mode  (in-list '(nearest up down zero))])
+  (parameterize ([bf-rounding-mode mode])
+    (for* ([i  (in-range -37 37)]
+           [j  (in-range 1 37)])
+      (define x (fl (/ i j)))
+      (define x0 (bigfloat->flonum (flonum->bigfloat x)))
+      (unless (equal? x x0)
+        (printf "bad flonum conversion: mode = ~v  x = ~v~n" mode x))
+      (check-equal? x x0))))
+
+;; String conversion
+
+(for ([mode  (in-list '(nearest up down zero))])
+  (parameterize ([bf-rounding-mode mode])
+    (for* ([i  (in-range -37 37)]
+           [j  (in-range 1 37)])
+      (define q (rational->bigfloat (/ i j)))
+      (define q0 (string->bigfloat (bigfloat->string q)))
+      (unless (equal? q q0)
+        (printf "bad string conversion: mode = ~v  q = ~v/~v~n" mode i j))
+      (check-equal? q q0))))
+
+(parameterize ([bf-precision 53])
+  (for ([f+xs+ys
+         (list
+          (list
+           bfsqr
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +inf.bf +inf.bf 1.bf 0.bf 0.bf 0.bf 0.bf 1.bf +inf.bf +inf.bf +nan.bf))
+          (list
+           bfsqrt
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +max.bf
+                 +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -0.bf 0.bf
+                 (bf "4.8811524304081624e-161614249") 1.bf (bf "1.4486472022088012e161614248")
+                 +inf.bf +nan.bf))
+          (list
+           bfcbrt
+           (list -inf.bf
+                 -max.bf -1.bf -min.bf
+                 -0.bf 0.bf
+                 +min.bf 1.bf +max.bf
+                 +inf.bf +nan.bf)
+           (list -inf.bf
+                 (bf "-1.2802902004094324e107742832") -1.bf (bf "-6.1993798416193228e-107742833")
+                 -0.bf 0.bf
+                 (bf "6.1993798416193228e-107742833") 1.bf (bf "1.2802902004094324e107742832")
+                 +inf.bf +nan.bf))
+          (list
+           bf-
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +inf.bf +max.bf 1.bf +min.bf 0.bf -0.bf -min.bf -1.bf -max.bf -inf.bf +nan.bf))
+          (list
+           bf/
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -0.bf (bf "-4.7651298097759032e-323228497") -1.bf -inf.bf -inf.bf
+                 +inf.bf +inf.bf 1.bf (bf "4.7651298097759032e-323228497") 0.bf +nan.bf))
+          (list
+           bfabs
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +inf.bf +max.bf 1.bf +min.bf 0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf))
+          (list
+           bflog
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf
+                 (bf #e-744261117.95489299) 0.bf (bf #e744261117.26174581) +inf.bf +nan.bf))
+          (list
+           bflog2
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf
+                 (bf -1073741824) 0.bf (bf 1073741823) +inf.bf +nan.bf))
+          (list
+           bflog10
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf
+                 (bf #e-323228496.62295526) 0.bf (bf #e323228496.32192528) +inf.bf +nan.bf))
+          (list
+           bflog1p
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf
+                 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf -inf.bf -min.bf -0.bf 0.bf +min.bf
+                 (bf #e0.69314718055994529) (bf #e744261117.26174581) +inf.bf +nan.bf))
+          (list
+           bfexp
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list 0.bf 0.bf (bf #e0.36787944117144233) 1.bf 1.bf
+                 1.bf 1.bf (bf #e2.7182818284590451) +inf.bf +inf.bf +nan.bf))
+          (list
+           bfexp2
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list 0.bf 0.bf (bf #e0.5) 1.bf 1.bf 1.bf 1.bf 2.bf +inf.bf +inf.bf +nan.bf))
+          (list
+           bfexp10
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list 0.bf 0.bf (bf #e0.1) 1.bf 1.bf 1.bf 1.bf 10.bf +inf.bf +inf.bf +nan.bf))
+          (list
+           bfexpm1
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -1.bf -1.bf (bf #e-0.63212055882855767) -min.bf -0.bf
+                 0.bf +min.bf (bf #e1.7182818284590453) +inf.bf +inf.bf +nan.bf))
+          (list
+           bfsin
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e-0.8414709848078965) -min.bf -0.bf
+                 0.bf +min.bf (bf #e0.8414709848078965) +nan.bf +nan.bf))
+          (list
+           bfcos
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e0.54030230586813977) 1.bf 1.bf
+                 1.bf 1.bf (bf #e0.54030230586813977) +nan.bf +nan.bf))
+          (list
+           bftan
+           (list -inf.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e-1.5574077246549023) -min.bf -0.bf
+                 0.bf +min.bf (bf #e1.5574077246549023) +nan.bf +nan.bf))
+          (list
+           bfcsc
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e-1.1883951057781212) -inf.bf -inf.bf
+                 +inf.bf +inf.bf (bf #e1.1883951057781212) +nan.bf +nan.bf))
+          (list
+           bfsec
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e1.8508157176809257) 1.bf 1.bf
+                 1.bf 1.bf (bf #e1.8508157176809257) +nan.bf +nan.bf))
+          (list
+           bfcot
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf (bf #e-0.64209261593433076) -inf.bf -inf.bf
+                 +inf.bf +inf.bf (bf #e0.64209261593433076) +nan.bf +nan.bf))
+          (list
+           bfasin
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf (bf #e-1.5707963267948966) -min.bf -0.bf
+                 0.bf +min.bf (bf #e1.5707963267948966) +nan.bf +nan.bf +nan.bf))
+          (list
+           bfacos
+           (list -inf.bf -max.bf
+                 -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf
+                 +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf
+                 (bf #e3.1415926535897931) (bf #e1.5707963267948966) (bf #e1.5707963267948966)
+                 (bf #e1.5707963267948966) (bf #e1.5707963267948966) 0.bf
+                 +nan.bf +nan.bf +nan.bf))
+          (list
+           bfatan
+           (list -inf.bf -max.bf
+                 -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf
+                 +max.bf +inf.bf +nan.bf)
+           (list (bf #e-1.5707963267948966) (bf #e-1.5707963267948966)
+                 (bf #e-0.78539816339744828) -min.bf -0.bf
+                 0.bf +min.bf (bf #e0.78539816339744828)
+                 (bf #e1.5707963267948966) (bf #e1.5707963267948966) +nan.bf))
+          (list
+           bfsinh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -inf.bf -inf.bf (bf #e-1.1752011936438014) -min.bf -0.bf
+                 0.bf +min.bf (bf #e1.1752011936438014) +inf.bf +inf.bf +nan.bf))
+          (list
+           bfcosh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +inf.bf +inf.bf (bf #e1.5430806348152437) 1.bf 1.bf
+                 1.bf 1.bf (bf #e1.5430806348152437) +inf.bf +inf.bf +nan.bf))
+          (list
+           bftanh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -1.bf -1.bf (bf #e-0.76159415595576485) -min.bf -0.bf
+                 0.bf +min.bf (bf #e0.76159415595576485) 1.bf 1.bf +nan.bf))
+          (list
+           bfcsch
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -0.bf -0.bf (bf #e-0.85091812823932156) -inf.bf -inf.bf
+                 +inf.bf +inf.bf (bf #e0.85091812823932156) 0.bf 0.bf +nan.bf))
+          (list
+           bfsech
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list 0.bf 0.bf (bf #e0.64805427366388535) 1.bf 1.bf
+                 1.bf 1.bf (bf #e0.64805427366388535) 0.bf 0.bf +nan.bf))
+          (list
+           bfcoth
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -1.bf -1.bf (bf #e-1.3130352854993312) -inf.bf -inf.bf
+                 +inf.bf +inf.bf (bf #e1.3130352854993312) 1.bf 1.bf +nan.bf))
+          (list
+           bfasinh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -inf.bf (bf #e-744261117.95489299) (bf #e-0.88137358701954305) -min.bf -0.bf
+                 0.bf +min.bf (bf #e0.88137358701954305) (bf #e744261117.95489299) +inf.bf +nan.bf))
+          (list
+           bfacosh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf
+                 +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf +nan.bf +nan.bf +nan.bf 0.bf
+                 (bf #e744261117.95489299) +inf.bf +nan.bf))
+          (list
+           bfatanh
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf -inf.bf -min.bf -0.bf 0.bf +min.bf +inf.bf +nan.bf +nan.bf +nan.bf))
+          (list
+           bfeint
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf
+                 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf (bf #e-744261117.37767732)
+                 (bf #e1.8951178163559368) +inf.bf +inf.bf +nan.bf))
+          (list
+           bfli2
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list -inf.bf (bf -276962305333851100) (bf #e-0.8224670334241132) -min.bf -0.bf
+                 0.bf +min.bf (bf #e1.6449340668482264) (bf -276962305333851100) -inf.bf +nan.bf))
+          (list
+           bfgamma
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf 1.bf +max.bf +inf.bf
+                 +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf +inf.bf +inf.bf 1.bf +inf.bf +inf.bf
+                 +nan.bf))
+          (list
+           bfpsi0
+           (list -inf.bf -1.bf -min.bf -0.bf 0.bf +min.bf
+                 1.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +inf.bf +inf.bf -inf.bf -inf.bf
+                 (bf #e-0.57721566490153287) +inf.bf +nan.bf))
+          (list
+           bfzeta
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list +nan.bf 0.bf (bf #e-0.083333333333333329) (bf #e-0.5) (bf #e-0.5)
+                 (bf #e-0.5) (bf #e-0.5) +inf.bf 1.bf 1.bf +nan.bf))
+          (list
+           bferf
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +max.bf +inf.bf
+                 +nan.bf)
+           (list -1.bf -1.bf (bf #e-0.84270079294971489) (bf "-2.6884366029284653e-323228497") -0.bf
+                 0.bf (bf "2.6884366029284653e-323228497") (bf #e0.84270079294971489) 1.bf 1.bf
+                 +nan.bf))
+          (list
+           bferfc
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 
+                 0.bf +min.bf 1.bf +max.bf +inf.bf +nan.bf)
+           (list 2.bf 2.bf (bf #e1.8427007929497148) 1.bf 1.bf
+                 1.bf 1.bf (bf #e0.15729920705028513) 0.bf 0.bf +nan.bf))
+          (list
+           bfbesj0
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list 0.bf (bf #e0.76519768655796661) 1.bf 1.bf
+                 1.bf 1.bf (bf #e0.76519768655796661) 0.bf +nan.bf))
+          (list
+           bfbesj1
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list 0.bf (bf #e-0.4400505857449335) -0.bf -0.bf
+                 0.bf 0.bf (bf #e0.4400505857449335) 0.bf +nan.bf))
+          (list
+           (λ (x) (bfbesj 0 x))
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list 0.bf (bf #e0.76519768655796661) 1.bf 1.bf
+                 1.bf 1.bf (bf #e0.76519768655796661) 0.bf +nan.bf))
+          (list
+           (λ (x) (bfbesj 1 x))
+           (list -inf.bf -1.bf -min.bf -0.bf
+                 0.bf +min.bf 1.bf +inf.bf +nan.bf)
+           (list 0.bf (bf #e-0.4400505857449335) -0.bf -0.bf
+                 0.bf 0.bf (bf #e0.4400505857449335) 0.bf +nan.bf))
+          (list
+           bfbesy0
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf
+                 (bf #e-473811343.56828988) (bf #e0.088256964215676956) 0.bf +nan.bf))
+          (list
+           bfbesy1
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf
+                 1.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf -inf.bf
+                 (bf #e-0.78121282130028868) 0.bf +nan.bf))
+          (list
+           (λ (x) (bfbesy 0 x))
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf
+                 +min.bf 1.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf
+                 (bf #e-473811343.56828988) (bf #e0.088256964215676956) 0.bf +nan.bf))
+          (list
+           (λ (x) (bfbesy 1 x))
+           (list -inf.bf -max.bf -1.bf -min.bf -0.bf 0.bf +min.bf
+                 1.bf +inf.bf +nan.bf)
+           (list +nan.bf +nan.bf +nan.bf +nan.bf -inf.bf -inf.bf -inf.bf
+                 (bf #e-0.78121282130028868) 0.bf +nan.bf)))])
+    (match-define (list f xs ys) f+xs+ys)
+    (for ([x  (in-list xs)]
+          [y  (in-list ys)])
+      (define y0 (f x))
+      (unless (equal? y y0)
+        (printf "f = ~a  x = ~v  y = ~v~n" f x y))
+      (check-equal? y0 y))))