16
$\begingroup$

In C (and many other programming languages), there is a function

double nextafter(double x, double y)

which takes two (IEEE 754) floating-point numbers and returns the next representable floating-point number after x in the direction of y. What is the Mathematica equivalent of this function for MachinePrecision numbers?


EDIT: This issue has been brought up in the comments, but for the benefit of future readers, note that this is a substantially more subtle task than simply adding or subtracting $MachineEpsilon. The problem is that the distance between one floating-point value and the next changes with magnitude. $MachineEpsilon, by definition, is the smallest positive floating-point value such that 1.0 + $MachineEpsilon > 1.0. The distance between 1.0e-300 and the next number up will be much smaller, while the distance between 1.0e+300 and the next number up will be much greater. In addition, there are issues raised by the transitions between one order of magnitude and the next. Observe that 1.0 + $MachineEpsilon/2 == 1.0, while 1.0 - $MachineEpsilon/2 < 1.0.

$\endgroup$
  • $\begingroup$ x + Sign[y] $MachineEpsilon should work, yes? $\endgroup$ – J. M. will be back soon May 24 '15 at 4:55
  • 1
    $\begingroup$ @J. M. Not quite. $MachineEpsilon is the smallest floating-point number such that 1 + $MachineEpsilon != 1, but for certain values, the distance to the next representable float is much smaller. For example, consider the distance between 1.0e-300 and the next number up. $\endgroup$ – David Zhang May 24 '15 at 4:56
  • 1
    $\begingroup$ @J. M. If I understand correctly, incrementing/decrementing a float is a rather nontrivial task using only floating-point arithmetic operations. It's much easier with bit-level access to the internal representation of a number, which I can't figure out how to get in Mathematica. $\endgroup$ – David Zhang May 24 '15 at 5:01
  • 2
    $\begingroup$ Have a look at Ulp et al. in the Computer Arithmetic package... $\endgroup$ – ciao May 24 '15 at 5:28
  • $\begingroup$ You're right, for numbers in $(-1, 1)$, adding an appropriate signed machine epsilon would not be applicable. But at least for numbers outside that range, the last one would work. $\endgroup$ – J. M. will be back soon May 24 '15 at 5:38
7
$\begingroup$

It looks like you have to program it yourself. At a boundary x == 2.^n, the distance to the next machine real is either x * $MachineEpsilon or x * $MachineEpsilon / 2. The documentation for MantissaExponent ambiguously states that the mantissa will be "between $1/b$ and $1$". It seems be the case that $1/b \le \mathtt{mantissa} < 1$.

nextafter[0., y_] := Sign[y] $MinMachineNumber;
nextafter[x_, y_] /; Precision[x] == MachinePrecision := 
  With[{mantexp = MantissaExponent[x, 2]},
   Piecewise[{
      {First[mantexp] + Sign[y] $MachineEpsilon/4., 
       Sign[x] Sign[y] < 0 && First[mantexp] == 1./2.}},
     First[mantexp] + Sign[y] $MachineEpsilon/2.] * 2.^Last[mantexp]
   ];

A couple of tests:

Block[{x = 128. - 128*$MachineEpsilon/2},
 Print @ RealDigits[{nextafter[x, -1], x, nextafter[x, 1]}, 2];
 Differences[{nextafter[x, -1], x, nextafter[x, 1]}]
 ]
(*
{{{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 0}, 7},
 {{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 1}, 7},
 {{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}, 8}}
{1.42109*10^-14, 1.42109*10^-14}
*)

Block[{x = 128.},
 Print @ RealDigits[{nextafter[x, -1], x, nextafter[x, 1]}, 2];
 Differences[{nextafter[x, -1], x, nextafter[x, 1]}]
 ]
(*
{{{1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
   1, 1, 1, 1, 1, 1, 1, 1, 1}, 7},
 {{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}, 8},
 {{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, 1}, 8}}
{1.42109*10^-14, 2.84217*10^-14}
*)

To deal with subnormal numbers, one has to use Compile, ASAIK.

nextafter = With[{min = $MinMachineNumber},
 Compile[{x, y},
  Block[{e, laste},
    If[x == 0,
       Sign[y] min/2.^52,
    If[x < 1,
       e = 2.^(Ceiling @ Log2 @ Abs[x] + 52);
       laste = e*2^-52,
       laste = e = 2.^Ceiling[Log2 @ Abs[x] - 52]];
    While[x + Sign[y] e != x,
          laste = e; e = e/2.];
    x + Sign[y] laste]],
  RuntimeOptions -> {"CompareWithTolerance" -> False}]]

You'll get an error on $MaxMachineNumber depending on the direction. I hope I haven't tripped over any other boundary traps.

$\endgroup$
  • 2
    $\begingroup$ Have you considered IEEE denormal numbers? $\endgroup$ – kirma May 24 '15 at 14:18
  • $\begingroup$ @kirma Nope. They're hard to get with M functions, but I believe they're possible in Compile. They tend to get converted to arb. prec. numbers, if you operate with them. Do you think the OP wants that? I'm not familiar with the C nextafter. $\endgroup$ – Michael E2 May 24 '15 at 14:52
  • 2
    $\begingroup$ I'm just doubtful of these nasty corner cases. I'm not really aware what typical usage scenarios of nextafter are - I want to stick to integer C programming. :) $\endgroup$ – kirma May 24 '15 at 15:16
  • 1
    $\begingroup$ @kirma, sometimes floating point arithmetic is fun; most of the time, it's a PITA to worry about. ;) $\endgroup$ – J. M. will be back soon May 24 '15 at 15:40
6
$\begingroup$

As it turns out, one can exploit the behavior of Interval[] when applied to a machine-precision number to obtain the previous and next representable machine-precision numbers (thanks to Szabolcs for the fix):

SetAttributes[nextafter, Listable];
nextafter[x_?MachineNumberQ, s_?NumericQ] /; s != 0 :=
         First[Interval[x]][[ -Sign[s - x] ]]

To obtain a result equivalent to Michael's two test cases:

With[{x = 128. - 128*$MachineEpsilon/2},
     {RealDigits[{nextafter[x, 127], x, nextafter[x, 129]}, 2], 
      Differences[{nextafter[x, 127], x, nextafter[x, 129]}]}]

and

With[{x = 128.},
     {RealDigits[{nextafter[x, 127], x, nextafter[x, 129]}, 2], 
      Differences[{nextafter[x, 127], x, nextafter[x, 129]}]}]

Additionally:

nextafter[0., {-1, 1}] === {-1, 1} $MinMachineNumber
   True

A caveat of this function is its inability to deal with subnormals.

$\endgroup$
  • $\begingroup$ Looks good but to do precisely what nextafter does, it should be Sign[s-x]. I realized this by playing a bit with C's nextafter. Michael's also uses s, and not s-x as a relative direction, so this function is equivalent to his. It's maybe interesting to note that this doesn't deal with subnormal numbers either (but C's nextafter does) $\endgroup$ – Szabolcs Aug 15 '15 at 8:17
  • 1
    $\begingroup$ I don't know much about denormal numbers. It appears that Mathematica generally switches to arbitrary precision at the point when it would get into the denormal range. Precision@Last@First@Interval[0.] === MachinePrecision but Precision[0.5 * Last@First@Interval[0.]] === 15.9546. However, if I return a denormal number from C, it can do calculations with it while keeping it as MachinePrecision. $\endgroup$ – Szabolcs Aug 15 '15 at 8:58
2
$\begingroup$

If you have a C compiler you could use C's nextafter directly.

Needs["CCompilerDriver`"]
ClearAll[nextafter]
"
#include \"WolframLibrary.h\"

DLLEXPORT mint WolframLibrary_getVersion() {
    return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
    return LIBRARY_NO_ERROR;
}
DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}

DLLEXPORT int nextafterM(
        WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
) {
    double x = MArgument_getReal(Args[0]);
    double y = MArgument_getReal(Args[1]);

    double result = nextafter(x, y);

    MArgument_setReal(Res, result);

    return LIBRARY_NO_ERROR;
}
";
CreateLibrary[%, "nextafter", "CompileOptions" -> "-Wall"]
nextafter = LibraryFunctionLoad[%, "nextafterM", {Real, Real}, Real]

Which gives:

x = 1.*^-300;
RealDigits[#, 2] & /@ {nextafter[x, -1.], x, nextafter[x, 1.]}
(* {{{1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 0}, -996},
    {{1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1}, -996}
    {{1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0}, -996}}} *)

x = 1.*^300;
RealDigits[#, 2] & /@ {nextafter[x, 0.], x, nextafter[x, 2 x]}
(* {{{1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1}, 997},
    {{1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0}, 997}
    {{1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1}, 997}}} *)
$\endgroup$
2
$\begingroup$

It seems that the most straightforward way goes through RealDigits:

nextAfter[x_Real] := 
 FromDigits[MapAt[1 + # &, MapAt[0*# &, RealDigits[x, 2], 1], {1, -1}], 2] + x
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.