Arbitrary-Precision Arithmetic behaves unexpectedly

8

2

I want to understand how small numerical errors propagate through a computation and may drastically change the final result. To this end, I consider a toy world in which every number is limited to exactly 5 significant digits.

The output of

$MaxPrecision = 5;
$MinPrecision = 5;

x = SetPrecision[10^(-10), 5];
x // FullForm
y = 1 + x;
y // FullForm
y - 1 // FullForm

is

1.`5.*^-10

1.`5.

1.`5.*^-10

The last line puzzles me. The FullForm of y shows that the Precision of 5 significant digits is not sufficient to resolve the small deviation x from 1. However, y - 1 still yields the exact result, instead of y - 1 = 0.

If one writes explicitly

1.`5. - 1 // FullForm

this yields, of course,

0

as it should.

Does this mean that FullForm[y] hides something?

user19095

Posted 2017-02-02T23:01:19.597

Reputation: 369

2I think you need NumericalMath\$NumberBits[y]in order to see the guard bits that I guess are otherwise hidden in e.g.FullFormandInputForm`. Maybe there is a cleaner way, I'm not sure. – Daniel Lichtblau – 2017-02-02T23:40:26.230

@Daniel Why not post an answer? – Mr.Wizard – 2017-02-02T23:52:14.847

@Mr.Wizard It's at best an answer/2 ... – Daniel Lichtblau – 2017-02-03T00:08:15.927

You might want to look at the Computer Arithmetic package, if you're interested in this sort of thing. – J. M.'s ennui – 2017-02-03T04:22:02.637

Answers

8

This is at best a partial answer.

I don't know in full detail what is going on with FullForm or InputForm but I will venture to guess they are influenced in some way by the $XXXPrecision settings.

$MaxPrecision = 5;
$MinPrecision = 5;
InputForm[x = SetPrecision[10^(-10), 5]]
(* Out[6]//InputForm=
1.`5.*^-10 *)
InputForm[y = 1 + x]
(* Out[7]//InputForm=
1.`5. *)
InputForm[y - 1]
(* Out[8]//InputForm=
1.`5.*^-10 *)

There do seem to be guard bits retained, and they can be found using NumericalMath`$NumberBits.

NumericalMath`$NumberBits[y]

(* Out[9]= {1, {1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1,
   1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 
  1, 1}, 1} *)

NumericalMath`$NumberBits[y - 1]

(* Out[10]= {1, {1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1}, {1, 
  1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0}, -33} *)

Why are bits considered "bad" in y now regarded as "good" in y-1? I do not know.

If you make the perturbation value sufficiently small, the precision tracking will discard it.

x = SetPrecision[10^(-20), 5];
y = 1 + x;
InputForm[{x, y, y - 1}]

(* Out[45]//InputForm=
{1.`5.*^-20, 1.`5., 0} *)

Let me add two quotes from Jerry Keiper, "Mathematica Numerics: Controlling the Effects of Numerical Errors in Computation", 1992 Mathematica Developer Conference:

When SetPrecision is applied to an exact number (i.e., an integer or a rational number), a bignum is created with a few more bits than are requested and the unknown bits are so indicated. These explicit, hidden bits are in fact used in the subsequent calculations, but the boundary between the known and unknown digits is calculated for each result.

and

By changing the values of MinPrecision and MaxPrecision to be the same, say 100 you can get fixed, high-precision arithmetic. You should be aware that there is a difference between this arithmetic and machine arithmetic (other than the precision). All bignum arithmetic in Mathematica uses hidden digits, i.e., digits that are included in the calculations but are neither displayed nor considered significant. When loss of digits occurs due to cancelation of leading digits a positive value of MinPrecision may cause some of the hidden digits to become exposed as the least significant digits of the result. If there are already at least MinPrecision digits in the number no hidden digits will be exposed.

Daniel Lichtblau

Posted 2017-02-02T23:01:19.597

Reputation: 52 368

3To give credit where due: The "Let me add..." with quotes from Jerry Keiper was added by another user. I only wish I had recalled that reference. – Daniel Lichtblau – 2017-02-03T16:38:56.870