Machine precision near zero: not fulfilled?

Question

I am puzzled by the behavior of Mathematica machine precision with numbers approaching zero. This manifests itself, e.g., with FixedPoint and the like.

In the examples below I will use the following "service" function to take only a few values of long lists:

ClearAll[take];
take[n_Integer][li_List] := If[Length[li] >= Abs[n], Take[li, n], li];

Consider finding the solution of two different equations:

FindRoot[Sin[2.0 x] == x, {x, 1}]

FindRoot[Sin[0.5 x] == x, {x, 1}]

They return respectively:

{x -> 0.9477471335169905}

{x -> 2.6312502524050707*^-15}

where the second one is actually zero to machine precision. edit: no: where the second can be considered as zero in many machine precision computations, even if it is several orders of magnitude bigger than the smallest floating point number.

Now consider doing the same with FixedPoint:

a = FixedPoint[Sin[2.0 #] &, 1.]

b = FixedPoint[Sin[0.5 #] &, 1.]

The behavior of the functions is such that the two should converge to the respective solutions above. Actually both never end.

An analysis of the values obtained (with FixedPointList) shows that the problem with the first function is an oscillation in the last bit(s):

take[-5]@ FixedPointList[Sin[2.0 #] &, 1., 1000] // Differences

{-2.220446*^-16, 2.220446*^-16, -2.220446*^-16, 2.220446*^-16}

so that a real fixed point is never achieved. In this case, SameTest -> Equal is of help:

FixedPoint[Sin[2.0 #] &, 1., SameTest -> Equal]

0.9477471335169858

On the other hand, using SameTest -> Equal does not help with the other example, since:

FixedPoint[Sin[0.5 #] &, 1., SameTest -> Equal]

also never ends.

The same analysis as before, though, shows this:

take[-5]@ FixedPointList[Sin[0.5 #] &, 1., 10000] // Differences

{-3.7962812653512530891065876917712356`15.477468515471298*^-3010,

-1.8981406326756265445532938458856178`15.477468515471298*^-3010,

-9.490703163378132722766469229428089`15.477468515471298*^-3011,

-4.745351581689066361383234614714044`15.477468515471298*^-3011}

Actually there is no bottom! The negative exponents become arbitrarily large and two successive numbers are never seen to be equal to zero.

Now the question: is it reasonable that, having given as input MachinePrecision numbers, arbitrarily large exponents are not confined into the expected range? This is what Chop would do, but it is not applied and must be forced with a user-given SameTest, which seems a little unexpected here.

Why can MachinePrecision numbers have an arbitrary negative exponent?

Answer 1

M11.3 has changed the way machine number underflows are handled. Previously, when a function was given a machine number input and produced a result that was so small that it could not be represented by a machine number (machine underflow), then Mathematica switched to using extended precision numbers. In M11.3, these small numbers now just evaluate to a machine number 0. (with a message):

Exp[-900.]

General::munfl: Exp[-900.] is too small to represent as a normalized machine number; precision may be lost.

0.

You can choose to use extended precision inputs instead, in which case an extended precision result will be returned:

Exp[-900`20]

1.3644772123656828*10^-391

It is not possible to change M11.3 so that it handles machine underflow in the same way that M11.2 does.

At any rate, this means that the output of M11.3 for your question has changed:

FixedPoint[Sin[.5 #] &, 1.]

General::munfl: 0.5 4.21308*10^-308 is too small to represent as a normalized machine number; precision may be lost.

0.

Answer 2

According to $MinMachineNumber (version <= 11.2), MachinePrecision computations are automatically converted to arbitrary precision numbers (of precision $MachinePrecision) when a computation results in a value less than $MinMachineNumber.

Precision[$MinMachineNumber]
Precision[$MinMachineNumber/2.]
(*
  MachinePrecision
  15.9546
*)

Update for 11.3+:

Now according to $MinMachineNumber, "machine numbers smaller than \$MinMachineNumber are represented as subnormal machine numbers." One also can get underflow if a result is too small to represent as a subnormal floating-point number.

Answer 3

Based on the reference to $MinMachineNumber from Michael E2 I found the practical solution to the FixedPoint never-ending search:

SetSystemOptions["CatchMachineUnderflow" -> False];

FixedPoint[Sin[0.5 #] &, 1.]

0.

SetSystemOptions["CatchMachineUnderflow" -> True]; (* the default *)

FixedPoint[Sin[0.5 #] &, 1.] (* will never stop *)

Consequently I will start using this option that I have been never aware of.

It was a little of a trap to me the fact that you may start with floating point numbers and be carried away from them because of the fact that:

"The Wolfram Language uses arbitrary-precision to represent numbers that would be denormalized".

So the reply to my question "is it reasonable that, having give as input MachinePrecision numbers, arbitrarily large exponents are not confined into the expected range?" is that this behavior can be controlled with the SystemOption above so that it can become reasonable whatever point of view one may have.

EDIT: The fact that zero is treated in a very special way by Mathematica (and I leave to the experts of numerical analysis to discute if this is good or not) can be highlighted by the following cases:

SetSystemOptions["CatchMachineUnderflow" -> True] (* default value *)

CatchMachineUnderflow->True

FixedPoint[Sin[#/2] &, 1.] // InputForm

$Aborted

FixedPoint[Sin[#/2] + 10.^-308 &, 1.] // InputForm

2.00000000000000140941884798286`15.954589770191005*^-308

FixedPoint[Sin[#/2] + 10.^-308 &, 1., SameTest -> Equal] // InputForm

2.00000000000001614340458780004`15.954589770191005*^-308

SetSystemOptions["CatchMachineUnderflow" -> False];

FixedPoint[Sin[#/2] &, 1.] // InputForm

0.

FixedPoint[Sin[#/2] + 10.^-308 &, 1.] // InputForm

2.*^-308

FixedPoint[Sin[#/2] + 10.^-308 &, 1., SameTest -> Equal] // InputForm

2.0000000000000146*^-308

When the fixed point value of the function is even slightly shifted from zero (+10.^-308) MMA finds the solution immediately. As soon as the value is zero, the default action no longer works, and one must be much more careful not to become stuck.