When I can 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. –  Michael E2 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. –  Mr.Wizard 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]. –  Vladimir Reshetnikov Jun 24 '13 at 23:36
Probably, you should never expect that NIntegrate returns only provably correct digits. –  TauMu Jun 27 '13 at 7:09
show 1 more comment

Control the Precision and Accuracy of Numerical Results

This is an excellent question.

Of course everyone could claim highest accuracy for her 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 computional 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 loose 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[1000000000000000000000000000000000000000000000000000000000000]. >> Out[105]= 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 looses precision, but in your defintion of a you'll loose 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:$(x-\frac{10^{-a}}{2},\frac{10^{-a}}{2})=(x-\frac{10^{-p} x}{2},\frac{10^{-p} x}{2}+x)$According to these definitions the following relation holds between precision and accuracy:$precision(x)=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 looses 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 high-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 yedi's join the party to give even more insights on that matter, than I would ever been able to.

-
@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. –  Stefan 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.) –  Michael E2 Jun 24 '13 at 15:54
@MichaelE2 nice job :) i appreciate your input indeed. –  Stefan Jun 24 '13 at 19:42
@Stefan Thank you for your deep analysis of the subject! The bounty is yours. –  Vladimir Reshetnikov 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... – Stefan Jul 2 '13 at 18:36 show 1 more comment 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.}  - add comment In your example when you enter a = 17 Derivative[0, 1][StieltjesGamma][0, a]  the derivative is being computed numerically. If you compute the derivative analytically Clear[a]; D[StieltjesGamma[0, a], a]  and then substitute in the value of$a\$

% /. a -> 17


no incorrect decimal digits are returned.

-
The whole question is about arbitrary-precision numeric computations. –  Vladimir Reshetnikov Jun 26 '13 at 3:00
@TheDoctor It appears your answer is not for the question being asked. Please consider changing it or removing it. –  R Hall Jun 30 '13 at 16:07