Is it possible to generate random floating point numbers including subnormals, with rand() of math.h?

Question

I would like to generate floating point numbers that include subnormal floating point numbers. Can we use the routine rand() of math.h to achieve this The programming language should be C99. I want this random number to be uniformly distributed from [-1e308 to 1e308].

Note that using `rand` this will always be pseudo-random number, not a true random. — Michał Szydłowski, Jun 30 '15 at 20:08
Generate the sign, fraction and the exponent part using `rand`, and recompose them to `float`. — Eugene Sh., Jun 30 '15 at 20:11
See [glibc](http://www.gnu.org/software/libc/manual/html_node/SVID-Random.html#SVID-Random) — Ross Presser, Jun 30 '15 at 20:13
@this I think you do as the documentation says “uniformly distributed.” Denormal must have a positive probability of appearing. — fuz, Jun 30 '15 at 20:25
@FUZxxl Are you sure about that? Do you exactly mean that subnormals can be generated by rand() of math.h? — zell, Jun 30 '15 at 20:27
@FUZxxl You are misusing uniform distribution in this context. — this, Jun 30 '15 at 20:27
@this the glibc manpage says that `drand48` generates uniformly distributed numbers in the range [0, 1.0). How am I misusing this word? — fuz, Jun 30 '15 at 20:29
@FUZxxl Becase uniform distribution doesn't mean every possible value will be generated, when dealing with floating point. — this, Jun 30 '15 at 20:30
@zell Generate a random floating point value between 0 and 1, and then divide by a large power of 2 to get denormals. — this, Jun 30 '15 at 20:50
Please specify a distribution and a range, otherwise this question is unanswerable. — Oliver Charlesworth, Jun 30 '15 at 20:52
Let's also think pragmatically about this. Assuming you mean uniform in [0,1), you would only expect to see a denormal once every 1E21 years (approximately) even if you were generating 1 billion randoms per second. Are you sure you want this? — Oliver Charlesworth, Jun 30 '15 at 21:06
@this I do not thing one can get subnormal using your approach. One would get 0. — zell, Jun 30 '15 at 21:24
@zell: If you want uniformly-distributed, why have you accepted an answer that very definitely does not give you a uniform distribution? — Oliver Charlesworth, Jun 30 '15 at 21:33
@zell Certainly you want an approximate _logarithmic_ distribution from [-1e308 to 1e308] and not uniform linear one. IOWs there are about as many random numbers generated between 1.0 to 2.0 as there are between 1e6 to 2e6 and well as 1e-6 to 2e-6. — chux - Reinstate Monica, Jun 30 '15 at 21:45
@zell No, you would get a denormal if you pick the power correctly. — this, Jul 01 '15 at 10:20
@this. Would you give more details, if you can, maybe as an answer to the question? Sorry but I do not see a viable way you can get subnormals with your claimed apporach. — zell, Jul 01 '15 at 23:40
@zell Divide a float by a larger float, until you get a denormal. Use a while loop, division operator and printf to see the value of the float.. — this, Jul 02 '15 at 09:12

score 1 · Answer 1 · edited May 23 '17 at 12:21

Writing to one field of a union and reading from another often has troubles, but the following works well-enough on many systems:

This will generate not-a-numbers. NAN

double random_double(void) {
  union {
    unsigned char uc[sizeof (double)];
    double d;
  } u;
  unsigned i;
  for (i=0; i<sizeof u.uc; i++) {
    u.uc[i] = rand();
  }
  return u.d;
}

A way to not return NAN. Also see IEEE 754 specifies an isNaN() predicate which does not raise an exception even when used with signaling NaN

double rand_double_NotNAN(void) {
  union {
    double d;
    unsigned char uc[sizeof(double)];
  } u;
  do {
    for (size_t i = 0; i < sizeof(double); i++) {
      u.uc[i] = rand();
    }
  } while (isnan(u.d));
  return u.d;
}

score 1 · Answer 2 · answered Jul 19 '15 at 16:29

It won't be efficient, but to properly answer the question:

first define a mapping from the integer range (0, MAX_RAND) to the reals. In defining this mapping, decide exactly what distribution you want: uniform over all reals, or logarithmically uniform, or whatever. Write a function that given an integer in the range, returns a double according to this mapping. Now use rand() and feed it into your function.

score 0 · Accepted Answer · answered Jun 30 '15 at 20:16

0

union {
  int a;
  float b;
} Number;

Number.a = rand();
printf("%f\n", Number.b);

this uses the bitmask of a random integer to create a float number with the same bitmask.

answered Jun 30 '15 at 20:16

mch

9,424
2
28
42

2

It seems a reasonable answer. Someone downvoted however. Why? – zell Jun 30 '15 at 20:19
4

I'm not yet a downvoter, but strictly speaking, the behaviour of this program is undefined. – Bathsheba Jun 30 '15 at 20:21
This approach is unportable. Writing to one union member and then reading from another is unspecified behaviour. – fuz Jun 30 '15 at 20:23
1

@Bathsheba can you explain why? it is correct in c99 if `sizeof(int) == sizeof(float)`. It is okay to write to one member of an union and read from another if they occupy the same bytes since c99. – mch Jun 30 '15 at 20:25
2

@mch Because a member can have a trap representation. – this Jun 30 '15 at 20:29
1

What is the distribution of the random numbers generated with this approach? – Oliver Charlesworth Jun 30 '15 at 20:54
Note: as `rand()` generates the value `0...RAND_MAX` which does not include negative `int`, likely half of the possible `float` will not be generated. – chux - Reinstate Monica Jun 30 '15 at 21:15

Is it possible to generate random floating point numbers including subnormals, with rand() of math.h?

3 Answers3