Uniformly random doubles from [0, 1) in C

A place to discuss the implementation and style of computer programs.

Moderators: phlip, Moderators General, Prelates

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Uniformly random doubles from [0, 1) in C

Postby Qaanol » Tue Apr 28, 2015 5:01 am UTC

Could someone please check that I’m doing the right thing (or recommend a better way) to rapidly fill a contiguous array with high-quality double-precision values distributed uniformly in [0, 1)?

My approach is to generate a random mantissa using arc4random, set the sign and exponent bits so the result is in [1, 2), and subtract 1. This is on an Intel Mac running OS X, so sizeof(long) == sizeof(double) == 8:

Code: Select all

void randomDoubleArray(double * v, const int n) {
    // Fill the array with random bytes
    arc4random_buf(v, n * sizeof(double));
   
    // Twiddle the bits
    for (int i = 0; i < n; i++) {
        long * const p = (long *)&v[i];
       
        *p &= 0x000fffffffffffffl; // Clear the sign and exponent bits
        *p |= 0x3ff0000000000000l; // Set exponent so v[i] is in [1,2)
        v[i] -= 1.0;
    }
}
wee free kings

User avatar
hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby hotaru » Wed Apr 29, 2015 12:13 am UTC

i see one thing that might be an issue... there are 262 doubles in [0, 1), but only 252 in [1, 2).
I'd probably do something like this, but I'm not sure if it's actually better:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 assert(sizeof(double) * CHAR_BIT == 64);
  arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i;
    uint64_t e = 1022;
    while(&& arc4random_uniform(2)) --e;
    *>>= 12;
    *|= e << 52; }}
Last edited by hotaru on Wed Oct 26, 2016 12:57 am UTC, edited 1 time in total.

Code: Select all

factorial product enumFromTo 1
isPrime n 
factorial (1) `mod== 1

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby Qaanol » Wed Apr 29, 2015 3:49 am UTC

hotaru wrote:i see one thing that might be an issue... there are 262 doubles in [0, 2), but only 252 in [1, 2).
I'd probably do something like this, but I'm not sure if it's actually better:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 assert(sizeof(double) * CHAR_BIT == 64);
  arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i;
    uint64_t e = 1022;
    while(&& arc4random_uniform(2)) --e;
    *>>= 12;
    *|= e << 52; }} 


Yes, there are far more doubles in [0, 1) than [1, 2), but I don’t necessarily need to make it possible to hit every representable double. Instead I am making it equally-likely to hit all doubles of the form k/252 where k is an integer in [0, 252−1).

Your method is quite clever, and appears to make every double in [0, 1) possible with probability proportional to the length of the interval from it to the next double. However, it does so at the cost of additional calls to arc4random for every entry. The first up-to-12 of those calls could be avoided by using the high bits of *p as stand-ins for the initial dozen arc4random_uniform calls though, so extra calls are only needed 1 in 4096 times.

Let’s see if I’ve got this…

Code: Select all

#define TWELVE_HIGH_BITS_64 (0xfff0000000000000l)

void randomDoubleArray(double * v, const int n) {
    assert(sizeof(double) * CHAR_BIT == 64);
    arc4random_buf(v, n * sizeof(double));
   
    for(int i = 0; i < n; i++) {
        uint64_t * const p = (uint64_t *)(v+i);
        uint64_t h = *p & TWELVE_HIGH_BITS_64;
        *p &= ~TWELVE_HIGH_BITS_64;  // Clear sign and exponent
       
        // Each exponent is as likely as all lower exponents,
        // so flip coins until tails, decrementing each time.
        // Use the bits of h for the first twelve coin flips.
       
        uint64_t e = 1022;
        if (h == TWELVE_HIGH_BITS_64) {
            e = 1010;  // First 12 flips were heads
            while(arc4random_uniform(2) && e) { e--; }
        } else {
            // Count the leading 1s of h
            while (h & (1l << 63))   { h <<= 1; e--; }
        }
       
        *p |= (e << 52);  // Set e as exponent of v[i]
    }
}


For comparison, running all three operations on separate 10,000,000-element arrays, here are typical times (in seconds) that I’m seeing:
0.393 : My function in this post
1.018 : Hotaru’s function in the previous post
0.306 : My function in the original post

I’m still not completely certain which approach I prefer. For the project I’m working on right now, most of the time I’ll actually need a double in (0, 1] so I can take its natural log. With my original function I just have to replace “v[i] -= 1.0;” with “v[i] = 2.0 - v[i];” to get that range. With the new version, I could do something like “v[i] = nextafter(v[i], DBL_MAX);”.

However, even with that modification to generate in (0, 1], the new version can still produce subnormal numbers, which may be a problem for applications that need to take the reciprocal, because that overflows and returns “inf” for small-enough subnormal values. For my own purposes, the natural log of the smallest positive (subnormal) double, log(0x1p-1074), does not overflow and in fact equals -0x1.74385446d71c3p+9, or approximately -744.44.
wee free kings

User avatar
hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby hotaru » Wed Apr 29, 2015 11:07 am UTC

Qaanol wrote:Your method is quite clever, and appears to make every double in [0, 1) possible with probability proportional to the length of the interval from it to the next double. However, it does so at the cost of additional calls to arc4random for every entry. The first up-to-12 of those calls could be avoided by using the high bits of *p as stand-ins for the initial dozen arc4random_uniform calls though, so extra calls are only needed 1 in 4096 times.

it can actually be improved more than that, by getting 32 bits at a time from arc4random() instead of 1 at a time from arc4random_uniform(), so one call to arc4random is only needed 1 in 4096 times, a second is only needed 1 in 17592186044416 times, a third is only needed one in 75557863725914323419136 times, etc.:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 assert(sizeof(double) * CHAR_BIT == 64);
  arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i,
             e = 1022;
    uint32_t r = *>> 52;
    int i = 12;
    while(&& r & 1)
    { --e; --i; r >>= 1;
      if(!i)
      { r = arc4random();
        i = 32; }}
    *&= 0x000ffffffffffffful;
    *|= e << 52; }} 


Qaanol wrote:I’m still not completely certain which approach I prefer. For the project I’m working on right now, most of the time I’ll actually need a double in (0, 1] so I can take its natural log. With my original function I just have to replace “v[i] -= 1.0;” with “v[i] = 2.0 - v[i];” to get that range. With the new version, I could do something like “v[i] = nextafter(v[i], DBL_MAX);”.

i'm pretty sure using nextafter like that would make the results (very slightly) non-uniform.

Qaanol wrote:However, even with that modification to generate in (0, 1], the new version can still produce subnormal numbers, which may be a problem for applications that need to take the reciprocal, because that overflows and returns “inf” for small-enough subnormal values. For my own purposes, the natural log of the smallest positive (subnormal) double, log(0x1p-1074), does not overflow and in fact equals -0x1.74385446d71c3p+9, or approximately -744.44.

if generating subnormals is a problem, replacing "e" in the while statement with "e > 1" should fix that, and produce values in [2-1022, 1).

Code: Select all

factorial product enumFromTo 1
isPrime n 
factorial (1) `mod== 1

KnightExemplar
Posts: 5494
Joined: Sun Dec 26, 2010 1:58 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby KnightExemplar » Wed Apr 29, 2015 12:40 pm UTC

However, even with that modification to generate in (0, 1], the new version can still produce subnormal numbers, which may be a problem for applications that need to take the reciprocal, because that overflows and returns “inf” for small-enough subnormal values. For my own purposes, the natural log of the smallest positive (subnormal) double, log(0x1p-1074), does not overflow and in fact equals -0x1.74385446d71c3p+9, or approximately -744.44.


Can't you do a 2nd pass for subnormals?

Basically, a two pass algorithm that does:
Pass 1: What you did in your original post
Pass 2: Poisson distribution for normals. Basically the pseudocode would be:

for(i=0; i<array_length; i+= poisson_distribution(1/2^52)){
array[i] = random_subnormal();
}

And random_subnormal() would be just a simple random-bits with the exponent set to... whatever the magic number for subnormals was. Poisson distribution would not have to be "efficient", because it would be generated _very_ rarely. In fact, subnormals will only happen once every 32 Petabytes on the average (once every 2^52, and each double takes up 8 bytes). But if the Poisson distribution did its magic, even the subnormals will be generated at the correct rate with the above pseudocode.

The other thing you can do is say declare your code is correct "on the average" for all buffer sizes of less than 1TB. (Which roughly has a probability of generating a subnormal at 2^40 / 2^52 or 1/2^12).

------------------------

Something I would be worried about though is that smaller and smaller numbers have less precisions. One in eight numbers (ie: all numbers 0.125 and lower) generated by your example will only have 2^49 bits generated (bottom 3 bits would be left to the rounder). So your numbers get less and less precision the smaller they get. Maybe doing a 2nd pass with some extra random bits would do the trick? I think your methodology is very good at choosing the correct exponent efficiently. But tacking on some extra bits on the end can't hurt.
First Strike +1/+1 and Indestructible.

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby Qaanol » Wed Apr 29, 2015 5:21 pm UTC

hotaru wrote:it can actually be improved more than that, by getting 32 bits at a time from arc4random() instead of 1 at a time from arc4random_uniform(), so one call to arc4random is only needed 1 in 4096 times, a second is only needed 1 in 17592186044416 times, a third is only needed one in 75557863725914323419136 times, etc.:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 assert(sizeof(double) * CHAR_BIT == 64);
  arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i,
             e = 1022;
    uint32_t r = *>> 52;
    int i = 12;
    while(&& r & 1)
    { --e; --i; r >>= 1;
      if(!i)
      { r = arc4random();
        i = 32; }}
    *&= 0x000ffffffffffffful;
    *|= e << 52; }} 


Yeah, I realized that and implemented it shortly after I posted. Be careful with your two different variables “i” though. My version looks like this inside the “if” from my previous post:

Code: Select all

            e = 1010;  // First 12 flips were heads
            bool proceeding = true;
            while (proceeding && e) {
                uint32_t b = arc4random();
                proceeding = (== -1);
                while (& 1) { b >>= 1; e--; }
                e = (1024) ? e;
            } 


However, testing reveals no significant difference in speed from the previous. I mean, the difference can only possibly come up 1 in 8,192 times, and half of those times it’s just saving a single extra function call.

hotaru wrote:i'm pretty sure using nextafter like that would make the results (very slightly) non-uniform.

No more so than it already is. Currently, what we are doing has the same outcome as this process:

1. Partition [0, 1) into half-open subintervals bounded by the representable doubles.
2. Choose a real in [0, 1) uniformly at random.
3. Find the lower bound of the subinterval the chosen number is in.
4. Return that lower bound.

Inserting a single “nextafter” transforms it into:

1. <as above>
2. <as above>
3. Find the upper bound of the subinterval the chosen number is in.
4. Return that upper bound.

Essentially, as an approximation for a true uniform distribution, the “standard” version is slightly on the low side, and the “nextafter” version is slightly on the high side.

If we wanted to make the standard version as “perfect” as possible, we’d have to do something pretty subtle. I’m not even sure exactly what. If there were only one bit so the options were 0 and ½, then we’d always have to pick ½ to make the average correct, but that would blow up our variance. If the options were 0, ¼, ½, and ¾, then…I haven’t thought it through and I’m not about to.

hotaru wrote:if generating subnormals is a problem, replacing "e" in the while statement with "e > 1" should fix that, and produce values in [2-1022, 1).

True, though that prevents both 0 and 1 from being rolled.

KnightExemplar wrote:Can't you do a 2nd pass for subnormals?

Basically, a two pass algorithm that does:
Pass 1: What you did in your original post
Pass 2: Poisson distribution for normals. Basically the pseudocode would be:

for(i=0; i<array_length; i+= poisson_distribution(1/2^52)){
array[i] = random_subnormal();
}

And random_subnormal() would be just a simple random-bits with the exponent set to... whatever the magic number for subnormals was. Poisson distribution would not have to be "efficient", because it would be generated _very_ rarely. In fact, subnormals will only happen once every 32 Petabytes on the average (once every 2^52, and each double takes up 8 bytes). But if the Poisson distribution did its magic, even the subnormals will be generated at the correct rate with the above pseudocode.

The other thing you can do is say declare your code is correct "on the average" for all buffer sizes of less than 1TB. (Which roughly has a probability of generating a subnormal at 2^40 / 2^52 or 1/2^12).

I’m not quite sure I follow. If my original function doesn’t produce any doubles with granularity finer than 1/252, why would I add to it the ability to generate subnormals specifically?

In any case, the reason I am taking the natural log of numbers in (0, 1] in the first place is to generate exponentially-distributed random variates (either for their own use, or as a step in the Box-Muller transform to generate normal variates). So it seems a little silly to use a Poisson distribution in the production of the initial uniform variates.

I suppose I could switch to a ziggurat algorithm for more complex distributions, but frankly I’d rather not at this point.

KnightExemplar wrote:Something I would be worried about though is that smaller and smaller numbers have less precisions. One in eight numbers (ie: all numbers 0.125 and lower) generated by your example will only have 2^49 bits generated (bottom 3 bits would be left to the rounder). So your numbers get less and less precision the smaller they get. Maybe doing a 2nd pass with some extra random bits would do the trick? I think your methodology is very good at choosing the correct exponent efficiently. But tacking on some extra bits on the end can't hurt.

While you are technically correct, the 49 bits in your example represent exactly the value that I want represented. There is no “rounding” when going from 1.125 to 0.125. Yes, the latter has more trailing zeros, but both represent exact values—namely some multiple of 2−52.

All multiples of 2−52 are equally likely, and none get skipped, so I don’t see the problem. That said, the more I think about it the more I like Hotaru’s approach. Especially for my purposes, where finer granularity at smaller values corresponds directly to a smoother distribution of output from the natural log, and thus a more representative mix of exponential variates.
Last edited by Qaanol on Wed Apr 29, 2015 9:45 pm UTC, edited 2 times in total.
wee free kings

KnightExemplar
Posts: 5494
Joined: Sun Dec 26, 2010 1:58 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby KnightExemplar » Wed Apr 29, 2015 7:05 pm UTC

Qaanol wrote:
KnightExemplar wrote:Something I would be worried about though is that smaller and smaller numbers have less precisions. One in eight numbers (ie: all numbers 0.125 and lower) generated by your example will only have 2^49 bits generated (bottom 3 bits would be left to the rounder). So your numbers get less and less precision the smaller they get. Maybe doing a 2nd pass with some extra random bits would do the trick? I think your methodology is very good at choosing the correct exponent efficiently. But tacking on some extra bits on the end can't hurt.

While you are technically correct, the 49 bits in your example represent exactly the value that I want represented. There is no “rounding” when going from 1.125 to 0.125. Yes, the latter has more trailing zeros, but both represent exact values—namely some multiple of 2−52.

All multiples of 2−52 are equally likely, and none get skipped, so I don’t see the problem. That said, the more I think about it the more I like Hotaru’s approach. Especially for my purposes, where finer granularity at smaller values corresponds directly to a smoother distribution of output from the natural log, and thus a more representative mix of exponential variates.


Yeah, my methodology was to try to combine the speed advantage to your methodology and hotaru's "more complete" coverage of the double precision area. I recognize that your methodology is uniform and extremely quick.

But I recognize Hotaru's methodology as "more uniform" but slower, probably due to the conditional statement (the while loop) in his procedure. So despite the fact that it barely adds any computation, you can see it in the benchmark that its dramatically slower (3x or so). So instead of doing hotaru's methodology, I was thinking maybe using your method as a baseline, and then fixing it up "post-hoc" with a 2nd pass based off of the Poisson distribution.

There's something inelegant about using the Poisson distribution here. So you're right from that perspective.

------------------

Another thing I was thinking was potentially "parallelizing" Hotaru's method. He's essentially counting the number of zeros in a row. So if you can generate a large number of exponents in parallel, it'd effectively be a faster version of his method.

In essence, use the popular "parallel count bits" method to generate a whole bunch of exponents from the random stream of bytes. And then perhaps do a 2nd pass where you set the exponents? But this isn't a counting problem... its more of a "count the number of zeros in a row until one shows up". Not "count the number of ones in a bit-sequence". Hmm... so I guess that won't work either.

I'm brainstorming of course. Maybe I'll code this up when I get back home. I think I'm onto something.
First Strike +1/+1 and Indestructible.

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby Qaanol » Wed Apr 29, 2015 7:27 pm UTC

KnightExemplar wrote:Yeah, my methodology was to try to combine the speed advantage to your methodology and hotaru's "more complete" coverage of the double precision area. I recognize that your methodology is uniform and extremely quick.

But I recognize Hotaru's methodology as "more uniform" but slower, probably due to the conditional statement (the while loop) in his procedure. So despite the fact that it barely adds any computation, you can see it in the benchmark that its dramatically slower (3x or so).

Um, did you see the first line in my benchmark, the one referring to the code immediately above it in the same post as that benchmark?

The only reason Hotaru’s original implementation was slow, is that the 12 high bits were discarded instead of being used to adjust the exponent, so there were extra calls to arc4random (on average, 2 extra calls per result).

Using the high bits in place of the first up-to-12 such calls cut that down to an average of 1/2,048 extra arc4random calls per result.

The more-recent adjustment to pull 32 bits at a time reduces it further to something just barely above 1/4,096 extra calls per result.

However, Hotaru’s version of the 32-bits-at-a-time approach ended up slightly slower than mine because it had to increment (and test!) an extra variable even when fewer than 12 coin-flips were required. My version has a single test to short-circuit the entire “difficult” branch when it’s not needed.

I really like where this function is at now:

Code: Select all

#define TWELVE_HIGH_BITS_64 (0xfff0000000000000l)

void randomDoubleArray(double *v, const int n) {
    assert(sizeof(double) * CHAR_BIT == 64);
    
    
// Fill entire array with random bytes
    arc4random_buf(v, n * sizeof(double));
    
    for
(int i = 0; i < n; i++) {
        // Integer variables for bit-twiddling
        uint64_t *const p = (uint64_t *)(v+i);
        uint64_t h = *& TWELVE_HIGH_BITS_64;
        
        
*&= ~TWELVE_HIGH_BITS_64;  // Clear sign, exponent.
        uint64_t e = 1022;           // Exponent for [0.5, 1)
        
        
// Each exponent is as likely as all lower exponents,
        // so flip coins until tails, decrementing each time.
        // Use the bits of h for the first twelve coin flips.
        
        if 
(== TWELVE_HIGH_BITS_64) {
            // First 12 flips are heads
            bool proceeding = true;
            e = 1010;
            
            while 
(proceeding && e) {
                // Keep on flipping coins.
                uint32_t b = arc4random();
                proceeding = (== -1);
                while (& 1) { b >>= 1; e--; }
                e = (< 1024) ? e : 0;
            }
        } else {
            // Count the number of leading 1's in h.
            while (& (1l << 63)) { h <<= 1; e--; }
        }
        
        
*|= (<< 52);  // Apply e as the exponent
    }
}
 
wee free kings

User avatar
hotaru
Posts: 1045
Joined: Fri Apr 13, 2007 6:54 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby hotaru » Thu Apr 30, 2015 1:21 am UTC

Qaanol wrote:However, Hotaru’s version of the 32-bits-at-a-time approach ended up slightly slower than mine because it had to increment (and test!) an extra variable even when fewer than 12 coin-flips were required. My version has a single test to short-circuit the entire “difficult” branch when it’s not needed.

actually, I think using __builtin_ctzl() to count 0s instead of 1s is faster, the or at the end can be eliminated, and the and to clear the 12 high bits isn't needed if they're already 0, so:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i;
    uint32_t r = *>> 52;
    if(r) { *&= (1ul << 52) - 1; *|= 1022ul << 52; }
    else
    
{ *|= 1010ul << 52;
      for(int j = 31; j; --j)
      { if((= arc4random())) goto done;
        *-= 32ul << 52; }
      r |= 1 << 18; }
    done:
    *-= (uint64_t)__builtin_ctzl(r) << 52; }} 

Code: Select all

factorial product enumFromTo 1
isPrime n 
factorial (1) `mod== 1

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby Qaanol » Thu Apr 30, 2015 5:05 am UTC

hotaru wrote:actually, I think using __builtin_ctzl() to count 0s instead of 1s is faster, the or at the end can be eliminated, and the and to clear the 12 high bits isn't needed if they're already 0, so:

Code: Select all

void randomDoubleArray(double *v, const int n)
{
 arc4random_buf(v, n * sizeof(double));
  for(int i = 0; i < n; ++i)
  { uint64_t *const p = (uint64_t *)+ i;
    uint32_t r = *>> 52;
    if(r) { *&= (1ul << 52) - 1; *|= 1022ul << 52; }
    else
    
{ *|= 1010ul << 52;
      for(int j = 31; j; --j)
      { if((= arc4random())) goto done;
        *-= 32ul << 52; }
      r |= 1 << 18; }
    done:
    *-= (uint64_t)__builtin_ctzl(r) << 52; }} 

Very nice!

Seems to have everything one could ask for: magic numbers, floating point bit hacks, poorly documented built-in functions, backwards loops, even a GOTO statement! What’s not to love?

Heck, even the underlying “ARC4” random number generator itself was originally a trade secret that got leaked some 20 years ago!

Anyway, I’m seeing run-times around 0.308 seconds for n=10,000,000, so basically as fast as my original code that missed lots of values.

Do we have a way to guarantee that ctz will be available though? It looks like on older architectures it falls back to bsr which is undefined for 0 input. Not a huge deal since I’m only targeting new CPUs, but worth considering. (Edit: oh wait, code flow guarantees r will be nonzero at that line.)

Here’s essentially the same code (edit: I added an extra arc4random call before setting the sentinel bit because I think it’s necessary) with some comments, and formatted in my preferred style ;-)

Code: Select all

void randomDoubleArray(double *v, const int n) {
    // Fill entire array with random bytes
    arc4random_buf(v, n * sizeof(double));
    
    for 
(int i = 0; i < n; ++i) {
        // Integer variables for bit-twiddling
        uint64_t *const p = (uint64_t *)+ i;
        uint32_t r = *>> 52;
        
        
// Use trailing bits of r to decrement
        // the exponent of v[i] once per zero.
        if (r) {
            *&= (1ul << 52) - 1;
            *|= 1022ul << 52;
        } else {
            // More than 11 consecutive zeros.
            *|= 1010ul << 52;
            for (int j = 31; j-->0;) {
                r = arc4random();
                if (r) goto done;
                *-= 32ul << 52;
            }
            r = arc4random();
            r |= 1 << 18;
        }
        
    done
: *-= (uint64_t)__builtin_ctzl(r) << 52;
    }
}
 


I also wrote a “demagicked” version with a bunch of #define statements, but it ended up being even less readable. I did succeed in getting rid of every numeric literal beside two 0’s in loop bounds and three 1’s in bitwise operations though.

I won’t post the whole thing, but here are the macros I used:

Code: Select all

#define uint_base uint64_t
#define uint_rand uint32_t
#define RAND_BITS 32
#define DBL_EXP_DIG (sizeof(double) * CHAR_BIT - DBL_MANT_DIG)
#define MANT_BITS (DBL_MANT_DIG - 1)     // 52
#define HIGH_BITS (DBL_EXP_DIG  + 1)     // 12
#define BASE_EXP  (DBL_MAX_EXP  - 2)     // 1022
#define ALT_EXP   (BASE_EXP - HIGH_BITS) // 1010


Well, I guess I’ll at least mention that the “31” in the inner loop is “ALT_EXP / RAND_BITS”, and the “18” for the sentinel bit is the same but with “%”.

Edit: On second thought, here’s the demagicked version. Makes it easier to switch over to floats, and (somewhat) more resilient against changes to the platform.
Spoiler:

Code: Select all

#define uint_base uint64_t
#define uint_rand uint32_t
#define RAND_BITS 32
#define DBL_EXP_DIG (sizeof(double) * CHAR_BIT - DBL_MANT_DIG)
#define MANT_BITS (DBL_MANT_DIG - 1)     // 52
#define HIGH_BITS (DBL_EXP_DIG  + 1)     // 12
#define BASE_EXP  (DBL_MAX_EXP  - 2)     // 1022
#define ALT_EXP   (BASE_EXP - HIGH_BITS) // 1010

void randomDoubleArray(double *v, const int n) {
    // Fill entire array with random bytes
    arc4random_buf(v, n * sizeof(double));
    
    for 
(int i = 0; i < n; ++i) {
        // Integer variables for bit-twiddling
        uint_base *const p = (uint_base *)+ i;
        uint_rand r = *>> MANT_BITS;
        
        
// Use trailing bits of r to decrement
        // the exponent of v[i] once per zero.
        if (r) {
            *&= ((uint_base)<< MANT_BITS) - 1;
            *|= (uint_base)BASE_EXP << MANT_BITS;
        } else {
            // More than 11 consecutive zeros.
            *|= (uint_base)ALT_EXP << MANT_BITS;
            for (int j = ALT_EXP / RAND_BITS; j-->0;) {
                r = arc4random();
                if (r) goto done;
                *-= (uint_base)RAND_BITS << MANT_BITS;
            }
            r = arc4random();
            r |= (uint_rand)<< (ALT_EXP % RAND_BITS);
        }
        
    done
: *-= (uint_base)__builtin_ctzl(r) << MANT_BITS;
    }
}
 
wee free kings

User avatar
Jplus
Posts: 1721
Joined: Wed Apr 21, 2010 12:29 pm UTC
Location: Netherlands

Re: Uniformly random doubles from [0, 1) in C

Postby Jplus » Sat Aug 08, 2015 7:04 am UTC

Out of curiosity, what is wrong with <random>?
"There are only two hard problems in computer science: cache coherence, naming things, and off-by-one errors." (Phil Karlton and Leon Bambrick)

coding and xkcd combined

(Julian/Julian's)

User avatar
Qaanol
The Cheshirest Catamount
Posts: 3069
Joined: Sat May 09, 2009 11:55 pm UTC

Re: Uniformly random doubles from [0, 1) in C

Postby Qaanol » Sat Aug 08, 2015 10:53 pm UTC

The C++ part.
wee free kings


Return to “Coding”

Who is online

Users browsing this forum: No registered users and 7 guests