What's the closest double to 1.0, that isn't 1.0?

Question

Is there a way to programmatically get the double that is closest to 1.0, but isn't actually 1.0?

One hacky way to do this would be to memcpy the double to a same-sized integer, and then subtract one. The way IEEE754 floating-point formats work, this would end up decreasing the exponent by one while changing the fractional part from all zeros (1.000000000000) to all ones (1.111111111111). However there exist machines where integers are stored little-endian while floating-point is stored big-endian, so that won't always work.

You can't assume that +1 is the same distance (from 1.0) as -1. The interleaving of base 10 and base 2 floating point representations means that the gaps are uneven. — Richard Critten, Aug 06 '16 at 08:35
@Richard: you're right. It is very unlikely that subtracting one ULP will get the, er, "nextbefore" value, because I guess the exponent would have to be adjusted as well. `nextafter()` is the only proper way of achieving what he wants. — Rudy Velthuis, Aug 06 '16 at 10:37
FYI have a read of this blog (not mine): http://www.exploringbinary.com/17-digits-gets-you-there-once-youve-found-your-way/ — Richard Critten, Aug 06 '16 at 11:40
@RudyVelthuis: It does work on every IEEE754 binary floating point format. — Edgar Bonet, Aug 07 '16 at 17:46
@RudyVelthuis: No, I just believe what IEEE 754 says. If you don't, _you_ should try: that's how you build confidence in your understanding of the spec. — Edgar Bonet, Aug 08 '16 at 08:27
Ok, then tell me: what "does work on every IEEE754 floating point format"? It is simply not true that if you decrement the significand that you get the "firstbefore()" value, especially not for 1.0, which has a significand which is a power of two. That means that `1.0000...` binary is decrement to `0.111111....` and to normalize it, you must shift it to the left: `1.11111...` which requires you to decrement the exponent. And then you are 2 ulp away from 1.0. So no, subtracting one from the integral value does NOT give you what is asked here. — Rudy Velthuis, Aug 08 '16 at 09:14
@Edgar: You might want to show your confidence in your reading skills and quote that part of the IEEE 754 spec that says that "it does work on every IEEE754 binary floating point". — Rudy Velthuis, Aug 08 '16 at 09:19
@RudyVelthuis: The actual standard is behind a paywall, and I will gladly read it for you if you buy me a copy. The data format, however, is available from many secondary sources. It is really common knowledge, so you have no excuse for telling me wrong before checking the facts. I'll just give you a hint: 1.0 = 0–01111111111–000...000, 1−ε/2 = 0–01111111110–111...111. I'll let you do the integral subtraction (it's not that hard). — Edgar Bonet, Aug 08 '16 at 20:46
@RudyVelthuis: For all positive, non-∞ FP values, the result of comparing the actual bit representation, as integer, of two FP numbers is the same as the result of comparing their FP values. Consequently, the closest two values to any positive floating-point value, other than DBL_MIN and DBL_MAX, can be computed by writing the FP value to memory, doing an integer increment/decrement, and then reading it back. At least, in IEEE754. Trust me, I've done this. Trouble is, you can't write PORTABLE code to do that because on some architectures, the endian-ness is different, between integers and FP. — jorgbrown, Sep 17 '16 at 01:22
If you don't believe me, run the code: http://ideone.com/kewPwu — jorgbrown, Sep 17 '16 at 01:29

Jarod42 · Accepted Answer · 2016-08-06T08:06:05.523

157

Since C++11, you may use nextafter to get next representable value in given direction:

std::nextafter(1., 0.); // 0.99999999999999989
std::nextafter(1., 2.); // 1.0000000000000002

Demo

edited Aug 06 '16 at 08:06

answered Aug 06 '16 at 07:59

Jarod42

203,559
14
181
302

14

This is also a nice way to increment a double to the next representable integer: `std::ceil(std::nextafter(1., std::numeric_limits::max()))`. – Johannes Schaub - litb Aug 06 '16 at 08:26
45

The next question's going to be "how is this implemented in the stdlib" :P – Lightness Races in Orbit Aug 06 '16 at 11:59
1

My man page for `nextafter` says "CONFORMING TO C99, POSIX.1-2001, POSIX.1-2008. This function is defined in IEC 559 (and the appendix with recommended functions in IEEE 754/IEEE 854)." So you may be able to get it via C even if you're 12 years behind C++11 :) – Ben Millwood Aug 06 '16 at 14:47
19

After reading @LightnessRacesinOrbit's comment I got curious. [This is how glibc implements `nextafter`](http://bazaar.launchpad.net/~vcs-imports/glibc/master/view/head:/math/s_nextafter.c), [and this is how musl implements it](https://git.musl-libc.org/cgit/musl/tree/src/math/nextafter.c) in case anyone else wants to see how it's done. Basically: raw bit twiddling. – Cornstalks Aug 07 '16 at 03:02
3

@Cornstalks: I am not surprised it's down to bit twiddling, the only other option would be having CPU support. – Matthieu M. Aug 07 '16 at 23:33
6

Bit twiddling is the only way to do it properly, IMO. You could do a lot of test tries, trying to slowly approach it, but that could be very slow. – Rudy Velthuis Aug 08 '16 at 10:57
1

@Cornstalks For positive normal doubles it's just the bit pattern, cast to int, plus one. Which is pretty cool, actually. – Dan Aug 14 '16 at 14:55

score 26 · Answer 2 · answered Aug 06 '16 at 16:02

In C and C++, the following gives the closest value to 1.0:

#include <limits.h>

double closest_to_1 = 1.0 - DBL_EPSILON/FLT_RADIX;

Note however that in later versions of C++, limits.h is deprecated in favour of climits. But then, if you are using C++ specific code anyway, you can use

#include <limits>

typedef std::numeric_limits<double> lim_dbl;
double closest_to_1 = 1.0 - lim_dbl::epsilon()/lim_dbl::radix;

And as Jarod42 writes in his answer, since C99 or C++11 you can also use nextafter:

#include <math.h>

double closest_to_1 = nextafter(1.0, 0.0);

Of course in C++ you can (and for later C++ versions should) include cmath and use std::nextafter instead.

score 23 · Answer 3 · answered Aug 06 '16 at 08:01

23

In C, you can use this:

#include <float.h>
...
double value = 1.0+DBL_EPSILON;

DBL_EPSILON is the difference between 1 and the least value greater than 1 that is representable.

You'll need to print it to several digits in order to see the actual value.

On my platform, printf("%.16lf",1.0+DBL_EPSILON) gives 1.0000000000000002.

answered Aug 06 '16 at 08:01

barak manos

29,648
10
62
114

10

So that wont work for some other values than `1.` as `1'000'000` [Demo](http://ideone.com/xdeu9e) – Jarod42 Aug 06 '16 at 08:14
9

@Jarod42: You're right, but OP asks specifically about `1.0`. BTW, it also gives the closest value greater than 1, and not the absolute closest value to 1 (which is possibly smaller than 1). So I agree that this is a partial answer, but I thought it could nevertheless contribute. – barak manos Aug 06 '16 at 08:17
@LưuVĩnhPhúc: I give precision about the restriction of the answer, and the closest in the other direction. – Jarod42 Aug 06 '16 at 08:20
8

This does not give the *closest* double to 1.0, as (assuming base 2) the double right *before* 1.0 is only half as far away as the double right *after* 1.0 (which is the one you calculate). – celtschk Aug 06 '16 at 15:38
@celtschk: You're right, I've explained that in the comment above. – barak manos Aug 06 '16 at 15:42
The closest value to one is indeed smaller than one, as @barakmanos sugests, since [double](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) values between 0.5 and 1.0 have twice the precision of values between 1.0 and 2.0. – HelloGoodbye Aug 06 '16 at 23:41
@HelloGoodbye: barakmanos is the same user who wrote this answer. – barak manos Aug 07 '16 at 05:10
The question is very explicit: “closest double to 1.0, that isn't 1.0”, and this answer is just wrong. – Edgar Bonet Aug 07 '16 at 17:49
How does this work at all since 1.0 is not representable with absolute accuracy? Or is 1.0 representable in memory with absolute accuracy? - Isn't that the point of the question, that 1.0 is *not* representable in memory with absolute accuracy? Does this even make sense? What? – FreelanceConsultant Aug 09 '16 at 18:51
@user3728501: Sorry, I do not understand your question. How does **what** work? – barak manos Aug 09 '16 at 18:56
@barakmanos `double value = 1.0 + DBL_EPSILON`; - `DBL_EPSILON` has some binary value in memory, probably something like `0 00000000000 0000000000000000000000000000000000000000000000000001` (0 sign, 0 exponent, all fraction set to zero except last bit? - maybe this isn't quite right but you see the principle of the idea)... but essentially my question can be answered by considering *is 1.0 representable exactly in floating point format?*... At first I thought "no", however I now think "yes", and I *think* that binary representation I just gave is 1.0? - In which case what is the representation – FreelanceConsultant Aug 09 '16 at 19:46
of `DBL_EPSILON` in memory as a 64 bit float? – FreelanceConsultant Aug 09 '16 at 19:47
@user3728501: Why don't you give it a try and see exactly how it is represented in memory? Here is a piece of code that breaks strict aliasing rule but should still give you the answer: `double x = 1.0; printf("%X",*(uint64_t*)&x);`. – barak manos Aug 10 '16 at 04:53
Could do that I suppose, but like everyone else in the world, we look for the fastest solution - ie a simple yes/no answer to the question is 1.0 absolutely representable in memory will do. Having thought about it further I conclude the answer is yes. – FreelanceConsultant Aug 10 '16 at 10:11

score 4 · Answer 4 · answered Aug 06 '16 at 08:17

4

In C++ you can also use this

1 + std::numeric_limits<double>::epsilon()

answered Aug 06 '16 at 08:17

phuclv

37,963
15
156
475

1

Like barak manos' answer, this will not work for any value other than 1. – zwol Aug 06 '16 at 16:32
2

@zwol tehnically for typical binary floating-point implementations it will work for any value between 1 and 2-epsilon. But, yes, you're right that it's only guaranteed to apply to 1. – Random832 Aug 06 '16 at 21:24
9

Technically, it doesn't work for 1, since the closest number to 1 is the number right before 1, not the one right after it. [double](https://en.wikipedia.org/wiki/Double-precision_floating-point_format)'s precision between 0.5 and 1 is twice as high as its precision between 1 and 2, hence the number right before 1 ends up closer to 1. – HelloGoodbye Aug 06 '16 at 23:46

What's the closest double to 1.0, that isn't 1.0?

4 Answers4

Linked