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

Mathematica works with exact numbers and with two different types of approximate numbers: machine-precision numbers that take advantage of specialized hardware for fast arithmetic on your computer, and arbitrary-precision numbers that are correct to a specified number of digits.

To be sure of n correct digits, use N[expr, n].

When you do a computation, Mathematica keeps track of which digits in your result could be affected by unknown digits in your input. It sets the precision of your result so that no affected digits are ever included. This procedure ensures that all digits returned by Mathematica are correct, whatever the values of the unknown digits may be.

Mathematica automatically increases the precision that it uses internally in order to get the correct answer

Of course, this sounds very reassuring, but I still have some doubts that all decimal digits ever returned by Mathematica when working with arbitrary-precision numbers are always provably correct, no matter what functions I invoked.

What are those cases when I can be certainly sure all displayed digits are correct?

Update:
Here is an example when some incorrect decimal digits are returned when working with arbitrary-precision arithmetic:

a = 17
(* 1.000000 *)

a // Precision
(* 7. *)

d = Derivative[0, 1][StieltjesGamma][0, a]
(* -1.6450 *)

MachineNumberQ[d]
(* False *)

d // FullForm
(* 1.645015523910436949472512823780090832695.155856939311388 *)

d // Precision
(* 5.15586 *)


So, Mathematica claims that at least 5 (hence, all) decimal digits of the result -1.6450 are correct. But in fact, the exact result is -Pi^2 / 6 that is -1.644934..., so only 3 digit are correct.

I am also concerned that Precision[...] itself returns a machine-precision number, that is subject to uncontrolled error-accumulation that possibly can result in claiming more digits of precision in a number than there actually are. Can I assume that Mathematica always errs on the safe side when computing a precision?

Update 2:
Another (gross) example:

a = 26
(* 2.00000 *)

Derivative[0, 1][StieltjesGamma][0, a]
(* 0.324 *)

d // FullForm
(* 0.323995226098963375800273854568809784893.339102855094484 *)

Precision[d]
(* 3.3391 *)


Here, one would expect that at least 0.32 are correct digits. But in fact, the exact result is 1 - Pi^2/6 that is -0.644934.... No correct digits, even the sign is wrong.

• The authors of the book The SIAM 100-Digit Challenge gave Mathematica (v5) solutions to the 10 problems as well as other implementations. Code here. Jun 24 '13 at 16:09
• You'll probably find Oleksandr Pavlyks screencast on Mathematical Numerics and Special Functions interesting.
– ssch
Jun 24 '13 at 16:31
• +1 toward your first "Good question" badge. Jun 24 '13 at 20:14
• BTW, I just realized I could write StieltjesGamma[0, #]&'[a] instead of a more verbose Derivative[0, 1][StieltjesGamma][0, a]. Jun 24 '13 at 23:36
• Probably, you should never expect that NIntegrate returns only provably correct digits. Jun 27 '13 at 7:09

Control the Precision and Accuracy of Numerical Results

This is an excellent question.

Of course everyone could claim highest accuracy for their product.

To deal with this situation there exist benchmarks to test for accuracy.

One such benchmark is from NIST. This specific benchmark deals with the accuracy of statistical software for instance.

The NIST StRD benchmark provides reference datasets with certified computational results that enable the objective evaluation of statistical software.

In an old issue of The Mathematica Journal, Marc Nerlove writes elaborately about performing the linear and nonlinear regressions using the NIST StRD benchmark (and Kernel developer Darren Glosemeyer from WRI discussing results using Mathematica version 5.1).

Numerically unstable functions:

But this is only one part of story. OK. There exist benchmark for statistical software etc., but what happens if we take some functions that are numerically unstable?

Stan Wagon has several examples of inaccuracies and how to deal with them in his book Mathematica in Action, which I can only warmly suggest. I have it now for (the latest edition) several years and everytime there is something new to discover with Mr. Wagon.

Let's take, for instance a numerical unstable Maclaurin polynomial of $\sin x$:

poly = Normal[Series[Sin[x], {x, 0, 200}]];
Plot[poly, {x, 0, 100}, PlotRange -> {-2, 2},
PlotStyle -> {Thickness[0.0010], Black}]


The result this we can see that the result breaks down at ~40: If we take one value x = 60 and perform a division we get a result back:

N[poly /. x -> 60] ==> -0.304811


Inserting the approximate real number 60.; there occurs a roundoff error:

poly /. x -> 60. ==> -4.01357*10^9


But inserting the number 60 (without the period); there is no problem at all:

ply /. x -> 60 ==> -((3529536438455<<209>>9107277890060)/(1157944045943<<210>>4588491415899))


The use of machine precision (caused by the decimal point) leads to an error:

10^17 + 1./100 - 10^17 ==> 0.


Machine precision is $53 \log_{10}(2) = 15.9546$.

This is the exact moment where N comes into play. We have to increase the precision:

poly /. x -> N[60,20] ==> 0. x 10^7


Still not good enough, because this number has no precision at all. So, let's increase the precision again:

poly /. x -> N[60,200] ==> -0.9524129804151562926894023114775409691611879636573830381666715331536022870514582375567159979758451142049758239018693823215314740415313661058559273332324475257579234995809519


This looks much better. If we impose the precision in our prior plot:

Plot[poly, {x, 0, 100}, PlotRange -> {-2, 2},
PlotStyle -> {Thickness[0.0010], Black}, WorkingPrecision -> 200] Not ideal, since in order to get an accurate result, we need to know what precision we need. There are numerical results which tend to lose precision during several iterations. Luckily there is some salvation in form of the Lyapunov exponent (denoted $\lambda$), which can quantify the loss of precision.

Conclusion:

What I've learned from here is, that it is a bad idea to mix small numbers with big ones in a machine precision environment. This is where Mathematica's adaptive precision comes into play.

Mathematica precision handling

Let's investigate further about precision handling inside Mathematica.

If we want to calculate $\sin(10^{30})$ in Mathematica we get:

N[Sin[10^30]] ==> 0.00933147


Using WolframAlpha we get:

WolframAlpha["Sine(10^30", {{"DecimalApproximation", 1}, "Content"}] ==>  - 0.09011690191213805803038642895298733027439633299304...


The result we get from our numerical workhorse is simply the wrong answer and this is getting worse if we increase the exponent.

(The guys at WolframAlpha seem to do it somewhat differently...but what?)

If we take $10^{30}$ and put turn this into a software real with $MachinePrecision as the actual precision we get 0 as the result, with the precision 0. This result is useless. Luckily we do know that it is indeed. Here the adaptive precision comes into play. The adaptive precision is controlled through the system variable $MaxExtraPrecision (default value is 50).

Let's say we want to compute $\sin(10^{30})$ but with a precision of 20 digits:

N[Sin[10^30], 20] ==> -0.090116901912138058030


Ah! We're getting close to the WolframAlpha engine!

If we ask for $\sin(10^{60})$ the result is:

N[Sin[10^60], 20] ==> N::meprec: Internal precision limit
$MaxExtraPrecision = 50. reached while evaluation Sin. >> Out= 0.8303897652  We run into problems, since the adaptive algorithm only adds 50 digits for extra precsion. But, luckily, the extra precision is controlled through $MaxExtraPrecision, which we're allowed to change:

$MaxExtraPrecision = 200; N[Sin[10^60], 20] ==> 0.83038976521934266466  Addendum (Michael E2): Note that N[Sin[10^30]] does all the computation in MachinePrecision without keeping track of precision; however N[Sin[10^30], n] does keep track and will give an accurate answer to precision n. (WolframAlpha probably uses something like n = 50.) Also specifying the precision of the input to be, say, 100 digits,N[Sin[10^60100], 20] will use 100-digit precision calculations internally and return the same answer as above to 20 digits of precision, provided as in this case 100 digits is enough to give 20. (Added at the request of @stefan.) Conclusion Equipped with that knowledge we could define functions that use adaptive precision to get an accurate result. Precision and accuracy It is not that Mathematica loses precision, but in your definition of a you'll lose precision in the first place. Let's first talk about precision and accuracy. Basically the mathematical definition of precision and accuracy is as follows: Suppose representation of a number$x$has an error of size$\epsilon$. Then the accuracy of$x \pm \epsilon/2$is defined to be$-\log_{10}|\epsilon|$and its precision$-\log_{10}|\epsilon/x|$. With these definitions we can say that a number$z$with accuracy$a$and precision$p$will lie with certainty in the interval: $$\left(x-\frac{10^{-a}}{2},\frac{10^{-a}}{2}\right)=\left(x-\frac{10^{-p} x}{2},\frac{10^{-p} x}{2}+x\right)$$ According to these definitions the following relation holds between precision and accuracy:$\operatorname{precision}(x)=\operatorname{accuracy}(x)+\log_{10}(|x|),$where the latter is called the scale of the number$x$. We can check if this identity holds: Function[x, {Precision[x], Accuracy[x] + Log[10, Abs[x]]}] /@ {N[1, 100], N[10^100, 30]} ==> {{100.,100.},{30.,30.}} (* qed *)  Let's define a function for both precision and accuracy: PA[x_] := {Precision[x], Accuracy[x]}  Now let's look at your definition of a: a = 17 PA[a] ==> {7., 7.} d = Derivative[0, 1][StieltjesGamma][0, a] ==> -1.6450 PA[d] ==> {5.15586, 4.93969}  You've lost precision! You defined a to have a precision and an accuracy of 7. But what is the precision and accuracy if you turn a into a symbol using machine precision: a = 1. PA[a] ==> {MachinePrecision, 15.9546}  This is a gain in precision obviously. Now let's call your canonical examples: d = Derivative[0, 1][StieltjesGamma][0, a] ==> -1.64493  Which is the exact result of$-\frac{\pi ^2}{6}$. The precision and accuracy of d is: PA[d] ==> {MachinePrecision, 15.7384}  Perfect. Now let's redefine your a to be 2. instead of 26: a = 2. PA[a] ==> {MachinePrecision, 15.6536} d = Derivative[0, 1][StieltjesGamma][0, a] ==> -0.644934  Which is the exact result of$1 - \frac{\pi ^2}{6}$PA[d] ==> {MachinePrecision, 16.1451}  Conclusion Dealing with numerical computing is dealing with loss of precision. It seems that Mathematica varies the Precision depending on the numerical operation being performed and the Precisions are more pessimistic than optimistic, which is actually quite good. In most calculations, one typically loses precision, but with an appropriate starting value you can gain precision as well. The general rule for the usage of high-precision numbers is: If you want to gain high precision you need to use high-precision numbers in your expression to be calculated. Consequently, every time you need a high-precision result you must take care that the starting expression has sufficient precision. There exists an exception to the above rule. If you use machine-precision arithmetic in expressions and the numbers are getting bigger than $MaxMachineNumber, Mathematica will switch automatically to high-precision numbers. If this is the case the rules apply as described in my Edit 2.

P.S.:

This was one of the questions I really like, since I know now more about that topic than before. Maybe one of the WRI/SE jedi can join the party to provide even more insights on this, possibly more than I would ever be able to provide.

• @MichaelE2 Thank you for your helpful input. I'd like to invite you, if you can afford the time, to update the post and enter your findings with your attribution. If this is ok for you. Or post it in a new post, if this is to circuitous for you. But thanks again for your input. Jun 24 '13 at 15:14
• Thanks. I added it, with correct formatting :), as an add-on to "Edit 2". You may wish to roll back, or re-edit it yourself, if you would like to incorporate it more effectively. That would be fine with me. (I'll delete my initial comment soon, I think.) Jun 24 '13 at 15:54
• @MichaelE2 nice job :) i appreciate your input indeed. Jun 24 '13 at 19:42
• @Stefan Thank you for your deep analysis of the subject! The bounty is yours. Jul 2 '13 at 18:32
• @VladimirReshetnikov thank you for accepting it :) I was about to write an addendum on NumericalMath$NumberBits function. This function shows how Mathematica simulates interval arithmetic by constantly maintaining a few more digits than needed... Jul 2 '13 at 18:36 This isn't an answer (yet) but it was too long for a comment. Here is an extended example where the quoted precision does not appear to be true: f1 = Derivative[0, 6][StieltjesGamma][0, #] &; f2 = {Accuracy@#, Precision@#, InputForm@#} &; f1 @ 115 f2 @ % N[f1 @ 1, 10] f2 @ %  725.59 {2.29961, 5.1603, 725.592784177191488027966503357541438517565.1603018938964} 726.0114797 {7.13906, 10., 726.0114797147721516339424224809075397033310.}  Update As @MichaelE2 points out, perhaps I didn't directly answer your question. I think the short answer is that you can rely on precision of the output if you are using adaptive precision control, and you can't rely on the precision of the output if you are not using adaptive precision control. In your examples you are not using adaptive precision control, so your results are not surprising. Extended response If you have a function, and you insert an inexact number, then Mathematica has very little control over the precision/accuracy of the output. Here is an extremely simple example: Precision[1.00015 - 1]  0.999957 You started off with a precision of 5, and then Mathematica returned a result with a precision of 1. The loss of precision is due to subtractive cancellation. If instead you insert an exact number and numericize, then Mathematica can use adaptive precision control to insure that the output has the precision requested during the numericization process. Adaptive precision control basically means that while numericizing an input, exact numbers can be numericized to much higher precision when needed, controlled by the global variable $MaxExtraPrecision. More schematically:

Precision[f[N[exact, prec]]] ->  ?
Precision[N[f[exact], prec]] -> prec


So, in your question, you are doing the former, and hence you are circumventing Mathematica's adaptive precision control, and so, just like with MachinePrecision computations, the result may not have any correct digits. Here's another simple example. Suppose your function is:

f[x_] := 1/1000 - Sin[x];


If we feed f an inexact number close to 1/1000, then due to subtractive cancellation, the output will have much less precision than the input:

Precision @ f[N[1/1000, 10]]


3.22185

On the other hand, if you feed an exact number, and then numericize, the evaluation will use adaptive precision control:

Precision @ N[f[1/1000], 10]


10.

We can use the following trick to see adaptive precision control in action:

Clear[f]
f[x_Real] := 1/1000 - Sin[x]

TracePrint[
N[f[1/1000], 10],
_f,
TraceInternal->True,
TraceAction->(Print[InputForm[#]]&)
]
Precision[%]


HoldForm[f[1/1000]]

HoldForm[f[1/1000]]

HoldForm[f[0.00110.]]

HoldForm[f[0.00129.93398985691138]]

1.666666583*10^-10

10.

You can see that N first feeds f an inexact number with precision 10, and as we saw earlier, this is not sufficient to produce an answer with the requested precision. Then, N feeds f an inexact number with precision 30, and this was sufficient.

The moral of the story is N[f[exact], prec] should produce a result with precision prec, while f[N[exact, prec]] need not.

• Concerning calculations like your first example, isn't the OP's question, when can we be sure the first digit is correct if the Precision is at least 1? In the OP's examples, the Precision is greater than the actual precision of the result. My read of the OP is that the question is about the accuracy of Precision`, not about how to achieve a certain precision. Feb 8 '18 at 18:22