I begin by defining a large integer n
:
Prelude> let n = 5705979550618670446308578858542675373983
Prelude> n :: Integer
5705979550618670446308578858542675373983
Next I looked at the behavior of s1
and s2
:
Prelude> let s1 = (sqrt (fromIntegral n))^2
Prelude> let s2 = (floor(sqrt(fromIntegral n)))^2
Prelude> s1 == fromIntegral n
True
Prelude> s1 == fromIntegral s2
True
Prelude> (fromIntegral n) == (fromIntegral s2)
False
Since any fractional part might be discarded, equality on the last 2 expressions was not expected. However, I didn't expect equality to be intransitive (e.g. n == s1, s1 == s2
, but n != s2
.)
Furthermore, floor
appears to lose precision on the integer part, despite retaining 40 significant digits.
Prelude> s1
5.70597955061867e39
Prelude> s2
5705979550618669899723442048678773129216
This lost precision becomes obvious when testing subtraction:
Prelude> (fromIntegral n) - s1
0.0
Prelude> (fromIntegral n) - (fromIntegral s2)
546585136809863902244767
Why does floor
lose precision, and how is this violating transitivity of equality (if at all)?
What is the best approach to computing floor . sqrt
without loss of precision?
It’s not floor
that is losing precision, but the conversion from Integer
(an arbitrary-precision integer) to Double
(a floating-point value, which has a limited precision). Accordingly, fromIntegral n :: Double
is no longer the same value as n
.
Double
has a 53-bit mantissa (52 explicitly stored, the leading one implicit), which is approximately equivalent to 16 decimal digits. Thus, only the (approx.) 16 most significant digits of the result are valid. The rest is just noise.
Finally, your first two comparisons compare Double
s; and n
converted into Double
, s2
converted into Double
, and s1
are all equal. In the third comparison, however, n
and s2
are both Integer
s; they can be compared as Integer
s, so calling fromIntegral
on them is a no-op, and their non-converted integer values are different. If you force conversion to Double
, the values become equal again:
Prelude> ((fromIntegral n) :: Double) == ((fromIntegral s2) :: Double)
True