**NumericalMath****`****ComputerArithmetic****`**

The arithmetic used by Mathematica is a mixture of variable precision software arithmetic and whatever is provided by the manufacturer of the floating-point hardware (or the designer of the compiler, if there is no floating-point hardware). If you want to examine the machine arithmetic on your machine you can use the functions defined in NumericalMath`Microscope`. If you want to learn about the basic ideas of computer floating-point arithmetic in general you can use NumericalMath`ComputerArithmetic`. This allows you to examine arithmetic with various bases, precisions, and rounding rules.

Computer numbers and nonnumbers in NumericalMath`ComputerArithmetic`.

Much of the information carried around in the data object that makes up a computer number is redundant. In particular, the first three arguments contain exactly the same information as the fourth argument. The redundancy exists partly for the sake of efficiency and partly to allow the user access to the various fields. The fifth argument has nothing to do with the computer number itself. It instead represents what the value of the number would be without the cumulative effects of all the roundoff errors that went into the computer number. It is computed using Mathematica's high-precision arithmetic and can generally be regarded as the "correct" value of the number. Comparing the computer number with this number gives the error in the computer number.

This loads the package.
In[1]:= **<<NumericalMath`ComputerArithmetic`**

Here is the computer number representing the ordinary number 2.
In[2]:= **a = ComputerNumber[2]**

Out[2]=

This gives the computer number representation of

.
In[3]:= **b = ComputerNumber[Pi]**

Out[3]=

Arithmetic works with computer numbers just as with ordinary numbers.
In[4]:= **c = a + b**

Out[4]=

Here is the complete structure of c.
In[5]:= **InputForm[c]**

Out[5]//InputForm=

You can also enter just the sign, mantissa and exponent.
In[6]:= **ComputerNumber[-1, 1234, -6]**

Out[6]=

But if your input doesn't make sense in the arithmetic currently in effect you get NaN. Here the problem is that the mantissa is only three digits long.
In[7]:= **ComputerNumber[-1, 123, 7]**

Out[7]=

Changing the type of arithmetic to be used.

The default arithmetic is four digits in base 10 with a rounding rule of RoundToEven. Only numbers between and

are allowed. Mixed-mode arithmetic is not allowed and division is not done with correct rounding. (It is performed as two operations: multiplication by the reciprocal. Each operation involves rounding errors.)
In[8]:= **Arithmetic[ ]**

Out[8]=

Now the arithmetic is set to be six digits in base 8.
In[9]:= **SetArithmetic[6, 8]**

Out[9]=

The result is displayed in octal (with trailing zeros suppressed).
In[10]:= **ComputerNumber[Pi]**

Out[10]=

There are several options that can be used when setting the arithmetic. In addition to changing the precision and the base, you can control the type of rounding used, the magnitude of the numbers allowed, whether mixed-mode arithmetic is to be allowed, and whether division is to be done with a single rounding operation.

Options for SetArithmetic.

It should be noted that correctly rounded division is implemented by the function IdealDivide. This function can be used whether or not the arithmetic is set to automatically use correct division. It is difficult to get x/y to give correctly rounded division since the x/y is converted to xy^(-1) before the package ever sees it. The way we get around this is to define $PreRead to convert x/y into x~IdealDivide~y before the parser has a chance to change it to xy^(-1). If you want to use $PreRead for your own purposes this will interfere with your definition.

Types of rounding available.

Mixed-mode arithmetic (any arithmetic operation involving an integer and a computer number) is not allowed.
In[11]:= **3 ComputerNumber[Pi]**

Out[11]=

This turns on mixed-mode arithmetic (and sets the arithmetic to six digits in hexadecimal).
In[12]:= **SetArithmetic[6, 16, MixedMode -> True]**

Out[12]=

Now the product is computed.
In[13]:= **3 ComputerNumber[Pi]**

Out[13]=

There are many things about computer arithmetic that are very different from ordinary arithmetic. Most of these differences are quite easy to demonstrate.

This sets the arithmetic back to the default.
In[14]:= **SetArithmetic[4, 10]**

Out[14]=

Expressions are evaluated numerically before they become computer numbers.
In[15]:= **ComputerNumber[Pi - 22/7]**

Out[15]=

If the separate parts are converted to computer numbers first, catastrophic cancellation of digits can occur.
In[16]:= **ComputerNumber[Pi] - ComputerNumber[22/7]**

Out[16]=

Summing the terms from smallest to largest gives one answer.
In[17]:= **(sum = 0;**

Do[sum += ComputerNumber[i]^(-2),

{i, 200}];

FullForm[sum])

Out[17]//FullForm=

As a general rule it is better to sum from smallest to largest. You can see what the error is by comparing the mantissa (the second argument) to the "correct" value (the last argument).
In[18]:= **(sum = 0;**

Do[sum += ComputerNumber[i]^(-2),

{i, 200, 1, -1}];

FullForm[sum])

Out[18]//FullForm=

Here is an example where summing from largest to smallest gives the better result.
In[19]:= **(sum = 0;**

Do[sum += 1/ComputerNumber[i], {i, 300}];

FullForm[sum])

Out[19]//FullForm=

The difference is slight, and such examples are rare.
In[20]:= **(sum = 0;**

Do[sum += 1/ComputerNumber[i],

{i, 300, 1, -1}];

FullForm[sum])

Out[20]//FullForm=

Basic arithmetic is all that is implemented in the package. We could easily extend things to include elementary functions.
In[21]:= **Sin[ComputerNumber[N[Pi]/7]]**

Out[21]=

Here is the square root of 47.
In[22]:= **sq = ComputerNumber[Sqrt[47]]**

Out[22]=

It is a theorem that correctly rounded square roots of small integers will always square back to the original integer if the arithmetic is correct.
In[23]:= **sq sq**

Out[23]=

But a similar theorem for cube roots does not exist.
In[24]:= **cr = ComputerNumber[3^(1/3)]; cr cr cr**

Out[24]=

This changes the arithmetic to seven digits in base 10 with a rounding rule of RoundToEven and an exponent range of

50 to 50.
In[25]:= **SetArithmetic[7]**

Out[25]=

This number rounds down.
In[26]:= **ComputerNumber[.9999999499999999999999999]**

Out[26]=

This number rounds up because it rounds toward the mantissa that is even (1000000) rather than the one that is odd (9999999).
In[27]:= **ComputerNumber[.9999999500000000000000000]**

Out[27]=

Again it rounds toward the even mantissa (1000000 rather than 1000001).
In[28]:= **ComputerNumber[1.000000500000000000000000]**

Out[28]=

The reciprocal of the reciprocal is not the original number; in fact it may be quite different.
In[29]:= **(x = ComputerNumber[9010004];**

y = 1/x;

z = 1/y)

Out[29]=

This is multiplication by the reciprocal. It involves two rounding operations.
In[30]:= **ComputerNumber[2]/x**

Out[30]=

This is true division. It uses a single rounding operation. You can get this better type of division to work with "/" by setting the option IdealDivide to True when you use SetArithmetic.
In[31]:= **IdealDivide[ComputerNumber[2], x]**

Out[31]=

Normal converts a computer number back to an ordinary rational number with exactly the same value.
In[32]:= **Normal[ComputerNumber[Pi]]**

Out[32]=

Now we change the arithmetic again.
In[33]:= **SetArithmetic[3, 2,**

RoundingRule -> Truncation,

ExponentRange -> {-3,3}]

Out[33]=

It is easy to plot the error in the computer number representation of each number.
In[34]:= **Plot[Normal[ComputerNumber[x]] - x,**

{x, -10, 10}, PlotPoints -> 193]

ComputerNumber::ovrflw: Overflow occurred in computation. The exponent is 4.

ComputerNumber::ovrflw: Overflow occurred in computation. The exponent is 4.

ComputerNumber::ovrflw: Overflow occurred in computation. The exponent is 4.

General::stop: Further output of ComputerNumber::ovrflw will be suppressed during this calculation.

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -10..

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -9.89858.

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -9.78798.

General::stop: Further output of Plot::plnr will be suppressed during this calculation.

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

General::stop: Further output of ComputerNumber::undflw will be suppressed during this calculation.

You can zoom in to see the hole at zero.
In[35]:= **Plot[Normal[ComputerNumber[x]] - x,**

{x, -1, 1}, PlotPoints -> 47]

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -0.122946.

ComputerNumber::undflw: Underflow occurred in computation. The exponent is -3.

General::stop: Further output of ComputerNumber::undflw will be suppressed during this calculation.

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -0.124379.

Plot::plnr: Normal[ComputerNumber[x]] - x is not a machine-size real number at x = -0.0814126.

General::stop: Further output of Plot::plnr will be suppressed during this calculation.