Split "mpfr.rkt" into "gmp.rkt" and "mpfr.rkt"

Fixed bigfloat functions that assumed (fixnum? x) means x fits in a _long
(not true on Win64)

Hopefully fixed dangling pointer errors that broke `math/bigfloat' on Win64.
It apparently had no _long/_int mismatches, but GC on Win64 will run between
creating an `_mpz' and using its value after passing it as an output argument
to MPFR functions. That doesn't seem to happen on 64-bit Linux or Mac. No
idea why, but Win64 exposed the problem so... that's good, I guess.

Rewrote `rational->bigfloat' to not use GMP's rationals

More/better bigfloat tests

Added bigfloat stress test w/ weak leak detection

Reenabled custodian shutdown callback that clears MPFR constants, because it
seems to work now

Removed `mpfr-available?' because it would only return non-#f
This commit is contained in:
Neil Toronto 2012-12-03 18:41:12 -07:00
parent bd1a63a43a
commit d935bc0643
7 changed files with 728 additions and 258 deletions

View File

@ -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

View File

@ -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))

View File

@ -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))

View File

@ -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 ()

View File

@ -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)

View File

@ -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))

View File

@ -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))))