1

In the ghci terminal, I was computing some equations with Haskell using the sqrt function.

I notice that I would sometimes lose precision in my sqrt result, when it was supposed to be simplified.

For example,

sqrt 4 * sqrt 4 = 4 -- This works well!
sqrt 2 * sqrt 2 = 2.0000000000000004 -- Not the exact result.

Normally, I would expect a result of 2.

Is there a way to get the right simplification result?
How does that work in Haskell?

BenMorel
  • 34,448
  • 50
  • 182
  • 322
kolam
  • 731
  • 4
  • 17
  • http://stackoverflow.com/q/16878251/2348704 has some good answers to what you are looking for – bhaskarc Jun 28 '14 at 12:54
  • 12
    What every programmer/computer scientist should know about float: http://www.validlab.com/goldberg/paper.pdf – Random Dev Jun 28 '14 at 13:07
  • I think this is a decent question. Anyone who downvoted it probably didn't read it thoroughly. – wvdz Jun 28 '14 at 22:34
  • You could do a little bit with rewrite rules if you really wanted to, but you might not get consistent results. – dfeuer Jun 29 '14 at 05:49
  • 1
    Actually we have talked a lot about this with GHC, and it's why implementing some of these optimisations are so hard. The *right* thing to do is to emulate this loss of precision when simplifying, so the optimised programs behave the same as unoptimized programs. This is a hard thing to do when cross compiling, and that's one reason why GHC doesn't do float constant folding yet (I actually claimed the ticket)! – kvanbere Jun 29 '14 at 06:40

3 Answers3

7

There are usable precise number libraries in Haskell. Two that come to mind are cyclotomic and the CReal module in the numbers package. (Cyclotomic numbers don't support all the operations on complex numbers that you might like, but square roots of integers and rationals are in the domain.)

>>> import Data.Complex.Cyclotomic 
>>> sqrtInteger 2
e(8) - e(8)^3
>>> toReal $ sqrtInteger 2 
Just 1.414213562373095     -- Maybe Double
>>> sqrtInteger 2 * sqrtInteger 2
2
>>> toReal $ sqrtInteger 2 * sqrtInteger 2
Just 2.0
>>> rootsQuadEq 3 2 1
Just (-1/3 + 1/3*e(8) + 1/3*e(8)^3,-1/3 - 1/3*e(8) - 1/3*e(8)^3)
>>> let eq x  = 3*x*x + 2*x + 1
>>> eq (-1/3 + 1/3*e(8) + 1/3*e(8)^3)
0
>>> import Data.Number.CReal 
>>> sqrt 2 :: CReal
1.4142135623730950488016887242096980785697 -- Show instance cuts off at 40th place
>>> sqrt 2 * sqrt 2 :: CReal
2.0
>>> sin 3  :: CReal
0.1411200080598672221007448028081102798469
>>> sin 3*sin 3 + cos 3*cos 3 :: CReal
1.0
Michael
  • 2,889
  • 17
  • 16
3

You do not lose precision. You have limited precision.

The square root of 2 is a real number but not a rational number, therefore it's value cannot be represented exactly by any computer (except representing it symbolically, of course).

Even if you define a very large precision type, it will not be able to represent the square root of 2 exactly. You may get more precision, but never enough to represent that value exactly (unless you have a computer with infinite memory, in which case please hire me).

Analog File
  • 5,280
  • 20
  • 23
  • I don't agree that this is a good answer. OP assumes that Haskell would simplify the statement `sqrt 2 * sqrt 2`. This implies that he understands that `sqrt 2` cannot be expressed precisely as a float. – wvdz Jun 28 '14 at 22:30
  • A compiler should always and only perform optimizations when the result would not change. `sqrt` has type `Floating a => a -> a`, `sqrt 2` has type `Floating a => a` and `(*)` has type `Num a => a -> a -> a` therefore `sqrt 2 * sqrt 2`also has type `Floating a => a`. This is a polymorphic result. The actual type used is decided when you instantiate it (for example by printing). The compiler has no (and cannot have) knowledge of exactly what code will be executed. And certainly not about what level of precision is required. It _cannot_ optimize it because the result does not have to be `2`. – Analog File Jun 29 '14 at 07:58
  • 1
    Not only the result could be anything, but in the specific case of a _double precision floating point_ type the result __must not__ be `2`. The result you get is the one you are guaranteed to get by the IEEE 754 standard that defines the floating point type. Any other result would be wrong. – Analog File Jun 29 '14 at 08:02
  • "except representing it symbolically..." the dichotomy floating-point vs. symbolic is common enough, but it's false. There are plenty of Haskell libraries which use efficient representations somewhere in between; the cyclotomic numbers as brought forth by Arthur are a beautiful example. (Oh, I see you've already been discussing that with Lennart and Laurent below. I disagree, cyclotomics are not symbolic, not in a fundamentally different sense than floats are.) – leftaroundabout Jun 29 '14 at 11:20
  • @leftaroundabout any lossless encoding is a symbolic representation, in many encodings the symbols are implicit but they still are there. Floating point encoding is (very) lossy and cannot be called symbolic in any way. Cyclotomic numbers can. And again, I said "can be argued", not exactly a false statement. – Laurent Giroud Jun 29 '14 at 19:23
  • @LaurentGiroud: floating-point is lossless as well if you use it only for the supported set of integers and binary fractions. Cyclotomic numbers are of course much more general and in particular have a much _nicer specification_ of what they can express exactly, but it's still just a subset. Dense in ℂ, but that doesn't mean it's feasible to express any complex number to arbitrary high precision. As to why, consider the rather simpler `Rational`s: would you call these symbolic? Expressing, say, π as "exact" by taking a ridiculously large denominator is not really a good idea. – leftaroundabout Jun 29 '14 at 19:46
  • Of course, this raises the question "what _is_ symbolic then". I'd call it symbolic if the expression is evaluated in its raw form, "outside to inside", trying to make sense of the thing as a whole – only this way can you e.g. solve integrals symbolically, but the price is that it's not really possible to have well-typed, total functional subexpressions. Whereas all these nice Haskell number types will be evaluated "inside to outside", i.e. any subexpression makes the same sense on its own as in the big picture. – leftaroundabout Jun 29 '14 at 19:53
0

The explanation for these results lies in the type of the values returned by the sqrt function:

> :t sqrt
sqrt :: Floating a => a -> a

The Floating a means that the value returned belongs to the Floating type class. The values of all types belonging to this class are stored as floating point numbers. These sacrifice precision for the sake of covering a larger range of numbers.

Double precision floating point numbers can cover very large ranges but they have limited precision and cannot encode all possible numbers. The square root of 2 (√2) is one such number:

> sqrt 2
1.4142135623730951
> sqrt 2 + 0.000000000000000001
1.4142135623730951

As you see above, it is impossible for double precision floating point numbers to be precise enough to represent √2 + 0.000000000000000001, it is simply rounded to the closest approximation which can be expressed using floating point encoding.

As mentioned by another poster, √2 is an irrational number which can be simplified to mean that it requires an infinite number of digits to represent correctly. As such it cannot be represented faithfully using floating point numbers. This leads to errors such as the one you noticed when multiplying it with itself.

You can learn about floating points on their wikipedia page: http://en.wikipedia.org/wiki/Floating_point.

I especially recommend that you read the answer to this other Stack Overflow question: Floating Point Limitations and follow the mentioned link, it will help you understand what's going on under the hood.

Note that this is a problem in every language, not just Haskell. One way to get rid of it entirely is to use symbolic computation libraries but they are much slower than the floating point numbers offered by CPUs. For many computations the loss of precision due to floating points is not a problem.

Community
  • 1
  • 1
Laurent Giroud
  • 207
  • 2
  • 9
  • 1
    Well, symbolic numbers are not the only answer as we have seen. – augustss Jun 29 '14 at 00:19
  • @augustss, it could be argued that cyclotomic numbers are a symbolic representation of complex numbers but I have removed the "only" nevertheless. – Laurent Giroud Jun 29 '14 at 01:27
  • @augustss true. But cyclotomic numbers are also not a general solution. They just happen to be capable of representing the square root of 2, but they still are a subset of the irrationals (reals that are not also rationals). They would have handled the `sqrt 2 * sqrt 2` case correctly in the same sense that plain doubles can handle `sqrt 4 * sqrt 4`: it just happens to be a computation where all intermediate values are finitely representable by the type. Plus I agree with @LaurentGiroud that in a very loose sense cyclotomics are like a symbolic representation of some complex numbers. – Analog File Jun 29 '14 at 08:08
  • @AnalogFile If you call cyclotomic numbers symbolic then you'd also have to call regular complex numbers symbolic. They both extend a number field with constants (sqrt(-1) in the case of complex numbers and roots on unity for cyclotomic numbers) and compute with that. But I wasn't just referring to cyclotomic numbers, I was also referring to constructive real numbers. Depending on your philosophical view these are exactly the real numbers or just a subset. Of course, some operations, like equality test, become semidecidable for constructive reals. – augustss Jun 29 '14 at 12:06
  • @augustss a + ib is a symbolic representation of a complex number, a.e^(i.theta) is another just as √2 is a symbolic representation of "the real number which square is 2". Finally I said "can be argued" and toned down my remark to leave room for people who disagree with my affirmation. – Laurent Giroud Jun 29 '14 at 19:31