Arbitrary-Precision Floating-Point Numbers

7.2. Arbitrary-Precision Floating-Point Numbers#

Just as BigInt handles oversized integers, the BigFloat type manages floating-point numbers with far more precision than the standard Float64.

This is essential because Float64 numbers are stored in binary, and many decimal fractions (like 0.1 or 2.1) cannot be represented perfectly, leading to tiny rounding errors. BigFloat allows us to increase the number of bits used for storage, minimizing these errors and enabling high-precision calculations.

7.2.1. The Float64 Conversion Trap#

Like with BigInt, creating a BigFloat using big() on a standard number can be tricky. The input number is first interpreted as a Float64 before being converted. If the number can be represented perfectly in binary, the result is what you’d expect.

# This works perfectly because 1.75 is exactly 1 + 1/2 + 1/4,
# which has a finite binary representation.
big(1.75)

1.75

However, if the number cannot be represented exactly as a Float64, the initial representation error is carried over and exposed by BigFloat’s high precision.

# The number 2.1 is first rounded to the nearest Float64.
# BigFloat then faithfully represents that small inaccuracy with many digits.
big(2.1)

2.100000000000000088817841970012523233890533447265625

The solution, once again, is to use the big"..." string macro. This tells Julia to parse the decimal string directly into a BigFloat at maximum precision, bypassing the Float64 intermediate step entirely.

# The string macro correctly parses the decimal "2.1" to the highest possible precision.
big"2.1"

2.099999999999999999999999999999999999999999999999999999999999999999999999999986

7.2.2. Controlling Precision with setprecision#

You can control the precision (the number of bits used for the significand) of BigFloat calculations using the setprecision function.

While you can set this globally, the recommended practice is to use a do block. This temporarily sets the precision for only the code inside the block, ensuring your high-precision code doesn’t have unintended side effects elsewhere.

# Temporarily set the precision to 512 bits for the duration of this block.
setprecision(512) do
    big(1) / big(3)
end

0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333346

Thanks to multiple dispatch, many of Julia’s built-in mathematical functions and constants have specialized methods for BigFloat. When you use them, they automatically operate at the currently set precision.

# big(pi) computes π to the current BigFloat precision.
display(big(pi))

# The result of sin(π) should be 0. The tiny non-zero result here is the numerical error
# at this high level of precision—still incredibly close to zero!
display(sin(big(pi)))

3.141592653589793238462643383279502884197169399375105820974944592307816406286198

1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77

7.2.3. Example: Visualizing Cubic Convergence with Halley’s Method#

Halley’s method is a root-finding algorithm that extends Newton’s method. It’s known for its cubic convergence, which means the number of correct digits roughly triples with each iteration.

\[ x_{n+1} = x_n - \frac{2f(x_n) f'(x_n)}{2[f'(x_n)]^2 - f(x_n)f''(x_n)} \]

This rapid convergence is often impossible to observe with standard Float64 because you hit the precision limit almost instantly. But with BigFloat, we can use it like a mathematical microscope to watch this phenomenon in action. Let’s try to compute \(x = \sqrt[3]{a}\) by finding the root of the equation \(f(x) = x^3 - a = 0\).

# Note: This function is written generically. It will work with Float64, BigFloat,
# or any other number type that supports the required arithmetic operations.
function halley_cuberoot(a, ε)
    # Initial guess is 1, but cast to the same type as the input 'a'.
    x = one(a)
    
    # Initialize the change in x to be larger than the tolerance to ensure the loop runs.
    Δx = 2ε
    
    while Δx ≥ ε
        # Apply Halley's formula for f(x) = x^3 - a
        # f'(x) = 3x^2, f''(x) = 6x
        fx = x^3 - a
        f_prime_x = 3x^2
        f_double_prime_x = 6x
        x_new = x - (2*fx*f_prime_x) / (2*f_prime_x^2 - fx*f_double_prime_x)
        
        Δx = abs(x_new - x)
        x = x_new
        
        # This print statement is our "microscope" to observe the error shrinking.
        println("Error = ", abs(x - cbrt(a)))
    end
    return x
end

halley_cuberoot (generic function with 1 method)

# Let's find the cube root of 2.
# By passing in big(2), we trigger the BigFloat methods throughout our function.
# The high-precision tolerance (1e-60) is only meaningful because we're using BigFloat.
halley_cuberoot(big(2), 1e-60)

Error = 0.009921049894873164767210607278228350570251464701507980081975112155299676513956221
Error = 4.149742382441322899723575934299353308297808730594470544772346647558790369568435e-07
Error = 3.001136168965729242876474785354871654493618340778791794290950778604895701742477e-20
Error = 1.135217767765559947507625958119369715609728496052297379270293469389407869634748e-59
Error = 0.0
Error = 0.0

1.259921049894873164767210607278228350570251464701507980081975112155299676513956
Contents