# Why is 0 vs 0.0 causing a precision problem?

I wasn't sure what to title this Question. I tried calculating different forms of a continued fraction, and I found that I get a correct result when I set d to 0 but a wrong result (in the final two terms) if I set d to 0.0

ClearAll["Global*"];
t = Sqrt[17];
d = 0.0;
n = 11;
a = ConstantArray[0, n];
Do[a[[i]] = Floor[t + d]; t = 1/(t - Floor[t + d]), {i, 1, n}];
a

{4, 8, 8, 8, 8, 8, 8, 8, 8, 7, 4}


I only caught this because I happen to know that this series should end in infinite 8's.

When I try to calculate the form of continued fractions that allows negative terms, I get the correct series from d=1/2 but wrong terms if d=.5. So I gather that this problem has something to do with precision and using rational numbers instead of floats but why is this happening for zero? Is Mathematica adding 0.0 incorrectly?

• The presence of any machine precision number (including 0.0) will result in the calculation being done with machine precision. – Bob Hanlon Feb 4 '19 at 21:33
• By using an extended precision version of 0 you can get a small speedup from using numbers instead of full symbolic expressions but without sacrificing accuracy. Try d=01. – b3m2a1 Feb 4 '19 at 21:47
• Also in terms of programming style this is probably faster: d + NestList[1/(# - Floor[# + d]) &, t, 200] // Floor // RepeatedTiming – b3m2a1 Feb 4 '19 at 21:47

## 2 Answers

Your algorithm is probing the limits of machine precision numbers.

When you add Sqrt[17]+0 you get $$\sqrt{17}$$ with infinite precision:

Sqrt[17] + 0
(* Sqrt[17] *)


When you add Sqrt[17]+0.0 you get much less precision:

Sqrt[17] + 0.0 // FullForm
(* 4.123105625617661 *)


The number you're effectively dealing with in the d=0.0 case is not $$t=\sqrt{17}$$ but

SetPrecision[Sqrt[17] + 0.0, ∞]
(* 4642204239785223/1125899906842624 *)


which is the closest available machine-precision number to $$\sqrt{17}$$.

If you change your code to

t = 4642204239785223/1125899906842624;
d = 0;
n = 11;
a = ConstantArray[0, n];
Do[a[[i]] = Floor[t + d]; t = 1/(t - Floor[t + d]), {i, 1, n}];
a


you get a similar answer to the d=0.0 case:

(* {4, 8, 8, 8, 8, 8, 8, 8, 8, 7, 2} *)


Another way of seeing this is that you're doing a continued-fraction representation, which differs when you convert a number to its machine-precision representation:

ContinuedFraction[Sqrt[17]]
(* {4, {8}} *)
ContinuedFraction[Sqrt[17] // N]
(* {4, 8, 8, 8, 8, 8, 8, 8} *)


Notice that the builtin ContinuedFraction command stops before spitting out "wrong" digits due to machine-precision issues.

Two comments rolled into one.

Try using extended precision numbers like:

d = 01


Also try a different programming style to get things cleaner and a bit faster:

t = Sqrt[17];
d = 01;
d + NestList[1/(# - Floor[# + d]) &, t, 200] // Floor // RepeatedTiming

{0.0055, {4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
8}}