IF the underlying statistical properties of the built in PRNG are good enough for your application (and they may or may not be), then you can get pseudo-random numbers of any size you want just by concatenating repeated calls.

my15bitRand = rand();
my30bitRand = rand()<<15 + rand();
my45bitRand = (rand()<<15 + rand())<<15 + rand();
my60bitRand = ((rand()<<15 + rand())<<15 + rand())<<15 + rand();

I was thinking of basically the same idea, but I think your idea is better as what is going on is clearer. Nice thanks!

 

Your approach doesn't produce anywhere close to the probability you want because the distribution of x*y is not uniform -- it's highly non-uniform.

Imagine that RAND_MAX was 99 (i.e., 100 possible values). The RAND_MAX*RAND_MAX would be 9801 and there would be 10,000 possible pairings of p1 and p2. Let's say that you wanted a probability of 0.1%. That would mean that your threshold would be 9.801, which means you would require p1*p2 to be 9 or less.

How many ways can this be achieved, out of the 10,000 possible outcomes? You are looking for there to only be about 10.

If p1 = 0, then it doesn't matter what p2 is, so there's 100 right off the bat.
If p1 = 1, then p2 can be {0,1,2,3,4,5,6,7,8,9}, so there's 10 more.
If p1 = 2, then p2 can be {0,1,2,3,4}, so another 5.
If p1 = 3, then p2 can be {0,1,2,3}, another 4.
If p1 = 4, then p2 can be {0,1,2}, another 3
If p1 = 4, 5, 6, 7, 8, or 9, then p2 can be {0,1}, another 12
If p1 > 9, then p2 can be {0}, which adds a final 90.

So the probability that p1*p2 < 9.801 is not 0.1%, but rather 2.24%

At the other end of the scale, let's say that you want a probability of 99%. There you want there to be 100 pairings whose product is at least 9702.99. So 9703 or higher. How many are there?

If p1 = 99, then p2 can be {98,99}
If p1 = 98, then p2 can be {99}

That's it. So instead of 100 pairs being above the threshold, there are only 3.

 

Your approach doesn't produce anywhere close to the probability you want because the distribution of x*y is not uniform -- it's highly non-uniform.

Imagine that RAND_MAX was 99 (i.e., 100 possible values). The RAND_MAX*RAND_MAX would be 9801 and there would be 10,000 possible pairings of p1 and p2. Let's say that you wanted a probability of 0.1%. That would mean that your threshold would be 9.801, which means you would require p1*p2 to be 9 or less.

How many ways can this be achieved, out of the 10,000 possible outcomes? You are looking for there to only be about 10.

If p1 = 0, then it doesn't matter what p2 is, so there's 100 right off the bat.
If p1 = 1, then p2 can be {0,1,2,3,4,5,6,7,8,9}, so there's 10 more.
If p1 = 2, then p2 can be {0,1,2,3,4}, so another 5.
If p1 = 3, then p2 can be {0,1,2,3}, another 4.
If p1 = 4, then p2 can be {0,1,2}, another 3
If p1 = 4, 5, 6, 7, 8, or 9, then p2 can be {0,1}, another 12
If p1 > 9, then p2 can be {0}, which adds a final 90.

So the probability that p1*p2 < 9.801 is not 0.1%, but rather 2.24%

At the other end of the scale, let's say that you want a probability of 99%. There you want there to be 100 pairings whose product is at least 9702.99. So 9703 or higher. How many are there?

If p1 = 99, then p2 can be {98,99}
If p1 = 98, then p2 can be {99}

That's it. So instead of 100 pairs being above the threshold, there are only 3.

You are correct, I am testing it at the moment, here is the modified one:

 

Your approach doesn't produce anywhere close to the probability you want because the distribution of x*y is not uniform -- it's highly non-uniform.

Imagine that RAND_MAX was 99 (i.e., 100 possible values). The RAND_MAX*RAND_MAX would be 9801 and there would be 10,000 possible pairings of p1 and p2. Let's say that you wanted a probability of 0.1%. That would mean that your threshold would be 9.801, which means you would require p1*p2 to be 9 or less.

How many ways can this be achieved, out of the 10,000 possible outcomes? You are looking for there to only be about 10.

If p1 = 0, then it doesn't matter what p2 is, so there's 100 right off the bat.
If p1 = 1, then p2 can be {0,1,2,3,4,5,6,7,8,9}, so there's 10 more.
If p1 = 2, then p2 can be {0,1,2,3,4}, so another 5.
If p1 = 3, then p2 can be {0,1,2,3}, another 4.
If p1 = 4, then p2 can be {0,1,2}, another 3
If p1 = 4, 5, 6, 7, 8, or 9, then p2 can be {0,1}, another 12
If p1 > 9, then p2 can be {0}, which adds a final 90.

So the probability that p1*p2 < 9.801 is not 0.1%, but rather 2.24%

At the other end of the scale, let's say that you want a probability of 99%. There you want there to be 100 pairings whose product is at least 9702.99. So 9703 or higher. How many are there?

If p1 = 99, then p2 can be {98,99}
If p1 = 98, then p2 can be {99}

That's it. So instead of 100 pairs being above the threshold, there are only 3.

True, that became apparent to me when I considered the first random number being 0 would force the entire expression to zero being it is just a multiplication, needs to be like you said a shift followed by an addition.

 

One slight tweak, which you will likely never notice the difference from.

If you want a probability of success of 1.0, you won't quite get it because if both of your rand() call return RAND_MAX, your p will be equal to probability range and you will return 0.

You could change the relation to <= instead of <. But then if your request probability is 0.0, you will still have a chance of getting a 1 returned.

You actually do just want Range = (1 << 30) combined with using strict inequality.

 

True, that became apparent to me when I considered the first random number being 0 would force the entire expression to zero being it is just a multiplication, needs to be like you said a shift followed by an addition.

There is still something wrong with the following though, because my testing is coming out with 35 times to many successes on a test bench:

 

What probabilityOfSuccess were you using?

 

One slight tweak, which you will likely never notice the difference from.

If you want a probability of success of 1.0, you won't quite get it because if both of your rand() call return RAND_MAX, your p will be equal to probability range and you will return 0.

You could change the relation to <= instead of <. But then if your request probability is 0.0, you will still have a chance of getting a 1 returned.

You actually do just want Range = (1 << 30) combined with using strict inequality.

If I can figure out how to attach code to the thread, I could send you the test code. I am doing something wrong.

 

What probabilityOfSuccess were you using?

0.01 (i.e. 1%)

 

Here is the test code I was testing with:


unsigned int GetSeedFromSystemTime();

int GenerateLowProbabilityEvent(double probabilityOfSuccess);

int main()
{
const double lowProbabilityEvent = 0.01;
const double trials = 10000;
srand(GetSeedFromSystemTime());
int successCounter = 0;

for(int x=0;x<trials;x++)
{
successCounter += GenerateLowProbabilityEvent(lowProbabilityEvent);
}

printf("Percentage is %f and should be %f\n", 100.0*successCounter / trials,100.0*lowProbabilityEvent);
}

unsigned int GetSeedFromSystemTime()
{
SYSTEMTIME time;
::GetSystemTime(&time);
return time.wMilliseconds;
}

// Return 1 for success and 0 otherwise where probabilityOfSuccess is a number from 0.0 to 1.0
int GenerateLowProbabilityEvent(double probabilityOfSuccess)
{
const unsigned int Range = 1 << 30;
unsigned int probabilityRange = Range * probabilityOfSuccess;
unsigned int p = rand()<<15 + rand();
return p < probabilityRange ? 1 : 0;
}

 

Be sure to take into account the order of precedence of << verses addition. I wasn't considering that when I made by suggestion above and, off the top of my head, I'm not sure what the order is, but I highly suspect that addition is higher precedence. So you need to use:

(rand()<<15) + rand()

 

What probabilityOfSuccess were you using?

Here is the output from the code, that I attached, with 35 times to many successes? You would not think this should be having issues?

 

Be sure to take into account the order of precedence of << verses addition. I wasn't considering that when I made by suggestion above and, off the top of my head, I'm not sure what the order is, but I highly suspect that addition is higher precedence. So you need to use:

(rand()<<15) + rand()

Nice catch! You are correct! That is the problem here is the output after that fix. I always though << should have precedence over addition, apparently not.

 

Here is the output from the code, that I attached, with 35 times to many successes? You would not think this should be having issues?

View attachment 285331

How many decimal places (or bits) can your test code handle?

 

Here's my code:

#include <stdio.h>
#include <stdlib.h>

#define PMAX (1<<30)
#define TRIALS (100000000)

int lowP(double pSuccess)
{
    return ((rand()<<15) + rand()) < (PMAX * pSuccess);
}

int main(void)
{
    for (double pSuccess = 0.0000001; pSuccess <= 1.0; pSuccess *= 10)
    {
        int hits = 0;
        for (int i = 0; i < TRIALS; i++)
            hits += lowP(pSuccess);
        printf("Target: %10.6f%% Result %10.6f%%\n", 100.0*pSuccess, 100.0*( (double) hits / (double) TRIALS));
    }
    
    return 0;
}

Here's what this produces:

Target:   0.000010% Result   0.000009%
Target:   0.000100% Result   0.000111%
Target:   0.001000% Result   0.001065%
Target:   0.010000% Result   0.010049%
Target:   0.100000% Result   0.100094%
Target:   1.000000% Result   0.998634%
Target:  10.000000% Result  10.001597%
Target: 100.000000% Result 100.000000%

 

How many decimal places (or bits) can your test code handle?

An int on this system is 4 bytes long.

 

Be sure to take into account the order of precedence of << verses addition. I wasn't considering that when I made by suggestion above and, off the top of my head, I'm not sure what the order is, but I highly suspect that addition is higher precedence. So you need to use:

(rand()<<15) + rand()

I do a lot of bit shifting and adding, but I think the reason I never an into that before is because I typically do it this way:

x=y<<16 | z;

as to:

x = y<<16 + z;

Thanks for your help.

 

Nice catch! You are correct! That is the problem here is the output after that fix. I always though << should have precedence over addition, apparently not.

I'm pretty sure that the reason that it is lower precedence is to make it easier to do shifts using calculated amounts, but most people only use it to do hardcoded fixed shifts.

 

Radioactive decay and detector noise are very good HRNG (Hardware Random Number Generator) but my absolute favorite HRNG is this one.

 

Using the noise on an unconnected analog input to seed the random number generator also works well.