# Understanding simple numerical calculation

I am trying to understand why N[0.1]//FullForm returns 0.1' (A). Indeed, $0.1$ in base two is $0.000110011001100110011...$. If I truncate this number to the first 52 digits and convert it back to base ten, I get:

Table[1/2^(4*i) + 1/2^(4*i + 1), {i, 1, 12}] // Total // FullForm // N
(* 0.09999999999999964 denoted by (B) *)


which is different from the obtained result (A).

However, computing the following operation returns the expected result (B)!

N[10.1 - 10)] // FullForm
(* 0.09999999999999964 *)


So, what is the difference between the N[0.1] // FullForm and N[10.1 - 10] // FullForm and why does only the second one return (B)?

• I think this question boils down to the difference between e.g. N[0.1] and N[1.1 - 1]. The division itself is not the source of the problem. – Marius Ladegård Meyer Jan 17 '17 at 22:11
• @MariusLadegårdMeyer Right, edited. – anderstood Jan 17 '17 at 22:16
• How did you calculate the binary representation? There is no such thing as a dot in binary formats. The binary representation of floating point numbers is via significant and exponent. For 0.1 it is 1x10^-1, that is, Both significant (1) and exponent (-1) are exactly represented in a 64 bit float. – Felix Jan 17 '17 at 22:16
• @Felix Consider $1100110011...11\times 2^{-48}$ if you prefer (and note that BaseForm[0.1, 2] returns "0.00011001100110011001101" in MMA ;)). – anderstood Jan 17 '17 at 22:21
• ah.. you can also do N[ 1.1 - 1. ] // FullForm to see the issue. (You might want to put such into the question. ) – george2079 Jan 17 '17 at 22:41

## Introduction

This seems to be a question more about IEEE 754 binary64 format than about Mathematica per se.

The significand is 53 bits (52 stored, since the leading bit is assumed to be 1 in a normal number). When the input "0.1" is converted to a number, presumably, at some point 1. is divided by 10. The OP used the term "truncate," but more precisely, the result should be rounded to the nearest 53-bit floating-point number.

Some tools one can use to explore the binary representation of machine precision numbers are

SetPrecision[x, Infinity]
RealDigits[x, 2, 53]


SetPrecision[x, Infinity] converts the machine real to the exact fraction represented by the floating point number x. RealDigits[x, 2, 53] shows the values of the bits. If you replace the 53 with 54 or higher, the last "bits" will be returned as Indeterminate. Below I plot the bits so that a 1 is red, 0 is white, and Indeterminate is gray. The 53-bit limit is indicated with a vertical grid line.

## OP's examples

In the OP's example Table[] command, the computed fraction is truncated too soon. Here is a table showing what is going on. The first three rows are the OP's table carried out to 12 iterations (as in the OP) up to 14. One has to go up to 14 to get all the bits needed to approximate 0.1 to machine precision. The last three rows show the bits of the machine real 0.1, the result of SetPrecision[0.1, Infinity], and the exact fraction 1/10. One can see that rounding 1/10 results in a carry at the last bit. ClearAll[tenth, mantissaplot, frexp];
tenth[n_] := Table[1/2^(4*i) + 1/2^(4*i + 1), {i, 1, n}] // Total;
mantissaplot[{digits_, exp_}, ref_: 0.1] := {ArrayPlot[
{digits}, ColorRules -> {1 -> Red, Indeterminate -> Gray},
Mesh -> True, MeshStyle -> Directive[Thin, Black],
ImageSize -> 450, Axes -> {True, False},
GridLines -> {{exp - Ceiling@Log2[ref], 53}, None}],
exp};
frexp[x_, bits_: 54] := mantissaplot@RealDigits[x, 2, bits];
TableForm[
Join[
Table[frexp[tenth[n], 57], {n, 12, 14}],
{frexp[0.1, 57],
frexp[SetPrecision[0.1, Infinity], 57],
frexp[1/10, 57]}
],
TableHeadings -> {{12, 13, 14, 0.1, Subscript[0.1, Infinity],
1/10}, {"significand", "exp"}}]


Now let's take up the OP example 10.1 - 10. In the table below, we can see that to store 10.1 the part representing 0.1 is shifted over 7 bits, so that the last 7 bits in 0.1 above are dropped. When 10. is subtracted, the result is shifted back, but the 7 bits are already lost. (This what is meant by precision loss due to subtractive cancellation.) TableForm[{
frexp[10.1],
frexp[10.1 - 10.],
frexp[0.1]
},
TableHeadings -> {{10.1, HoldForm[10.1 - 10.], 0.1}, {"significand", "exp"}}]


Update:

My original answer is based on the apparently wrong assumption that the numbers are represented by $s\times 10^e$ with significand $s$ and exponent $e$. The correct binary representation is $s\times 2^e$ with $s$ having 52 bits plus one bit for the sign and $e$ being 11 bits long for 64 bit floating point numbers. Therefore, 0.1 is not exactly represented by this format.

All your numbers are exactly represented as floating point numbers. There are certain operations that preserve this precision. It seems not changing the exponent is one of them:

N[0.2 - 0.1] // FullForm
(* 0.1 *)


and only changing the exponent is one of them:

N[10/0.1] // FullForm
(* 100. *)


However, changing both in one operation may lead to numerical errors:

N[1.1-1.0] // FullForm
(* 0.10000000000000009 *)
`