NumericalMath`Microscope`
Mathematica uses both the machine arithmetic that is provided by the machine and its own arbitraryprecision arithmetic, which is the same on all machines (except for the variations in the size of the largest representable number and the amount of memory on the machine). This package deals with the floatingpoint arithmetic provided by the machine, not the arbitraryprecision arithmetic of Mathematica.
Numbers on a computer comprise a discrete set. There are gaps between the numbers and when you do arithmetic the result often is not representable. When the result falls between two representable numbers the best that can be done is to use the closest representable number as the result instead of the correct result. The set of numbers that can be represented on a computer in floatingpoint format are commonly referred to as the set of machine numbers. If you want to investigate how arithmetic is done on computers in a general sense you should use the package NumericalMath`ComputerArithmetic`. That package allows you to vary certain fundamental parameters that control the arithmetic.
The difference between two consecutive machine numbers is called an ulp (unit in the last place, i.e., one digit in the least significant place). The size of an ulp varies depending on where you are in the set of machine numbers. Between and an ulp is equal to $MachineEpsilon. Between and it is equal to 2 $MachineEpsilon. Ideally no function should ever return a result with error exceeding half of an ulp since the distance from the true result to the nearest machine number is always less than half of an ulp, the worst case being when it is exactly halfway between two machine numbers. It is relatively easy to achieve this ideal for the four arithmetic operations and for the square root function. For the elementary transcendental functions it becomes more difficult due to the fact that they cannot be evaluated in a single operation and the combined effects of all the rounding errors would have to be made less than half an ulp. Nevertheless, it is quite possible to design algorithms for these functions such that the error is never more than one ulp.
The package NumericalMath`Microscope` provides four functions that are useful in examining the machine floatingpoint arithmetic on your machine.
Functions for examining machine floatingpoint arithmetic.
The default neighborhood of the point x = a used in Microscope and MicroscopicError is from 30 ulps to the left of a to 30 ulps to the right of a, where the ulp used is that defined at x = a, that is, Ulp[a]. Near powers of 2 and near 0 it is not clear what should be done since the size of an ulp changes. Near powers of 2 it is not much of a problem in the abscissa since an ulp is chosen to be the smaller of the two values and the resulting machine numbers just get included multiple times due to rounding effects. You can change the default size of the neighborhood by including a third value in the neighborhood specification.
A problem arises with Microscope and MicroscopicError near places where the function value becomes infinite. Since the function is changing by orders of magnitude very rapidly it is difficult to choose a good scale that displays the information in which you are interested. In such cases you may have to examine the function several hundreds of ulps away from the point you really want.
Controlling the size of the neighborhood.
The option PlotJoined controls the way the plot is drawn. It cannot be overemphasized that the functions that we are examining map the machine numbers to the machine numbers. For our purposes the real numbers in between do not even exist. Thus the plots produced by Microscope and MicroscopicError should be nothing but a set of plotted points. However, it is easier to see what is happening if the points are joined with straight line segments. Finally, if you want to see what a function does as a map from the real numbers to the real numbers, that is, rounding the real argument to the nearest machine number and applying the function to this result, you can do so, but when you see large errors, don't conclude that the function is to blame; it can't do any better than deal with the machine numbers that it gets as arguments.
Options for Microscope and MicroscopicError.
This loads the package.
In[1]:= <<NumericalMath`Microscope`
Here is the size of an ulp at 1.
In[2]:= Ulp[1]
Out[2]=
The size of an ulp changes at powers of 2.
In[3]:= {Ulp[7.999999], Ulp[8], Ulp[8.000001]}
Out[3]=
There is a large hole near 0.
In[4]:= {Ulp[0], Ulp[$MinMachineNumber]}
Out[4]=
Here are the combined errors of rounding to the nearest machine number and then taking the square root.
In[5]:= MachineError[Sqrt[x], x > 49 Pi/50]
Out[5]=
Here is the error of just the square root.
In[6]:= MachineError[Sqrt[x], x > N[49 Pi/50]]
Out[6]=
For another number close to N[49 Pi/50], the error in the square root function is less than half an ulp.
In[7]:= MachineError[Sqrt[x], x > N[49 Pi/50]  6 Ulp[49 Pi/50]]
Out[7]=
The function is not provided by the computer manufacturer so it is evaluated as a nested sequence of operations that are provided. This shows how much error is involved in the sequence, but this information is rarely useful. For testing purposes you normally want to look at an individual function evaluated at a machine number to avoid the effect of an illconditioned function magnifying prior error.
In[8]:= MachineError[Tan[Sqrt[Exp[Sin[Log[x]]]]], x > N[Pi]]
Out[8]=
Here is the function evaluated at machine numbers near 7.
In[9]:= Microscope[Log[x], {x, 7}, PlotJoined > False]
Out[9]=
Here is the function evaluated at machine numbers near 7 with the points joined by straight lines.
In[10]:= Microscope[Log[x], {x, 7, 10}]
Out[10]=
And here are the combined effects of rounding real numbers to the nearest machine number and taking the of the result. This is not a valid test of the function.
In[11]:= Microscope[Log[x], {x, 7, 10}, PlotJoined > Real]
Out[11]=
This shows the error in evaluating the sine function at machine numbers near 16. Note how the size of an ulp changes at 16. Note also that rounding has moved some points to coincide so that there appear to be fewer than points. The scale on the vertical axis represents ulps.
In[12]:= MicroscopicError[Sin[x], {x, 16, 15}, PlotJoined > False]
Out[12]=
Here is the error in evaluating sine at 21 machine numbers near 31 with the points joined by straight lines.
In[13]:= MicroscopicError[Sin[x], {x, 31, 10}, PlotJoined > True]
Out[13]=
And here are the combined effects of rounding real numbers to the nearest machine number and taking the sine of the result. This is not a valid test of the sine function. The error in evaluating the sine function at a machine number near 31 is much smaller than the combination of errors resulting from first rounding a real number to the nearest machine number and then taking the sine of that machine number.
In[14]:= MicroscopicError[Sin[x], {x, 31, 10}, PlotJoined > Real]
Out[14]=
Here you see how bad the square root function really is. It appears that near its error is bounded by about ulps.
In[15]:= MicroscopicError[Sqrt[x], {x, 3.3, 100}]
Out[15]=
You cannot examine this function exactly at 1 since it is singular there. You can, however, examine it near 1. The label for 1 + 300 Ulp[1] is printed as 1. because the label is limited to 6 characters.
In[16]:= MicroscopicError[1/(x  1)  1/(x + 1), {x, 1 + 300 Ulp[1]}]
Out[16]=
One point needs to be made clear here since rounding errors are often misunderstood and blame is incorrectly placed on an algorithm that has nothing to do with the error. Algorithms are designed to map the set of machine numbers into the set of machine numbers. To consider other numbers does not make sense. As far as the computer is concerned there are no other numbers. So, for example, if you want to find the square root of , you might do it as follows.
Start with a good numerical approximation to .
In[17]:= a = N[49 Pi/50, 30]
Out[17]=
This finds the nearest machine number to a.
In[18]:= b = N[a]
Out[18]=
This finds the square root of the machine number b using the machineprecision algorithm for square root.
In[19]:= c = Sqrt[b]
Out[19]=
Here the error in c is measured in ulps. This result gives the error for the particular machine on which this example was run. Your machine may give a very different result.
In[20]:= (SetPrecision[c,30]Sqrt[a])/$MachineEpsilon
Out[20]=
Here the error is less than half an ulp. However, this error is due to both the square root function and the rounding error incurred in converting the number a to the machine number b. In fact, the square root function has an error of less than 0.1 ulps.
Here is the correct value for the square root of b.
In[21]:= d = Sqrt[SetPrecision[b, 30]]
Out[21]=
This shows the error in the square root function.
In[22]:= (SetPrecision[c, 30]  d)/$MachineEpsilon
Out[22]=
