Arbitrary precision floating-point numbers
Contents
7.2. Arbitrary precision floating-point numbers¶
Similarly, the BigFloat data type can create floating-point numbers with higher precision than the regular Float64 type. It can also be created from an existing Float64 number:
BigFloat(1.75) # Accidentally OK, since 1.75 is exactly represented by Float64
1.75
But again this might give unintended results, since the number is first converted to Float64 with only about 16 decimal digits of accuracy:
BigFloat(2.1) # First rounded to Float64, low precision
2.100000000000000088817841970012523233890533447265625
Again this can be solved using the string syntax:
parse(BigFloat, "2.1") # High precision
2.099999999999999999999999999999999999999999999999999999999999999999999999999986
The precision, or the number of bits used to store the number, can be controlled with the setprecision function. It will set the precision for all later operations, unless used in a do block as shown below when it will only apply to that block of code.
setprecision(512) do
BigFloat(1) / BigFloat(3)
end
0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333346
Julia’s built-in functions support the higher precision:
display(BigFloat(pi))
display(sin(BigFloat(pi)))
3.141592653589793238462643383279502884197169399375105820974944592307816406286198
1.096917440979352076742130626395698021050758236508687951179005716992142688513354e-77
7.2.1. Example: Cubic convergence¶
Halley’s method is an extension of Newton’s method for solving \(f(x)=0\), which involves the second derivative but gives cubic convergence:
This cubic convergence is sometimes hard to observe, since the machine precision is achieved very quickly. But with arbitrary precision we can confirm the “tripling of the accurate digits” in each iteration. For example, consider computing \(x=\sqrt[3]{a}\) by solving the equation \(f(x) = x^3 - a\):
function halley_cuberoot(a, ε)
x = 1 # Initial guess
Δx = 1 # Need to initialize to get started
while Δx ≥ ε
xnew = x - 2(x^3 - a)*3x^2 / (2(3x^2)^2 - (x^3 - a)*(6x))
Δx = abs(xnew - x)
x = xnew
println("Error = ", abs(x - cbrt(a)))
end
x
end
halley_cuberoot (generic function with 1 method)
halley_cuberoot(BigFloat(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