18
$\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.

Improve this question
$\endgroup$
8
  • $\begingroup$ x + Sign[y] $MachineEpsilon should work, yes? $\endgroup$
    – J. M.'s missing motivation ♦
    Commented May 24, 2015 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
    Commented May 24, 2015 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
    Commented May 24, 2015 at 5:01
  • 2
    $\begingroup$ Have a look at Ulp et al. in the Computer Arithmetic package... $\endgroup$
    – ciao
    Commented May 24, 2015 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.'s missing motivation ♦
    Commented May 24, 2015 at 5:38

4 Answers 4

Reset to default
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.

Improve this answer
$\endgroup$
4
  • 2
    $\begingroup$ Have you considered IEEE denormal numbers? $\endgroup$
    – kirma
    Commented May 24, 2015 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
    Commented May 24, 2015 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
    Commented May 24, 2015 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.'s missing motivation ♦
    Commented May 24, 2015 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.

Improve this answer
$\endgroup$
2
  • $\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
    Commented Aug 15, 2015 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
    Commented Aug 15, 2015 at 8:58
5
$\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}}} *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ To compile with Visual Studio <mathl.h> is needed.. $\endgroup$
    – cesar guerra
    Commented Oct 1, 2020 at 17:47
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
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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