Inverse confusion about floating point precision

[Bengt Richter]

On 9 May 2005 11:06:22 -0700, "Dan Bishop" <danb_83 at yahoo.com> wrote:

>Skip Montanaro wrote:
>> I understand why the repr() of float("95.895") is
>"95.894999999999996".
>> What I don't understand is why if I multiply the best approximation
>to
>> 95.895 that the machine has by 10000 I magically seem to get the lost
>> precision back.  To wit:
>>
>>     % python
>>     Python 2.3.4 (#12, Jul  2 2004, 09:48:10)
>>     [GCC 3.3.2] on sunos5
>>     Type "help", "copyright", "credits" or "license" for more
>information.
>>     >>> 95.895
>>     95.894999999999996
>>     >>> 95.895 * 10000
>>     958950.0
>>
>> Why isn't the last result "958949.99999999996"?  IOW, how'd I get
>back the
>> lost bits?
>
>You were just lucky.
>
>The floating-point representation of 95.895 is exactly 6748010722917089
>* 2**-46.
>
>Multiplying by 10000 gives you 67480107229170890000 * 2**-46.  But
>floats can have only 53 significant bits, so this gets normalized to
>8237317776998399.658203125 * 2**-33 and rounded to 8237317776998400 *
>2**-33, which happens to be exactly equal to 958950.
>
>For analogy, consider a decimal calculator with only 3 significant
>digits.  On this calculator, 1/7=0.143, an error of 1/7000.
>Multiplying 0.143 by 7 gives 1.001, which is rounded to 1.00, and so
>you get an exact answer for 1/7*7 despite roundoff error in the
>intermediate step.
>

In bits, the above appears as

 >>> prb(95.895)
 '1011111.1110010100011110101110000101000111101011100001'
 >>> len(prb(95.895).split('.')[1])
 46
 >>> prb(95.895*2**46)
 '10111111110010100011110101110000101000111101011100001'
 >>> int(prb(95.895*2**46),2)
 6748010722917089L
 >>> int(prb(95.895*2**46),2)*10000
 67480107229170890000L
 >>> prb(int(prb(95.895*2**46),2)*10000)
 '111010100001111001011111111111111111111111111111111111010100010000'
 >>> prb(int(prb(95.895*2**46),2)*10000)[:53]
 '11101010000111100101111111111111111111111111111111111'
 >>> int(prb(int(prb(95.895*2**46),2)*10000),2)/2.**46
 958950.0
 >>> prb(int(prb(int(prb(95.895*2**46),2)*10000),2)/2.**46)
 '11101010000111100110'

Regards,
Bengt Richter