Very different results from evaluating same expression with different precisions

Question

When I evalute the following expression,

N[95881665812878 - 120000000000000 (121576521638975 - (321097753837557 Log[1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3))])/ Log[5]) + ( 321097753837557 Log[(-1 + 1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3)))/(-6 + 4/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3)))])/ Log[5]]

I obtain

-1.875*10^12

But, when I evalute it asking for 3-digits, I get

0.487

What is the reason for the different results?

Answer 1

As others have mentioned, the wrong result is given by N[expr] and the errors are due to cancellation. Let's discuss a bit why N[expr, 3] is able to give a good result.

Mathematica can do computations with inexact ( = floating point) numbers in two ways:

• Using the computers native floating point arithmetic, which is very fast, but has no precision tracking. This is what happens when you use N[expr] or the completely equivalent N[expr, MachinePrecision].

• Using Mathematica's own implementation of arbitrary precision arithmetic which does have precision tracking. This is what happens when you use N[expr, precision], where precision is a number. Note that N[expr, $MachinePrecision] will use Mathematica's arithmetic to give a result to about 16 digits of precision, the value of the numerical constant $MachinePrecision. This is different from N[expr, MachinePrecision]. MachinePrecision is a symbol and $MachinePrecision is a number.

Unlike most other similar programs, Mathematica is able to keep track of the precision of inexact numbers. Every inexact number has an attached property which indicates its precision:

In[17]:= Precision[1.0]
Out[17]= MachinePrecision

In[18]:= Precision[1.0`5]
Out[18]= 5.

If you use MachinePrecision, the precision is not tracked:

In[19]:= Precision[1.1 - 1.0]
Out[19]= MachinePrecision

If you use arbitrary precision arithmetic, then Mathematica computes an estimate of the precision of the result:

In[20]:= Precision[1.1`5 - 1.0`5]
Out[20]= 3.67778

Note the loss of precision when we subtract two numbers to obtain a result an order of magnitude smaller than both. We started with numbers specified to 5 digits and ended up with a result precise to only ~3.68 digits. Note that the precision of the result is only estimated, so don't rely on it being 100% accurate.

This feature allows Mathematica to increase the number of digits it uses for calculations to obtain as many precise digits in the result as you request:

expr = 95881665812878 - 
  120000000000000 (121576521638975 - (321097753837557 Log[
         1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/
               3))])/Log[
       5]) + (321097753837557 Log[(-1 + 
         1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/
               3)))/(-6 + 
         4/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/
               3)))])/Log[5]

As you noticed,

N[expr, 3]

(* ==> 0.487 *)

Let's do an experiment now and convert all exact numbers in the input expression to inexact numbers with only 16 digits of precision (which is what machine precision arithmetic uses):

In[22]:= expr /. x_?NumberQ :> N[x, 16]
Out[22]= 0.*10^13

This time the result seems to be "0", but notice the red highlight (if you try this in a notebook), and let's check its precision:

In[23]:= Precision[%]
Out[23]= 0.

Mathematica estimates that when doing the calculations with 16 digits only, there's no precision left in the result, so it's useless! This explains why the result obtained with machine precision arithmetic is wrong. Let's try with 50 digits now:

In[24]:= expr /. x_?NumberQ :> N[x, 50]
Out[24]= 0.48687991410140791970

In[25]:= Precision[%]
Out[25]= 20.3909

Now we got a result that's accurate to 20 digits.

When you request N[expr, 3], the number of digits the system uses for calculations is automatically increased so the result will be accurate to 3 digits.

---

If you come from, or switch to, a different system that can do arbitrary precision arithmetic, it's good to remember that in other systems the number of digits typically indicated the working precision, not the precision of the requested result. Here's an example with Maple:

expr := -14589182596676904118334187122+38531730460506840000000000000*ln((1/3)*(19+3*sqrt(33))^(1/3)+(1/3)*(19-3*sqrt(33))^(1/3)+1/3)/ln(5)+321097753837557*ln(((1/3)*(19+3*sqrt(33))^(1/3)+(1/3)*(19-3*sqrt(33))^(1/3)-2/3)/((4/3)*(19+3*sqrt(33))^(1/3)+(4/3)*(19-3*sqrt(33))^(1/3)-14/3))/ln(5);

evalf(expr, 20);

3.34187122499671*10^8 # incorrect result computed with (insufficient) 20 digits of working precision

evalf(expr, 50);

0.486879914101407919700386926722921333

---

Further reading:

• When I can assume that all decimal digits returned by Mathematica are provably correct?

Answer 2

My guess is that the wrong answer is the one given by N[expr] and not N[expr,3]. My mind is simple and I cannot manage big numbers, so, let's give 'em names:

aN = 95881665812878; bN = 120000000000000; cN = 121576521638975; 
dN = 321097753837557;
eN = Log[1/3 (1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3))];
fN = 321097753837557;
gN = Log[(-1 + 1/3 (1 + (
            19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3)))/(-6 + 4/3 (
                  1 + (19 - 3 Sqrt[33])^(1/3) + (19 + 3 Sqrt[33])^(1/3)))];

Let's define a substitution rule that will allow us to use the simple names a ... g

exactRule = {a -> aN, b -> bN, c -> cN, d -> dN, e -> eN, f -> fN, g -> gN};

Now, your expression is

expr = a - b (c - (d e)/Log[5]) + (f g)/Log[5]

where b c and b d e/Log[5] are veeeeeery similar one to the other, and this give rise to a cancellation error.

Please note that you can rewrite your expression as

expr1 = a + (b d e + f g - b c Log[5])/Log[5]

This expression seems to be less prone to suffer from cancellation errors, probably because of the different order in which the operations are performed.

Now, the real question is: why is N[expr] returning a wrong answer, when N[expr,3] is giving a correct approximation? Perhaps N[expr, nn] is computed in high precision and then in nn precision, while N[expr] simply forces machine precision, which gets fooled by cancellation errors.

EDIT^2: Yes, I believe it's a cancellation error problem. I tried on paper three forms for the above expression and there always is a double cancellation between two quantities*10^14 and two quantities*10^28. If you substitute each value of a, b, c,..., g with the corresponding N[a,nn], you can see that you need to increase nn to get consistent answers in all cases. But if you compute the expression with exact values and then apply N[expr,nn] you get the correct approximation with nn precision.

EDIT^3 : To exemplify that, let's use the following accessory procedures (with a caveat: the special case with parameter 0 is introduced to simplify the code for producing the table of values, but it does not produce the effect of N[expr,0]):

argN[expr_, rule_, n_] := expr /. (N[#, n] & /@ rule)
argN[expr_, rule_, 0] := expr /. (N[#] & /@ rule)

exprN[expr_, rule_, 0] := N[expr /. rule]
exprN[expr_, rule_, n_] := N[expr /. rule, n]

argN[expr, rule, n] evaluates expression expr by substituting the parameters specified by the substitution rule with precision n, while exprN[expr, rule, n] returns the numeric value of expr with precision n after it has been computed according to the substitution rule (the special case n=0 means it will use N[] without parameters).

Some code to create the correspondent text is here:

argNtext[n_] := If[n == 0, "expr[N[vars]]", StringJoin["expr[N[vars,", ToString[n], "]]"] ]
exprNtext[n_] := If[n == 0, "N[expr]", StringJoin["N[expr,", ToString[n], "]"] ]

And now, what follows will show how different the result is between the different cases, that is expr[N[vars]], N[expr[vars]] and expr[N[vars,n]], N[expr[vars],n] for different values of n when the parameters are replaced with the exact numbers specified by exactRule.

Table[{argNtext[k], argN[expr, exactRule, k],
      exprNtext[k], exprN[expr, exactRule, k]},
    {k, 0, 42, 3}] // TableForm

For Mathematica novices, try to figure out what you will get before evaluating the above table.

Answer 3

Here's a simpler example of a similar phenomenon where it might be easier to see what is happening. Consider:

10^100 - (10^100 - 1)
1

10.^100 - (10.^100 - 1.)
0.

The first of these operates on exact numbers and hence gives the answer as an exact number. The second (because of the 10.) operates in floating point numbers which have limited accuracy (in this case, less than one part in 10^100). The first case is like your your first, where all the numbers are exact (your answer is then turned into floating point at the end of the calculation by the N). The second (using N[stuff,3]) carries out all the calculations in limited precision arithmetic.