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].
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 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.
> 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.
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.
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.
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.
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;
}
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?
> 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`.
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.
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;
}
to
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.
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.
I feel the world (okay, maybe just hn :p) would be a better place if people generally kept the “I tried this with GPT and this is what it said” comments to GPT-related posts.
ChatGPT n that are currently "disruptive". Disruptive to the point that I already have colleagues who get close to palpitations if they can't access it. Within mere weeks or a few months their entire modus operandi has changed. There are major issues (1) with LLM, not least their inscrutability and habit of hallucinating.
Soon we might be getting all whizzed up over "hyperdimensional computing" (2) which, let's face it has at least sorted out how to sound like a cool kid. However it is possible to work out how a hyperdimensional ... thingie "works" because of the way it is built and that is a powerful thing. ChatGPT n that are a bunch of weights and who knows what that actually means. No doubt we will have more stuff and frankenAI and who knows what else but it will be exciting.
I think things have only just started to get weird. It looks to me like a real paradigm shift.
To understand what the hell is happening we need to discuss things and that means discussing the old with the new - hence my comment. Feel free to disagree but please engage.
That can be true, simultaneously with the repetitive “I copy/pasted this into chatGPT, here's what it said!” comments being boring. They don't add anything to the discussion.
If the grandparent provided some analysis of what was output, maybe tied it back into the sorts of errors humans make and how gpt does or does not make dealing with this sort of issue easier, that'd be fine. I may or may not be terribly interested in it, but I wouldn't downvote it. But it didn't.
(edit to add: fwiw, I didn't flag the comment, and upvoted your reply)
I failed to spell it out properly in my OP. How can I persuade (prompt) ChatGPT to derive the correct answer or perhaps I can't.
The French Canadian bloke whose blog is linked here fed an exam to ChatGBT and it apparently passed and he mentioned others too. However my personal experience is rather less convincing. I have got it to write scripts and more but it lacks somewhat yet I still find myself saying please and thank you to it 8)
I'm not sure we even know what we don't know with respect to these things and that is a sure sign for me of a proper paradigm shift.
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'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]
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.
69 comments
[ 4.8 ms ] story [ 163 ms ] thread"<=" 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."
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.
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.
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...
It is indeed, and one of the first things you learn in scientific computing is to make your values fit mostly within [-1, 1].
Floating point numbers have the same precision (relative precision) everywhere.
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.
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.
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 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.
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.
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.
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?
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.
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`.
https://github.com/openjdk/jdk/blob/caa841d9a52352a975394e55...
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.
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.
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.
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.
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:
to 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.
ChatGPT n that are currently "disruptive". Disruptive to the point that I already have colleagues who get close to palpitations if they can't access it. Within mere weeks or a few months their entire modus operandi has changed. There are major issues (1) with LLM, not least their inscrutability and habit of hallucinating.
Soon we might be getting all whizzed up over "hyperdimensional computing" (2) which, let's face it has at least sorted out how to sound like a cool kid. However it is possible to work out how a hyperdimensional ... thingie "works" because of the way it is built and that is a powerful thing. ChatGPT n that are a bunch of weights and who knows what that actually means. No doubt we will have more stuff and frankenAI and who knows what else but it will be exciting.
I think things have only just started to get weird. It looks to me like a real paradigm shift.
To understand what the hell is happening we need to discuss things and that means discussing the old with the new - hence my comment. Feel free to disagree but please engage.
(EDIT - fixup formatting and refs)If the grandparent provided some analysis of what was output, maybe tied it back into the sorts of errors humans make and how gpt does or does not make dealing with this sort of issue easier, that'd be fine. I may or may not be terribly interested in it, but I wouldn't downvote it. But it didn't.
(edit to add: fwiw, I didn't flag the comment, and upvoted your reply)
The French Canadian bloke whose blog is linked here fed an exam to ChatGBT and it apparently passed and he mentioned others too. However my personal experience is rather less convincing. I have got it to write scripts and more but it lacks somewhat yet I still find myself saying please and thank you to it 8)
I'm not sure we even know what we don't know with respect to these things and that is a sure sign for me of a proper paradigm shift.
Not a graph, but gets much closer to answering this.
https://www.researchgate.net/figure/The-distribution-of-repr...
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.
or is it 1,065,353,218 (comments) and why the difference?
The commenters counted also subnormals.
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
[0] https://github.com/camel-cdr/cauldron/blob/main/cauldron/ran...
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
works almost instantly (on average 4 iterations), the current code does indeed have a bottle neck for: 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.
How many floating-point numbers are in the interval [0,1]? - https://news.ycombinator.com/item?id=13759458 - Feb 2017 (153 comments)