Does Mathematica have an equivalent of C's nextafter?

Question

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.

@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. — David Zhang, May 24 '15 at 04:56
@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. — David Zhang, May 24 '15 at 05:01
Have a look at Ulp et al. in the Computer Arithmetic package... — ciao, May 24 '15 at 05:28
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. — J. M.'s missing motivation, May 24 '15 at 05:38
@J. M. I still don't believe it would; the gap between 1.0e300 and the next floating-point value is a lot larger than $MachineEpsilon. In general, the gaps between consecutive floating-point numbers change with the magnitude of the numbers. — David Zhang, May 24 '15 at 05:41
@ciao Ulp is definitely a step in the right direction, but I'm not sure how to use it to implement nextafter. In particuar, I'm not sure how to handle the transitions between one exponent value and the next; note that 1.0 + Ulp[1.0] == 1.0 while 1.0 - Ulp[1.0] != 1.0. — David Zhang, May 24 '15 at 05:42
You can also consider doing something like {IntegerPart[#1/$MachineEpsilon], #2} & @@ MantissaExponent[x, 2] to split the floating point number to integer-valued components and working further with those. — kirma, May 24 '15 at 05:43

score 7 · Answer 1 · edited Jun 10 '16 at 03:01

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.

@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. — Michael E2, May 24 '15 at 14:52
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. :) — kirma, May 24 '15 at 15:16
@kirma, sometimes floating point arithmetic is fun; most of the time, it's a PITA to worry about. ;) — J. M.'s missing motivation, May 24 '15 at 15:40

J. M.'s missing motivation · Answer 2 · 2015-08-15T12:57:44.047

6

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.

edited Aug 15 '15 at 12:57

answered Aug 15 '15 at 07:30

J. M.'s missing motivation

124,525
11
401
574

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) – Szabolcs Aug 15 '15 at 08:17
1

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. – Szabolcs Aug 15 '15 at 08:58

score 5 · Answer 3 · answered Jun 10 '16 at 07:05

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}}} *)

To compile with Visual Studio <mathl.h> is needed.. – cesar guerra Oct 01 '20 at 17:47 — cesar guerra, Oct 01 '20 at 17:47

score 2 · Answer 4 · answered Jun 10 '16 at 07:23

2

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

answered Jun 10 '16 at 07:23

Alexey Popkov

61,809
7
149
368

Does Mathematica have an equivalent of C's nextafter?

4 Answers4

Linked