Machine-precision numbers cause loss of precision

Question

I am trying to use the Mod function as follows.

Table[
  Nest[Mod[2 #, 1] &, FractionalPart[Pi], n], 
  {n, 100}] // N

Giving the following result:

{0.283185,0.566371,0.132741,0.265482,0.530965,0.0619298,0.12386,0.247719,0.495439,0.990877,0.981755,0.963509,0.927018,0.854036,0.708073,0.416146,0.832291,0.664583,0.329165,0.658331,0.316661,0.633322,0.266645,0.533289,0.0665783,0.133157,0.266313,0.532626,0.0652523,0.130505,0.261009,0.522018,0.0440369,0.0880737,0.176147,0.352295,0.70459,0.40918,0.818359,0.636719,0.273438,0.546875,0.09375,0.1875,0.375,0.75,0.5,0.,0.,0.,0.,0.,-1.,-2.,-4.,-8.,-17.,-35.,-70.,-141.,-282.,-564.,-1129.,-2259.,-4518.,-9036.,-18072.,-36145.,-72290.,-144580.,-289161.,-578323.,-1.15665*10^6,-2.31329*10^6,-4.62658*10^6,-9.25317*10^6,-1.85063*10^7,-3.70127*10^7,-7.40254*10^7,-1.48051*10^8,-2.96101*10^8,-5.92203*10^8,-1.18441*10^9,-2.36881*10^9,-4.73762*10^9,-9.47525*10^9,-1.89505*10^10,-3.7901*10^10,-7.5802*10^10,-1.51604*10^11,-3.03208*10^11,-6.06416*10^11,-1.21283*10^12,-2.42566*10^12,-4.85133*10^12,-9.70265*10^12,-1.94053*10^13,-3.88106*10^13,-7.76212*10^13,-1.55242*10^14}

Now when I do the following

ReleaseHold[
 Hold[Table[
   Nest[Mod[2. #1, 1.] &,     
     FractionalPart[3.14159], n], {n, 100.}]]]

I get the following

{0.283185,0.566371,0.132741,0.265482,0.530965,0.0619298,0.12386,0.247719,0.495439,0.990877,0.981755,0.963509,0.927018,0.854036,0.708073,0.416146,0.832291,0.664583,0.329165,0.658331,0.316661,0.633322,0.266645,0.533289,0.0665783,0.133157,0.266313,0.532626,0.0652523,0.130505,0.261009,0.522018,0.0440369,0.0880737,0.176147,0.352295,0.70459,0.40918,0.818359,0.636719,0.273438,0.546875,0.09375,0.1875,0.375,0.75,0.5,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}```

Clearly something to do with working precision. But could you explain why this discrepancy ?

Answer 1

The machine-precision series

A machine precision number has the exact form $a \times 2^{-b}$, where $a$ and $b$ are integers. [The exponent $-b$ is written with a minus sign for convenience in what follows.] Multiplying by the machine precision number 2. is an exact operation (no round-off error). If $a$ and $b$ are positive and $a$ has the binary expansion $+ a_{b}a_{b-1}\cdots a_{1}$, which would be the case if $a$ were the result of FractionalPart[] or Mod[x, 1.], then Mod[2. a, 1.] would result in $$ a_{b-1}\cdots a_{1} \times 2^{-(b-1)} $$ Note that the length of Mod[2. a, 1.] in bits and the exponent in $2^{-(b-1)}$ would be diminished by one. After $b-1$ more iterations, the result will be zero.

Computational evidence

res = ReleaseHold[
   Hold[Table[
     Nest[Mod[2. #1, 1.] &, FractionalPart[3.14159], n], {n, 100.}]]];

Count[res, x_ /; x > 0.]

(*  49  <-- number of nonzero results *)

SetPrecision[FractionalPart[3.14159], Infinity]
Log2@ Denominator[%]

(*

  159416167809847/1125899906842624
  50  <-- on the 50th step, we should get 0.
*)

The exact series

The step in the iteration Mod[2 #, 1] & either (1) multiplies by 2 or (2) multiplies by 2 and subtracts 1. In the exact series, the result is never zero, so the effect is that we'll have a sequence of multiplying by 2 following by subtracting one. In short, the current result is multiplied by the smallest power of 2 to make the result between 1 and 2 and then 1 is subtracted. Here's the result after 13 iterations:

-1 + 2  (
 -1 + 2  (
  -1 + 32 (
   -1 + 8  (
    -1 + 8  *
       (-3 + Pi)))))

Applying N introduces a rounding error in approximating Pi - 3. The rounding error is propagated until 0. is reached as in the machine-precision series. At that point 0. is multiplied by a power of 2. and 1. is subtracted. And so forth. The rounding error is larger than one might suspect because N[Pi - 3] seems to be computed as N[Pi] - N[3], and the subtractive cancellation yields a relative error of around 2^{-51} instead of 2^{-53}. Coincidentally, the absolute error itself is 2^{53}:

N[Pi - 3] - N[Pi - 3, 16]
Log2[Abs@%]
(*
  -1.11022*10^-16
  -53.
*)

Consequently, after 53 iterations, the -1 shows up. After that, multiplying by a power of 2 and subtracting one yields larger and larger negative results.