I was curious about this. On the one hand, comparing doubles with == is rarely a good idea but, on the other hand, your explanation seems valid.
After some testing I discovered a problem but not with the comparison. The problem is with calculating the halfway value. There are some doubles where their difference cannot be represented as a double:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <float.h>
double random_real(double low, double high, int (*random_bit)(void)) {
if (high < low)
return random_real(high, low, random_bit);
double halfway, previous = low;
while (1) {
halfway = low + (high - low) / 2;
if (halfway == previous)
break;
if (random_bit() & 1)
low = halfway;
else
high = halfway;
previous = halfway;
}
return halfway;
}
int main(int argc, char *argv[]) {
srand(time(NULL));
for (int i = 0; i < 1000000; i++) {
double r = random_real(-DBL_MAX, DBL_MAX, rand);
printf("%f\n", r);
}
}
Actually, the problem is that the algorithm is having to calculate DBL_MAX+DBL_MAX, which of course is going to exceed the maximum value for a double-precision number (by definition). That isn't a very realistic use-case either, but in any case you could just clamp the inputs like so:
double random_real(double low, double high, int (*random_bit)(void)) {
if (high < low)
return random_real(high, low, random_bit);
const double max = DBL_MAX / 2;
if (high > max)
high = max;
const double min = -max;
if (low < min)
low = min;
double halfway, previous = low;
while (1) {
halfway = low + (high - low) / 2;
if (halfway == previous)
break;
if (random_bit() & 1)
low = halfway;
else
high = halfway;
previous = halfway;
}
return halfway;
}
After some testing I discovered a problem but not with the comparison. The problem is with calculating the halfway value. There are some doubles where their difference cannot be represented as a double: