# How does computer store non-repeating decimal very accurately?

Maybe it is not accurate at all, but it looks very accurate for me. I'm making a scientific calculator for my portfolio, and I found an interesting phenomenon.

Because computers store decimals not accurately but approximately, I thought that asin(sin(30rad)) wouldn't give 30rad, and sqr(sqrt(7)) wouldn't give 7, etc..

However, it seems that every kind of calculations like above gives an accurate result. (I used cmath header of c++, and I checked that python math modules also gives the same)

This consequence seems even arcane to me. Does the computer really store those non-repeating deicmals? or is it just a matter of implementation of math libraries, like substituting a special symbol for them?

• Read the Wikipedia article on floating point. – Yuval Filmus Jul 28 '18 at 6:34
• The results you see are just an illusion. They are rounded, hiding from you small rounding errors. – Yuval Filmus Jul 28 '18 at 6:35
• The last sentence of your question hints at something called "symbolic computation". In your case c++/python are just using floating point and you're "lucky" to get the right answer, but there do exist systems that can evaluate such expression exactly. – Tom van der Zanden Jul 28 '18 at 8:26
• Here is also proof that C++ doesn't compute the exact answer. It only looks like it: https://ideone.com/ZwNbds – Jakube Jul 28 '18 at 15:43

It may be that sqr (sqrt (7)) is displayed as 7, but it isn't actually exactly equal to 7. That's something you need to check. What you see is not always what you get. It may be that sqr (sqrt (7)) is exactly 7, by "coincidence" (not really coincidence, more like "unpredictable with my limited knowledge").

Take any floating point number x, 1 ≤ x < 4. When the square root is calculated, you get another floating point number y, 1 ≤ y ≤ 2. There are fewer floating point numbers (about half as many) 1 ≤ y ≤ 2 as floating point numbers 1 ≤ x < 4. If you square two numbers y and y', that calculation doesn't know how y and y' were created. So there are about half as many floating point numbers sqr (sqrt (x)), 1 ≤ x ≤ 4, as original numbers in that interval. So for about half the numbers, sqr (sqrt (x)) ≠ x.

The other way round, calculating sqrt (sqr (x)), you can probably prove that the result must be equal to x if both operations are calculated with high quality. If the squaring operation produces a relative error epsilon, then the square root operation divides that relative error by two.

PS. Yuval checked that for x = 7, sqr (sqrt (x)) is not equal to x. If you checked all integers from 1 to 10,000 except those that are squares, I'd expect that sqr (sqrt (x)) equals x in about half the cases that you check.

• Indeed, 5105 of the integers from 1 to 10000 satisfy math.sqrt(n)**2==n. – Yuval Filmus Jul 28 '18 at 15:40
• @gnasher729 It probably would make sense to mention the term "pigeonhole principle" in addition to the current textual description. – njuffa Jul 28 '18 at 17:20
• @njuffa I wouldn't really like to do that. It feels inappropriate. – gnasher729 Jul 28 '18 at 23:12
• @gnasher729 I am not quibbling with the choice, but I am puzzled in which sense this would be inappropriate? Am I overlooking something? Trying to fit the roots of $2^{n}$ floating-point numbers into $2^{n-1}$ available encodings seems like a textbook example of that principle, i.e. we know there must be collisions. – njuffa Jul 28 '18 at 23:48

If you compute the square root of 7 and then square it, you won’t get exactly 7. Here is what I get in python:

>>> math.sqrt(7)**2

7.000000000000001

And here is what you get if you print the value instead:

>>> print math.sqrt(7)**2

7.0

What happened? When printing a floating point value, the default precision is 12 decimal digits (in python). Since math.sqrt(7)**2 differs from 7 in the 15th digit, when rounding to 12 digits we just get 7, giving you the illusion that there is no error.

Let's play with "bc -l"

• IEEE single : 1 (sign) + 8 (exp) + 23 (mantissa) [+1 hidden]

• IEEE double : 1 (sign) + 11 (exp) + 52 (mantissa ) [+1 hidden]

sqrt(7) = 10.10100101010011111111010100111010010111110001110100110110111100011\ 10011101010011111100110000111111100001101111010001000001101010101001... (binary)

With single precision, we store 23bits :

SQ7S = 10.1010010101001111111101 [next truncated bits are 010011...] (no additional rounding needed)

With double precision, we store 52 bits :

SQ7d = 10.101001010100111111110101001110100101111100011101001 [next trucated bits are 1011...]

Due to the default "round to nearest-even" rounding rule, this double precision number will be rounded up to :

SQ7D = 10.101001010100111111110101001110100101111100011101010

Calculating the square :

SQ7S * SQ7S = 110.111111111111111111111 [00100...] = 6.999999523... (decimal)

SQ7D * SQ7D = 111.00000000000000000000000000000000000000000000000000 [1011...]

rounded to : SQ7D * SQ7D = 111.00000000000000000000000000000000000000000000000001 = 7.000000000000000888... (decimal)

So : Neither single nor double precision IEEE float gives exact result to sqr(sqrt(7)), but the double precision result is closer, of course.

(I have also tried with x87-style long double, 64bits mantissa, and it is still wrong)