# Machine-precision numbers cause loss of precision

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 ?

• As you said, because In your second example you are working with machine precision numbers instead of arbitrary precision. The very small values close to zero In the middle of your run are rounded to zero at machine precision, but they retain their value at arbitrary precision. – MarcoB Jun 13 at 15:14
• As as aside, your Table[Nest[...], {n, 100}] expression is functionally equivalent to NestList[Mod[2#, 1]&, FractionalPart[Pi], 100]. The only difference is the fact that the latter expression includes your starting value; if you really don’t want that, you can use Rest@NestList[...]. – MarcoB Jun 13 at 15:18
• Calculate with arbitrary-precision then convert to machine precision for display. Look at Table[ Nest[Mod[2 #, 1] &, FractionalPart[Pi], n], {n, 100}] // N[#, 20] & // N – Bob Hanlon Jun 13 at 16:51
• How those negative values appear ? @BobHanlon – q than a Jun 14 at 14:03
• @qthana - The code that I posted in my comment does not produce any negative values. Values from code with known precision-induced errors are by definition meaningless. – Bob Hanlon Jun 14 at 14:18

### 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.

• How the negative values appear after zeros ? – q than a Jun 14 at 14:05
• @qthana The result of the iteration computing the Mod is either to multiply by 2. or to multiply by 2. and subtract 1.. In the exact series, the result is never zero, so the effect is that at each step, the current result is multiplied by the smallest power of 2 to make the result between 1. and 2. and then subtract one. 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. – Michael E2 Jun 15 at 20:53