Nitpicking, there is no such thing as “exactly 0.001” in double (or any) precision IEEE-754, only something whose string representation is exactly 0.001.
7 Likes
In Matlab
>> sprintf('%.20f', 10^-3)
ans =
'0.00100000000000000002'
6 Likes
Note that this is a different thing, it is the Int elevated to the Int. Fortran returns 0 in that case.
To me, 0.0010000000000000002 == 0.001 when I am working with Float64. So I don’t view it as “more” or “less” accurate. We all know about the quirks of floating point arithmetic…
3 Likes
No, it’s not. In Matlab the ints are flints
>> class(10)
ans =
'double'
3 Likes
But I’m confused now
# Julia 1.8beta
julia> @sprintf("%.20f", 10.0^-3)
"0.00100000000000000024"
# Juila 1.7.0
julia> @sprintf("%.20f", 10.0^-3)
"0.00100000000000000002"
Folks, let’s be precise. We are talking about the following binary64 numbers: 0x3F50624DD2F1A9FC (printed as 0.001) and 0x3F50624DD2F1A9FD (printed as 0.0010000000000000002)
julia> reinterpret(Float64, 0x3F50624DD2F1A9FC)
0.001
julia> reinterpret(Float64, 0x3F50624DD2F1A9FD)
0.0010000000000000002
8 Likes
Has an analysis been done on how often the approaches actually differ? I tried a few other numbers and it seems that they sometimes give the same result anyway.
The standard IEEE 754-2019 is quite clear on this question, in Subclause 9.2:
A conforming operation shall return results correctly rounded for the applicable rounding direction for all operands in its domain.
(emphasis added by me)
“Correctly rounded” here means according to the defined rounding mode. All rounding modes (Subclause 4.3) require that
the floating-point number nearest to the infinitely precise result shall be delivered
Beacuse the true result is actually smaller than 0x3F50624DD2F1A9FC, the next float (0x3F50624DD2F1A9FD) is completely invalid. The true result lies between 0x3F50624DD2F1A9FB and 0x3F50624DD2F1A9FC.
7 Likes
There is one digit more than expected there:
julia> @sprintf("%.20f", 10.0^-3)
"0.00100000000000000002"
julia> reinterpret(Float64, 0x3F50624DD2F1A9FD)
0.0010000000000000002
Matlab:
'0.00100000000000000002'
So Matlab is returning the same number other languages return.
1 Like
The accurate approach has average error of .25 ULP, and max error of .5 ULP. The fast approach has average error of .84 ULP, and maximal error of 2 ULP.
7 Likes
As such, after more testing, I’m fairly convinced that the slightly slower answer is the one we want.
22 Likes
I don’t understand how your post actually answers my question about a substantive comparison between the two approaches, but you do seem to suggest that the “less accurate method” is not just less accurate but actually incorrect (according to the standards). If this is indeed true, then we should not even be having this poll. I am not knowledgeable enough on the standards to have that discussion.
1 Like
Subclause 9.2 is specifically for correctly rounded functions (what the next version of C++ will call sqrt_cr, exp_cr etc.) Julia doesn’t guarantee you correct rounding for elementary functions (other than sqrt).
2 Likes
My most common use of x^-3 is in defining SI scaling factors:
julia> const m = 10^-3
0.0010000000000000002
julia> const u = 10^-6
1.0000000000000004e-6
# ...
Then I use 124.2m etc. in many places. So if the m value is less inaccurate then it makes all the other numbers less accurate. The value of 10^-3 is probably in my top 20 of common numbers I use so I really want those to be accurate by default.
4 Likes
Matlab and Julia1.7. Not Julia1.8
Do we have the same issue for other exponents? Using v1.7.2:
julia> 10^-1 == 0.1
true
julia> 10^-2 == 0.01
false
julia> 10^-3 == 0.001
false
julia> 10^-4 == 0.0001
false
julia> 10^-5 == 0.00001
false
julia> 10^-6 == 0.000001
false
# etc
Maybe it is just me being OCD but I feel let down when I get something like 1.0000000000000004e-6. Maybe I want too much.
Edit: I didn’t realize that this was the case. I thought ops were guaranteed rounding to the right number…from now on I guess I type it out the long way :(. … or just 1e-6 but there is a syntactic beauty with 10^-6.
To be fair, if you depend on those constants so much, you should probably be using Unitful.jl (or just typing out m = 0.001).
2 Likes
In Julia 1.7, literal pow (i.e. 10^-4) will fairly consistently give an inaccurate answer, but 10.0^-4 would give the correct answer. In 1.8, all negative powers of 10 other than negative 2 and negative 3 currently give the right answer. 10^-2 will likely remain wrong for somewhat complicated reasons, but 10^-3 can be fixed relatively easily.
Often I just want the SI constants I need (because it pollutes the namespace, taking up single character identifiers) and not actual units.