Can I use exp() and rand()?
I left the remaining as an exercise for the reader.
Regards,
Non-uniform random numbers
Ok... here goes... this is a modification to my general quadratic residue generator. Since everyone else is using rand() (which, btw, is a 3rd party library) I'm using GMP (mostly because I didn't feel like scaling my numbers down, it could be done, but there's really no reason to)
Output-
[jperalta@testarossa ~/src]$ ./biased.o 10 0
78967482
3532507987
769619069
1657007137
1493601390
318111018
1070407249
3036123774
2652192994
3264996435
[jperalta@testarossa ~/src]$ ./biased.o 10 1
0
4
12
0
6
0
0
2
2
2
[jperalta@testarossa ~/src]$ ./biased.o 10 2
0
0
0
4
0
0
2
2
2
2
[jperalta@testarossa ~/src]$ ./biased.o 10 3
2
2
0
0
0
0
0
2
0
0
QRG.H-#include <gmp.h>class QRG{ private: mpz_t n; mpz_t seed; public: QRG(); unsigned long GetNext(); unsigned long GetNext(int); unsigned int GetNextBit(); ~QRG();};QRG.CC-#include <gmp.h>#include <stdlib.h>#include "qrg.h"QRG::QRG(){ mpz_t p, q; mpz_init(p); mpz_init(q); mpz_init(n); mpz_init(seed); mpz_set_ui(p,2); mpz_pow_ui(p,p,127); mpz_sub_ui(p,p,1); // p is the Mersenne prime 2^127-1 mpz_set_ui(q,2); mpz_pow_ui(q,q,521); mpz_sub_ui(q,q,1); // q is the Mersenne prime 2^521-1 mpz_mul(n,p,q); srandomdev(); // Seed the system random number generator from /dev/random mpz_set_ui(seed,random()); // Use the system random number generator to generate our initial n mpz_powm_ui(seed,seed,2,n); // Calculate the initial seed value mpz_clear(p); // Delete p mpz_clear(q); // and q}QRG::~QRG(){ mpz_clear(n); mpz_clear(seed);}// This function returns a random 32-bit integer// and then calculates and stores the next seed valueunsigned long QRG::GetNext(){ unsigned long output = 0; mpz_t r; mpz_init(r); for (int i = 0; i < 32; i++) { unsigned long mod = mpz_mod_ui(r,seed,2); output = output + mod; if (i < 31) output = output << 1; mpz_powm_ui(seed,seed,2,n); } mpz_clear(r); return output;}unsigned long QRG::GetNext(int bias){ unsigned long output = 0; mpz_t r; mpz_init(r); for (int i = 0; i < 32; i++) { bool exiting = false; unsigned long mod = mpz_mod_ui(r,seed,2); output = output + mod; if (i < 31) output = output << 1; mpz_powm_ui(seed,seed,2,n); for (int i = 0; i < bias; i++) { if (GetNextBit() == 1) exiting = true; } if (exiting == true) break; } mpz_clear(r); return output;}// This function returns a random 0 or 1// and then calculates and stores the next seed valueunsigned int QRG::GetNextBit(){ unsigned int output = 0; mpz_t r; mpz_init(r); output = mpz_mod_ui(r,seed,2); mpz_powm_ui(seed,seed,2,n); mpz_clear(r); return output;}TEST.CC-#include <iostream>#include "qrg.h"int main(int argc, char** argv){ int trials = atoi(argv[1]); int bias = atoi(argv[2]); QRG* qrg = new QRG(); for (int i = 0; i < trials; i++) std::cout << qrg->GetNext(bias) << std::endl; delete qrg;}Output-
[jperalta@testarossa ~/src]$ ./biased.o 10 0
78967482
3532507987
769619069
1657007137
1493601390
318111018
1070407249
3036123774
2652192994
3264996435
[jperalta@testarossa ~/src]$ ./biased.o 10 1
0
4
12
0
6
0
0
2
2
2
[jperalta@testarossa ~/src]$ ./biased.o 10 2
0
0
0
4
0
0
2
2
2
2
[jperalta@testarossa ~/src]$ ./biased.o 10 3
2
2
0
0
0
0
0
2
0
0
Quote:Original post by Conner McCloud
As I recall [and its been a few years], F isn't the distribution, its the integral of the distribution.
I wasn't sure of the right term either, but distribution function is the integral of the probability function:
http://mathworld.wolfram.com/DistributionFunction.html
So it is the distribution function, which btw was called F on my math courses while f was the probability function.
Okay, one of the early xorshift randomizers, which make up the predecessors of the Merseen Twister. Amazingly fast. Uniform distribution. ( http://www.jstatsoft.org/v08/i14/xorshift.pdf )
Now, I want to assert that somewhere exist a fast interger square root function. I'm going to call it int fastsqrt(int). I don't know an implementation for it, and I don't feel like spending 5 minutes in google.
So, the probability of 0 through 65535 is equal. The probability of 65536 through 131071 is less than that of 0 through 65535, but equal among it's range. Follows up to 4294901760 through 4294967295, each number in the range is equivalent chance, but the whole range's chance is near impossible.
Now, the follow up.
For the sample, nuRand(100), 99 will be near impossible, but 0 will be the most occurant. My solution falls within the range of the given parameters, but shows faults when you test it's upper extremities.
static unsigned int seeda;static unsigned int seedb;unsigned int xorranda(){ seed^=(seed<<1); seed^=(seed>>3); seed^=(seed<<10); return seed;}unsigned int xorrandb(){ seed^=(seed<<4); seed^=(seed>>7); seed^=(seed<<19); return seed;}void xorsrand( unsigned int p ){ seeda=seedb=(p>>16)|(p<<16);}Now, I want to assert that somewhere exist a fast interger square root function. I'm going to call it int fastsqrt(int). I don't know an implementation for it, and I don't feel like spending 5 minutes in google.
unsigned int biased_rand(){ return (((65536 - fastsqrt(xorranda()))<<16)+(xorrandb&65535));}So, the probability of 0 through 65535 is equal. The probability of 65536 through 131071 is less than that of 0 through 65535, but equal among it's range. Follows up to 4294901760 through 4294967295, each number in the range is equivalent chance, but the whole range's chance is near impossible.
Now, the follow up.
int nuRand(const int nUpper){ return (int) ((float)biased_rand() / 4294967296.0 ) * nUpper;}For the sample, nuRand(100), 99 will be near impossible, but 0 will be the most occurant. My solution falls within the range of the given parameters, but shows faults when you test it's upper extremities.
Quote:Original post by jperalta
Since everyone else is using rand() (which, btw, is a 3rd party library)...
No its not. Its just as much part of the language as + or & is.
FlowingOoze: Oh, right. Nevermind [grin]
pragmaFury: That's why you never you never use magic numbers. I copied that straight from my random number class, but switched the actual random number generator call to rand() for simplicity. But I didn't bother changing the to float calculation to account for the new value of RAND_MAX. Oops. Should anybody care, to get that code to work replace the 'return rand()*(1.0/4294967295.0);' in getFloat() with 'return rand()*(1.0/RAND_MAX)'. I didn't submit that to compete, just to illustrate the generation of normal distributions, so don't worry about retrying my solution, I won't feel slighted.
CM
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
