How many floating-point numbers are in the interval [0,1]? (2017) | Hacker News Reader

How many floating-point numbers are in the interval [0,1]? (2017)

lemire.me

113 pointsburntcaramel66 comments

For 32 bit floats, you can skip the math and just test all of them. LLVM will vectorise and unroll this nicely.

    fn main() {
        let start = std::time::Instant::now();
        let total = (0..=u32::MAX)
            .filter(|&x| {
                let f = f32::from_bits(x);
                0. <= f && f <= 1.
            })
            .count();
        println!("total {total} in {:?}", start.elapsed());
    }

    total 1065353218 in 1.364751583s

Edit: Apparently if you move the sum to its own function it runs in 500ms. A bit temperamental.

Edit 2: it's the size of the sum accumulator that makes it slow. The version above is like `.fold(0usize, |a, _| a + 1)`. When I moved it to another function, I cast the return value to u32, so LLVM saw basically `.fold(0u32, |a, _| a + 1)` and could use u32 throughout. Godbolt says the usize version ends up with floats in xmm* registers on x86, which fit 4 32-bit floats, but the u32 version ends up with floats in ymm* registers (8 32-bit floats) and similar half-as-wide behaviour on ARM.

This is 0x3f800002, which should look pretty familiar to people who work with floating point in hex; it's 2 + the representation of 1.0. To understand the lowest order bits, you've got -0.0 (0x80000000) in there, as well as everything from 0.0 (0x00000000) to 1.0 (0x3f800000) inclusive.

This is a larger number than the one cited in the post because the latter only includes "normal" numbers. This includes the "denormals" as well.

It's very easy to generalize this to 64 bits, there are 4,607,182,418,800,017,410 of them.

I did the same thing, recalling "There are Only Four Billion Floats–So Test Them All!" [0]. Except I used Python, so it took over 10 minutes:

    >>> import struct
    >>> count = 0
    >>> for i in range(0x1_0000_0000):
    ...     f = struct.unpack('<f', struct.pack('<I', i))[0]
    ...     if 0 <= f <= 1: count += 1
    ...
    >>> count
    1065353218

[0] https://randomascii.wordpress.com/2014/01/27/theres-only-fou...

My very-slightly-different Python implementation runs in 9m56s minutes on Python 3.11, but in 2m54s using pypy3.9 (all on an M1 Max)

pypy never ceases to amaze me for low-effort performance gains on computational stuff.

A naive implementation in C compiled with Apple's clang (-O2) takes 0.62 seconds, of course.

This should be a bit faster, while still sticking to stock Python:

  import time, array

  n = 0
  start_time = time.time()
  values = bytearray(b"".join(i.to_bytes(4, "big") for i in range(2**24)))
  for prefix in range(256):
      values[::4] = prefix.to_bytes(1, "big") * (2**24)
      arr = array.array("f", values)
      n += sum(1 for x in arr if 0.0 <= x <= 1.0)

  print(f"Found: {n} Time: {time.time() - start_time:.1f} seconds")

It uses an array.array() to convert from a large number of concatenated 4 bytes to floats, rather than going through the struct module one-by-one.

Rather than build all 0x1_0000_0000 values in memory, I process 0x100_0000 at a time, and accumulate the partial sums.

On my laptop this reports "Found: 1065353218 Time: 228.5 seconds", so just under 4 minutes.

A slightly improved algorithm in NumPy takes 5.4 seconds.

    import time, numpy as np

    start_time = time.time()
    arr = np.arange(2**24, dtype=np.uint32)
    incr = np.full(2**24, 2**24, dtype=arr.dtype)

    values = np.frombuffer(arr, np.float32)
    n = ((0.0 < values) & (values <= 1.0)).sum()
    for i in range(1, 256):
        arr += incr
        n += ((0.0 < values) & (values <= 1.0)).sum()

    print(f"Found: {n} Time: {time.time() - start_time:.1f} seconds")

In standard C you have nextafterf(3) that gives the next float and allows to traverse them starting from 0:

    #include <math.h>   // nextafterf
    #include <stdio.h>  // printf

    int main()
    {
             long n = 1;
             float f = 0;
             while (f <= 1)
             {
                      f = nextafterf(f, 2);
                      n = n + 1;
             }
             printf("%ld\n", n);
             return 0;
    }

It runs in ~300ms in java

        long start = System.currentTimeMillis();
        int count = 1;
        for (int intBits = Float.floatToIntBits(0.0f); Float.intBitsToFloat(intBits) <= 1.0f; intBits ++) {
            count++;
        }
        System.out.println(count);
        System.out.println((System.currentTimeMillis() - start) + "ms");

Alternatively, it can be solved at compile time with the streams api

        long start = System.currentTimeMillis();
        long count = IntStream.rangeClosed(Float.floatToIntBits(0.0f), Float.floatToIntBits(1.0f)).count();
        System.out.println(count);
        System.out.println((System.currentTimeMillis() - start) + "ms");

  1065353217
  0ms

At compile time!?

Really cool trick!

I can't tell that'd be possible - I see the arguments could be evaluated at compile time, but knowing IntStream.rangeClosed will be evaluated at compile time is a leap, to me.

Did you know before you wrote the code that'd happen? Is it like C? You cross your fingers and hope the compiler unrolls?

The key is that the length of IntStream.rangeClosed(a, b) is just going to be equal to (b - a) + 1.

So there's actually no reason to even invoke the streams API there.

> So there's actually no reason to even invoke the streams API there.

While technically correct, in practice it is not the case.

First of all (and this relates to flerchin's comment), it's not actually done at "compile time" in the traditional sense (i.e. when you invoke `javac`). Rather, it is done at runtime when the method is compiled by the JIT.

As for this comment, the call to `IntStream.rangeClosed` will not simply be reduced to `(b - a) + 1` as you suggest (and one reason for this is that streams are far more complex internally than regular iterators). In reality, it will just be (potentially) inlined and then further optimized alongside the rest of the method's code.

Edit: the mentioned previous comment was from flerchin, not meese

Edit 2: I might have misunderstood what you were getting at there. You are sort of correct in that the `count` operation on this stream is optimized, but it is still technically going through the streams API.

It's not clear to me how the streams api did it, but it's clear to me now that the count of a range of integers is just the end minus the start. This is a fast operation, and I'm not sure if it was performed at compile time, by the JIT, or as a specification of the streams api.

The optimization is done in the JDK implementation (math shortcut) and then further taken advantage of by the JIT (inlining and constant folding).

How is this done?

The implementation of a `Spliterator` that is used internally by `IntStream.rangeClosed` overrides the `estimateSize` method to return `((long) upTo) - from + last`, and the `characteristics` method to return, among others, the `Spliterator.SIZED` characteristic. That characteristic indicates that `estimateSize` can be trusted to be exact, and when present, `estimateSize` will be used as a shortcut in the implementation of `Stream.count`, which can be seen in `ReduceOps.makeRefCounting`.

It's not actually done at compile time. See my reply to messe for details.

If you do this,

        int count = 1;
        int limit = Float.floatToIntBits(1.0f);
        for (int intBits = Float.floatToIntBits(0.0f); intBits <= limit; intBits++) {
            count++;
        }

it will start out with the same performance as your previous one, and then after a few tries the JVM will optimize it into a constant.

  1065353218, 0.91
  1065353218, 0.69
  1065353218, 0.00
  1065353218, 0.64
  1065353218, 0.64
  1065353218, 0.63
  1065353218, 0.66
  1065353218, 0.63
  1065353218, 0.63
  1065353218, 0.00
  1065353218, 0.00
  1065353218, 0.00
  ...
  1065353218, 0.00

  Process finished with exit code 0

Does 0 ms = "at compile time"? Seeing that the impl. just computes

    ((long) upTo) - from + last;

I would think it just takes less than 1 ms on first execution.

https://github.com/openjdk/jdk/blob/caa841d9a52352a975394e55...

It's not actually done at compile time. See my reply to messe for details.

Also, in the case of your test, it's not even letting the JIT kick in, but rather using a shortcut to calculate the length of the range.

> LLVM will vectorise and unroll this nicely

Curious to know the total compile + run time under different compiler optimizations. For code you only need to run once, I don't see how having the compiler unroll the loops actually saves you any time.

Unrolling loops reduces the cost of loop termination checking. A tight loop like `for (i=0; i < n; i++) { acc += 1; }` is mostly overhead, you are constantly checking to see if you’re done or not. If you were given a loop like this in real life, like “come back to me and check if you’ve made enough doughnuts after every single individual doughnut you make, I’ll tell you when you’re done”, you would just quit your job in frustration, because it is a huge waste of your time.

The compiler can unroll loops like that to `for (i=0; i < n; i+=8) { acc += 8; }` (ignoring the remainder of n/8), which does 8x fewer compare+conditional jump instructions. It has nothing to do with how many times you’ll be running the code; this kind of optimisation reduces the runtime of any single execution.

Thanks for the explanation! So how would that play out in optimizing the code you posted? It will run all 2^32 iterations without checking for completion? How does it know when it's finished, and is that the only "unrolling" optimization?

I had thought that in cases like this (where the output can already be determined at compile time) the compiler could just run the loop to get the answer, therefore avoiding running it every time the code runs, which was what I was picturing when I wrote my earlier comment.

Here is an annotated version of rustc's `--release --target aarch64-apple-darwin` assembly output, for the fast version of my code that runs in 500ms.

Annotated https://gist.github.com/cormacrelf/4b2b37d1d870377d548b3a32a... and full Godbolt https://rust.godbolt.org/z/aKK4TEs3x. You may be interested to see the LLVM Opt Pipeline, any of the green passes (non-zero diffs) that are called LoopUnrollPass or LoopVectorizePass.

You can see that it does a couple of LoopUnrollPasses and LoopVectorizePasses. In the end it results in a transformation from something like this:

    for (uint32_t i = 0; i <= 0xffff; i++) {
        if (in range) acc += 1;
    }

    u32x4 v1 = [1,2,3,4];
    u32x4 v0 = [0,0,0,0], v5 = [0,0,0,0];
    u32x4 v6, v7;
    uint32_t w8 = 0xfff8;
    do {
        v6 = v1 + [4,4,4,4];
        v7 = v1 >= 0 && v1 <= 1;
        v6 = v6 >= 0 && v6 <= 1;
        v0 -= v7; // add 1 to each lane if lane was true in v7
        v5 -= v6; // add 1 to each lane if lane was true in v6
        v1 += [8,8,8,8];
        w8 -= 8;
    } while (w8);
    return element_sum(v0 + v5);

That is basically unrolled x8, and then the loop body is grouped into two groups of four vector ops. You also save a lot of time by vectorising the sum operation over 8 lanes. The only non-vector operations in the loop body are on the w8 counter. The ratio of real work to "increment checking" is 8x better, and because we unrolled multiple ops into the loop body, we were able to vectorise it and run it in parallel.

Re "how does it know when it's finished?": it counts like normal, as you can see there, it uses a register to count iteration. Unrolling does not mean "completely eliminate the termination check" unless your heuristic/supplied unroll threshold is larger than the (known) number of iterations. It generally just does the termination check fewer times by taking bigger strides. It's like walking up stairs two at a time.

The answer re compile-time-full-evaluation: compilers use heuristics and pattern matching, they don't automatically try compile time evaluation for things they don't recognise. Nobody would have bothered to write a heuristic to detect "count floating point values in range" and eval it. There is no compiler pass to arbitrarily try to execute eveything in your program at compile time that's remotely possible.

Edit: [0,1,2,3], not [1,2,3,4]. You get the gist.

Makes a lot of sense, thanks for taking the time to explain!

In Turbo Pascal :-) you have to filter for NaN or else you will end up with a Runtime Error.

  const NaN = $FF shl 23;

  var x, t: longint;
      f: single absolute x;

  begin 
    t := 0;
    x := 0;
    repeat
      inc(t, ord((x and NaN <> NaN) and (0<=f) and (f<=1)));
      inc(x);
    until x = 0;
    WriteLn('total: ', t);
  end.

total: 1065353218 (in ca. 40 seconds)

If you write it as an SQL query in ClickHouse, it will take only 0.170 seconds:

    SELECT sum(reinterpretAsFloat32(number) BETWEEN 0.0 AND 1.0) FROM numbers_mt(0x100000000)

    1 row in set. Elapsed: 0.170 sec. Processed 2.39 billion rows, 19.08 GB (14.02 billion rows/s., 112.13 GB/s.)

> 1,056,964,610. There are 4,294,967,296 possible 32-bit words, so about a quarter of them are in the interval [0,1]. Isn’t that interesting? Of all the float-pointing point numbers your computer can represent, a quarter of them lie in [0,1]. By extension, half of the floating-point numbers are in the interval [-1,1].

I don't think it's that strange. It's a scientific notation with the signed exponent around 0, so (barring special cases like NaN and infinite), it will be divided roughly in half around the pivot of 2^0 for both signs.

Surely this was by design?

At least when solving physics problems, you'll frequently scale the equations so the actual values in the problem are around -1..1. Thus having extra precision around -1..1 when solving such problems would be beneficial.

For example, if you're trying to solve the motion of a mass on a spring that gets an initial kick, rather than representing the spring length in meters, you divide by some characteristic length[1] and get a non-dimensional length of order 1. If you want the actual length, you just multiply the non-dimensional length by the characteristic length.

[1]: https://hplgit.github.io/scaling-book/doc/pub/book/html/._sc...

>Surely this was by design?

It is indeed, and one of the first things you learn in scientific computing is to make your values fit mostly within [-1, 1].

> Thus having extra precision around -1..1 when solving such problems would be beneficial.

Floating point numbers have the same precision (relative precision) everywhere.

Yeah it can be hard to internalize what this means and when it matters, but basically you get about 7 decimal places of precision regardless of the numbers size.

So between [-1, 1] you can represent a number like 0.1234567 and then 0.1234568 for a min delta of 0.000001.

But around a billion, you can only represent 1,234,567,000 and then 1,234,568,000 for a min delta of 1000.

These are not exactly right numbers, just rough estimates to get the idea, but the point is if you are trying to do something like add 1 centimeter to a position that is 10,000 kilometers from an origin, you're adding 0, no matter how many times you do it, it will never increment.

Right, but my point is that you scale the physics problems because it makes them more general. In my example, if you changed the length of the spring you don't have to recalculate the whole solution, you just calculate the new characteristic length for your new spring length and "undo" the scaling with this new characteristic length using the previously computed solution.

Thus, since you know a lot of problems will have values around order 1 it makes sense to design the encoding such that you get extra absolute precision for order 1 numbers.

Somewhat desirable, as it means you can get the reciprocal of any number without a loss of precision, right?

Others have previously pondered this, and the related issue of how to randomly select floats in [0,1). I found this to be a helpful account:

https://mumble.net/~campbell/2014/04/28/uniform-random-float

and there's working C code too:

http://mumble.net/~campbell/2014/04/28/random_real.c

I've previously written an algorithm that generates random floats in any [a,b], which can generate all possible floating point values, including subnormals, with the proper probability, and does so quickly for any choice of a and b. [0]

[0] https://github.com/camel-cdr/cauldron/blob/main/cauldron/ran...

Thanks for sharing this. I was worried when I saw the rejection test at the end (around line 1749) that this could be slow, but I guess that for cases like a=0 and b=1, your code avoids that rejection ever happening, is that right? So the rejection only kicks in when the bounds are in the middle of range of one fixed exponent?

The rejection case should only ever be executed for ranges that differ in more than one exponent. (This is handles by the cases above)

The generated exponent starts of at the max exponent of a and b, and is logarithmically decremented (I'm not sure if this is the correct word for it, but decremented until a coin flip hits tails, would be a better way to explain it).

So the exponent and the fraction range can't ever get really small, since there is always at least one exponent that accepts the full range of fractions.

Edit: Actually, looking at it again, the second special case doesn't apply to subnormals, so while generating in the range

    float x = 2.3509886E-38f; // exp=00000001,frac=11111111111111111111111
    float min = nextafterf(nextafterf(x, 0), 0);
    float max = nextafterf(nextafterf(x, 1), 1);

works almost instantly (on average 4 iterations), the current code does indeed have a bottle neck for:

    float x = 1.1754942E-38f; // exp=00000000,frac=11111111111111111111111
    float min = nextafterf(nextafterf(x, 0), 0);
    float max = nextafterf(nextafterf(x, 1), 1);

It's not insanely-minutes-long-stall-bad, but takes about 0.02 seconds to generate one number in the above range with a good prng (opposed to 2e-08 seconds for the first case)

I'll look into fixing this.

Edit 2: I've added a commited a comment about this for now. I wonder why my test suite didn't catch this.

Thanks for the detailed info. I'm glad it's not me trying to craft this code. Good luck with the optimizations.

float32: 1 for exactly 1.0 (exponent=127,mantissa=0) + 127 * 2^23 for [0,1) (127 distinct exponent values (0 through 126) and any bit pattern in the mantissa; this includes +0 and denormals) + 1 for negative zero = 1065353218 = 0x3f800002 + 2, i.e. two more than the representation of 1.0, because everything up to and including that (as an unsigned integer) represents something in [0,1], plus the missing negative zero.

One could make an argument that negative zero is not in [0, 1]. Usually, negative zero is the result of operations like `(small negative number) / (big positive number)`, whose result at infinite precision would not be inside [0, 1]. Negative zero is equal to zero, though, so maybe it is inside [0, 1].

[0,1] = {x | 0 <= x and x <= 1}.

"<=" does not depend on where x came from or what it would have been had something else happened that didn't happen. There's no "maybe."

Are we talking about the interval [0, 1] in R or about [0.0f, 1.0f] (the interval in the space of floats)?

I'm pretty sure that your definition of an interval is correct even when the set doesn't have a total order, but if we're talking about floats whose values are in [0, 1] in the space of real numbers, that may not include -0.0f, since floating point numbers usually represent one number in a range of real values, and -0.0f represents -epsilon where epsilon is between 0 and the real value of the smallest denorm.

> I'm pretty sure that your definition of an interval is correct even when the set doesn't have a total order

Yes, and in the usual ordering on floats that's actually implemented by hardware, 0.0f == -0.0f, so -0.0f is in [-0.0f, 0.0f].

> since floating point numbers usually represent one number in a range of real values, and -0.0f represents -epsilon where epsilon is between 0 and the real value of the smallest denorm.

That's a valid interpretation but apparently the standard ordering wasn't specified from that perspective.

* I meant in [0.0f, 1.0f], of course

* 0x3f800000 + 2, of course

I apply the above theory for generating random numbers in JavaScript, based on given inputs (server seed, client seed, and nonce). JavaScript uses double precision floating point numbers, meaning 52 bits for the mantissa.

  const hash = crypto
    .createHmac('SHA256', serverSeed)
    .update(`${ clientSeed }:${ nonce }`)
    .digest('hex');
  
  const significant = hash.substring(0, 52 / 4);
  const integer = parseInt(significant, 16);
  const float = integer / (2 * 52);
  
  return { float, hash, };

Presumably an offshoot of the discussion around floating point gotchas in gaming: https://news.ycombinator.com/item?id=35539595

Tangentially, it may be interesting to think how this differs for other float formats. Such as: Google Brain's bfloat16, Nvidia's TensorFloat

Relevant: https://stackoverflow.com/questions/43153699/how-can-smalles...

Discussed at the time:

How many floating-point numbers are in the interval [0,1]? - https://news.ycombinator.com/item?id=13759458 - Feb 2017 (153 comments)

So is it 1,056,964,610 (TFA)

or is it 1,065,353,218 (comments) and why the difference?

TFA counted only normal numbers and zero.

The commenters counted also subnormals.

Can I see a graph to know how many of them are in any interval?

https://fabiensanglard.net/floating_point_visually_explained...

Not a graph, but gets much closer to answering this.

I always have trouble re-finding the one I remember, but here's a pretty good visualization:

https://www.researchgate.net/figure/The-distribution-of-repr...

As many as you want depending on how many bits. Floats are made to approximate real numbers. There is an uncountably infinite number of real numbers on any open interval of real numbers. Since (0,1)\in[0,1] there you go. Just keep adding bits.

The question in the title is obviously something explored in the article and to be understood in the context of the article, not to be answered with a glib non-answer by some random HN commenter who evidently didn’t even take a look at the article.

yes, poor headline writing is my fault. and it's not addressed in the article, so I think it's you that didn't read it. It's not even remotely addressed in the article and in fact, the what I posted is an answer to this comment at the end about changing to 64 bits. if he addressed this he would know the answer, and he didn't appear to.

So you chose to just decide I was being "glib" and slander me before understanding the context? Have a weekend. I think you unironically need to touch grass.

I really hate this website sometimes.

About 1/4 of all of them.

A grey link already for me xD

Reply on news.ycombinator.com