5 Arbitrary-Precision Floating-Point Numbers (Bigfloats)
(require math/bigfloat) | package: math-lib |
A C type of arbitrary-precision floating-point numbers.
Elementary and special functions that are efficient and proved correct.
Well-defined semantics that correspond with the latest IEEE 754 standard.
With a few noted exceptions, bigfloat functions regard their arguments as if they were exact, regardless of their precision. Conceptually, they compute exact results using infinitely many bits, and return results with (bf-precision) bits by rounding them using (bf-rounding-mode). In practice, they use finite algorithms that have been painstakingly proved to be equivalent to that conceptual, infinite process.
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 easy to install on major Unix-like platforms.
5.1 Quick Start
Set the bigfloat function result precision using (bf-precision <some-number-of-bits>).
Use bf to convert real values and well-formed strings to bigfloats.
Operate on bigfloats using bf-prefixed functions like bf+ and bfsin.
Convert bigfloats to real values using bigfloat->real, bigfloat->flonum, and bigfloat->integer. Format them for display using bigfloat->string.
For examples, continue through the FAQ.
5.2 Fictionally Asked Questions
Why use math/bigfloat?
There are a few reasons.
Reason: Flonums have either too much or too little precision for your application.
> (flsqrt 3.0) 1.7320508075688772
> pi 3.141592653589793
> (bf-precision 16) > (bfsqrt (bf 3)) (bf #e1.73206)
> (bf-precision 179) > pi.bf (bf #e3.141592653589793238462643383279502884197169399375105819)
A flonum has a 53-bit significand (we’ll say it has 53 bits of precision) and an 11-bit exponent. A bigfloat has an arbitrary precision of at least 2 bits and a 31-bit exponent.
Reason: To compute ridiculously large or small numbers with confidence.
> (bf-precision 128) > (bfexp (bfexp (bfexp (bf 3)))) (bf "2.050986436051648895105860942072054674579e229520860")
> (bflog (bflog (bflog (bfexp (bfexp (bfexp (bf 3))))))) (bf 3)
Reason: To verify your floating-point hardware.
IEEE 754-2008 stipulates that conforming implementations must correctly round the results of all operations. Roughly speaking, results can’t be more than half a bit off, where the bit in question is the least significant in the significand.
> (bf-precision 53) > (bigfloat->flonum (bfexp (bf 400))) 5.221469689764144e+173
Reason: To control rounding of the least significant bit.
IEEE 754 provides for different rounding modes for the smallest bit of a flonum result, such as round to even and round toward zero. We might use this to implement interval arithmetic correctly, by rounding lower bounds downward and upper bounds upward. But there isn’t a portable way to set the rounding mode!
MPFR allows the rounding mode to be different for any operation, and math/bigfloat exposes this capability using the parameter bf-rounding-mode.
When shouldn’t I use math/bigfloat?
When you need raw speed. Bigfloat functions can be hundreds to thousands of times slower than flonum functions.
That’s not to say that they’re inefficient. For example, bflog implements the algorithm with the best known asymptotic complexity. It just doesn’t run directly on hardware, and it can’t take fixed-precision-only shortcuts.
Why are there junk digits on the end of (bf 1.1)?
That’s approximately the value of the flonum 1.1. Use (bf #e1.1) or (bf "1.1") to make the junk go away. In general, you should prefer to convert exact rationals and strings to bigfloats.
Why is the last digit of pi.bf not rounded correctly?
All the bits but the last is exact, and the last bit is correctly rounded. This doesn’t guarantee that the last digit will be.
A decimal digit represents at most log(10)/log(2) ≈ 3.3 bits. This is an irrational number, so the decimal/bit boundary never lines up except at the decimal point. Thus, the last decimal digit of any bigfloat must represent fewer than 3.3 bits, so it’s wrong more often than not. But it’s the last bit that counts.
5.3 Type and Constructors
> (bf-precision 128) > (bf 4) (bf 4)
> (bf 1/7) (bf #e0.1428571428571428571428571428571428571426)
> (bf 41/10) (bf #e4.099999999999999999999999999999999999995)
> (bf "not a number") bf: expected a well-formed decimal number; given "not a number"
> (bf "15e200000000") (bf "1.499999999999999999999999999999999999998e200000001")
In the last example, the result of (bf "15e200000000") is displayed as a string conversion because the exact rational number would be very large.
* It can be a good idea if you’re testing a flonum implementation of a function against a bigfloat implementation.
> (require (only-in math/base random-bits)) > (bf-precision 64) > (bf (random-bits 64) -64) (bf #e0.416872969910248753552)
> (bf-precision 64)
> (bfcopy (parameterize ([bf-precision (+ (bf-precision) 10)]) (bf/ (bf+ 1.bf (bfsqrt 5.bf)) 2.bf))) (bf #e1.61803398874989484821)
5.4 Accessors and Conversion Functions
procedure
x : Bigfloat
procedure
(bigfloat-signbit x) → (U 0 1)
x : Bigfloat
> (bigfloat-signbit -1.bf) 1
> (bigfloat-signbit 0.bf) 0
> (bigfloat-signbit -0.bf) 1
> (bigfloat-signbit -inf.bf) 1
procedure
(bigfloat-significand x) → Integer
x : Bigfloat
procedure
(bigfloat-exponent x) → Integer
x : Bigfloat
To access the significand and exponent at the same time, use bigfloat->sig+exp.
procedure
(bigfloat->sig+exp x) → (Values Integer Integer)
x : Bigfloat
If (values sig exp) = (bigfloat->sig+exp x), its value as an exact rational is (* sig (expt 2 exp)). In fact, bigfloat->rational converts bigfloats to rationals in exactly this way, after ensuring that (bfrational? x) is #t.
This function and the two-argument variant of bf are mutual inverses.
procedure
(bigfloat->integer x) → Integer
x : Bigfloat
procedure
x : Bigfloat
procedure
(bigfloat->real x) → (U Exact-Rational Flonum)
x : Bigfloat
procedure
(bigfloat->flonum x) → Flonum
x : Bigfloat
bigfloat->integer, bigfloat->rational and bigfloat->real return values that can be converted exactly back to x using bf. For the first two, this is done by raising an error if x is not respectively integer or rational. On the other hand, bigfloat->real returns +inf.0, -inf.0 or +nan.0 when x is not a rational bigfloat.
bigfloat->flonum rounds x to 53 bits precision to fit the value into a flonum, using the current value of bf-rounding-mode.
> (bf-precision 64) > (bigfloat->integer (bf 21/10))
bigfloat->integer: contract violation
expected: bfinteger?
given: (bf #e2.09999999999999999991)
> (bigfloat->integer (bfround (bf 21/10))) 2
> (define x (bf 1/7)) > (bigfloat->flonum x) 0.14285714285714285
> (bigfloat->rational x) 10540996613548315209/73786976294838206464
> (rationalize (bigfloat->rational x) (expt 2 (- (bf-precision)))) 1/7
> (bf= x (bf (bigfloat->rational x))) #t
Be careful with exact conversions. Bigfloats with large exponents may not fit in memory as integers or exact rationals. Worse, they might fit, but have all your RAM and swap space for lunch.
procedure
(bigfloat->string x) → String
x : Bigfloat
procedure
(string->bigfloat s) → (U Bigfloat False)
s : String
The string returned by bigfloat->string includes enough digits that string->bigfloat can reconstruct the bigfloat precisely. In other words, string->bigfloat is a left inverse of bigfloat->string.
If s isn’t a well-formed decimal number with an optional exponent part, string->bigfloat returns #f. (In contrast, (bf s) raises an error.)
> (bf-precision 64) > (bigfloat->string (bf 4)) "4"
> (bigfloat->string (bf 1/10000)) "1.00000000000000000001e-4"
> (string->bigfloat "0.14285714285714285714") (bf #e0.142857142857142857141)
> (string->bigfloat "square root of two") #f
> (string->bigfloat (bigfloat->string pi.bf)) (bf #e3.14159265358979323851)
> pi.bf (bf #e3.14159265358979323851)
5.5 Parameters
parameter
(bf-precision) → Integer
(bf-precision bits) → void? bits : Integer
For nonzero, rational bigfloats, the number of bits bits includes the leading one bit. For example, to simulate 64-bit floating point, use (bf-precision 53) even though flonums have a 52-bit significand, because the one bit is implicit in a flonum.
This parameter has a guard that ensures (bf-precision) is between bf-min-precision and bf-max-precision.
parameter
(bf-rounding-mode) → (U 'nearest 'zero 'up 'down)
(bf-rounding-mode mode) → void? mode : (U 'nearest 'zero 'up 'down)
5.6 Constants
Most bigfloat “constants” are actually identifier macros that expand to the application of a zero-argument function. This allows, for example, pi.bf to depend on the current value of bf-precision, and allows all of them to be constructed lazily. Most constants are memoized, possibly at multiple precisions.
value
value
value
value
value
> (bf-precision 10) > pi.bf (bf #e3.1406)
> (bf-precision 179) > pi.bf (bf #e3.141592653589793238462643383279502884197169399375105819)
> phi.bf bf #e1.618033988749894848204586834365638117720309179805762863)
> gamma.bf (bf #e0.5772156649015328606065120900824024310421593359399235988)
> catalan.bf (bf #e0.9159655941772190150546035149323841107741493742816721343)
> log2.bf (bf #e0.6931471805599453094172321214581765680755001343602552545)
value
value
value
value
value
value
value
value
value
value
The constants -inf.bf, -0.bf, 0.bf, +inf.bf, and +nan.bf have fixed precision.
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
value
5.7 Predicates
procedure
x : Bigfloat
procedure
(bfpositive? x) → Boolean
x : Bigfloat
procedure
(bfnegative? x) → Boolean
x : Bigfloat
procedure
(bfinteger? x) → Boolean
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
(bfrational? x) → Boolean
x : Bigfloat
procedure
(bfinfinite? x) → Boolean
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat y : Bigfloat
5.8 Rounding
procedure
(bftruncate x) → Bigfloat
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
5.9 Mathematical Operations
When given no arguments, bfmin returns +inf.bf, and bfmax returns -inf.bf.
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
When bf+ and bf- are given more than two arguments, they compute the answers in a way that incurs rounding error only once.
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
y : Bigfloat x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
(bfsinh+cosh x) → (Values Bigfloat Bigfloat)
x : Bigfloat
procedure
(bffactorial x) → Bigfloat
x : Integer
procedure
(bflog-gamma x) → Bigfloat
x : Bigfloat
procedure
(bflog-gamma/sign x) → (Values Bigfloat (U -1 1))
x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
n : Integer x : Bigfloat
procedure
x : Bigfloat
procedure
x : Bigfloat
procedure
n : Integer x : Bigfloat
A “j” in the name indicates that a function computes a Bessel function of the first kind. A “y” indicates the second kind.
The “j” or “y” is followed by the order: zero, one, or n (user-specified).
5.10 Low-level Functions
procedure
(bigfloat->ordinal x) → Integer
x : Bigfloat
procedure
(ordinal->bigfloat n) → Bigfloat
n : Integer
procedure
(bigfloats-between x y) → Integer
x : Bigfloat y : Bigfloat
procedure
x : Bigfloat n : Integer
procedure
x : Bigfloat
procedure
x : Bigfloat
The major difference is that these operate using (bf-precision) bits. Additionally, unlike other bigfloat functions, all of these convert their bigfloat arguments to (bf-precision) bits.
procedure
(bfcanonicalize x) → Bigfloat
x : Bigfloat
For zero or non-rational x, returns -inf.bf, -0.bf, 0.bf, +inf.bf, or +nan.bf, depending on the value of x.
Two nonzero, rational bigfloats are equal? if and only if their canonicalized significands and exponents are equal. Two zero or non-rational bigfloats are equal? if and only if their canonicalizations are eq?.
Canonicalizing bigfloats won’t change answers computed from them.
> (bf-precision 64) > (define x (bf 1 -2)) > x (bf #e0.25)
> (bfcanonicalize x) (bf #e0.25)
> (bigfloat-precision x) 64
> (bigfloat-precision (bfcanonicalize x)) 2