# Machine precision near zero: not fulfilled?

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.796281265351253089106587691771235615.477468515471298*^-3010,

-1.898140632675626544553293845885617815.477468515471298*^-3010,

-9.49070316337813272276646922942808915.477468515471298*^-3011,

-4.74535158168906636138323461471404415.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?

• Why do you say 2.6312502524050707*^-15 is "actually zero to machine precision"? Sep 19, 2015 at 12:33
• Sep 19, 2015 at 12:50
• @MichaelE2 I just meant that Chop[2.6312502524050707*^-15] == 0, so that it is effectively treated as zero Sep 19, 2015 at 12:51
• Chop replaces approximate real numbers closer to zero than the specified tolerance (default 10^-10) by zero. There's no connection to machine-precision zero. You might like to consider the results of these for a = 2.6312502524050707*^-15: a == 0., 1 + a == 1., 1 + a === 1., {31, 32} + a - {31, 32}. Note Equal and SameQ also have tolerances, Internal$EqualTolerance and Internal$SameQTolerance. Sep 19, 2015 at 13:05

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[-90020]


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.

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. • Yes, but FindRoot is not behaving like this and converges to a reasonable solution. FixedPoint doesn't. Sep 19, 2015 at 12:56 • @user8074 That's because of the difference in the way the two functions work. But the question didn't mention those functions but asked how a machine-number computation end up being outside the range of machine numbers. (FindRoot searches until the criterion set by PrecisionGoal and AccuracyGoal is met, provided such goals are achievable at the specified, or default, WorkingPrecision. Typically the goal is roughly to get at least 8 digits, when working at MachinePrecision.) Sep 19, 2015 at 13:17 • A "megafunction" like FindRoot can use precision control and (usually) maintain a fixed precision. Basic arithmetic operations will instead underflow to equivalent software bignums of$Machineprecision. Sep 20, 2015 at 22:32
• @DanielLichtblau by “megafunction” do you mean an extremely top level function? Apr 5, 2018 at 14:26
• @QuantumDot Basically yes, large functions that deal with options such as involving method or precision control. Also note my second sentence is no longer true as of version 11.3. Apr 5, 2018 at 15:53

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.0000000000000014094188479828615.954589770191005*^-308

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


2.0000000000000161434045878000415.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.

• Note that underflow is not the reason this runs forever: FixedPoint[0.5 - Sin[#] &, 1.]. The real solution, IMO, is to realize how FixedPoint works and programmatically include an appropriate tolerance when using it in approximate numerics. Sep 19, 2015 at 13:35