[spyrja]

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;
  }

[WantonQuantum]

Fixed it!

  #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/2 - 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);
    }
  }

[spyrja]

Yep, that works too. You likely won't ever need a random number that large, of course, but if you did that would be the way to do it.