Velocity Reviews > Re: Random 64 bit number in stdc

Re: Random 64 bit number in stdc

James Dow Allen
Guest
Posts: n/a

 08-19-2013
I've polished up my own random library. If you're
very masochistic^H^H^H^H^H^H^H^H^H^H^H eager to be a "beta tester",
e-mail me. PRNG isn't novel; it's based on Marsaglia etc.
Only "novelty" is avoiding unnecessary waste of random bits.

Here's a summary of the external specification:

/* Seed or reseed the random generator */
void rand_seed(int seed);

/* Return cnt random bits. 0 < cnt < 33 */
uint32_t rand_bits(unsigned int cnt);

/* Return random uniform 0 < x < 1 */
double rand_prob(void);

/* Return 1 with probability prob */
int rand_decide(double prob);

/* Return random uniform 0 <= x < range <= 65535 */
unsigned int rand_index(unsigned short range);

/* Return a uniform random variable from `minv' to `maxv' */
double rand_unif(double minv, double maxv);

/* Return Poisson random variable with mean and variance = ln(e_to_p) */
int rand_poisson(double e_to_p);

/* Return gaussian random variable */
double rand_gaussian(double meanv, double sdevv);

/* Set a random point on the surface of a N-dim hypersphere */
void rand_surfpoint(double coord[], double radius, int n);

/* Set a random point in interior of a N-dim hypersphere */
double rand_inball(double coord[], double radius, int n);

struct walkerent {
int x1, x2;
double prob1;
};

/* Given histogram wp[0...(n-1)], prepare for rand_walker() */
void rand_setwalker(struct walkerent *wp, unsigned short n);

/* Return 0 <= x < n according to set distribution */
int rand_walker(struct walkerent *wp, unsigned short n);

but abandoned that when I saw makefile complexity.
.... In a program which really only needs printf! ::whack::
(Am I the only one who prefers simplicity and standaloneness?)

James

James Dow Allen
Guest
Posts: n/a

 08-19-2013
Rosario1903 <(E-Mail Removed)> might have writ, in
news:(E-Mail Removed):

> On Mon, 19 Aug 2013 04:19:39 +0000 (UTC), James Dow Allen wrote:
>
>>I've polished up my own random library.

> why not
> void rand_seed(uint32_t* seed, uint32_t size);
> where seed is one array of unsigned of size: size

Your point is valid. I'll offer four excuses:

(1) Big seed may be needed for cryptography and casinos.
Small seed is good enough for simulations, which is my background.

(2) I did not find a good write-up on seeding for Marsaglia nor
Mersenne Twister. Correctness (i.e. copying faithfully so as not to
accidentally destroy randomness) is more important than functionality !

(3) My added value, if any, isn't the generator. It's the methods for
conserving bits, implementing Walker algorithm, etc.

(4) I wrote this for my own use and amusement. If it has any future at
all, that future will come *with the help* of other C programmers.

James

Eric Sosman
Guest
Posts: n/a

 08-19-2013
On 8/19/2013 12:19 AM, James Dow Allen wrote:
> I've polished up my own random library. If you're
> very masochistic^H^H^H^H^H^H^H^H^H^H^H eager to be a "beta tester",
> e-mail me. PRNG isn't novel; it's based on Marsaglia etc.
> Only "novelty" is avoiding unnecessary waste of random bits.
>
> Here's a summary of the external specification:
>
> /* Seed or reseed the random generator */
> void rand_seed(int seed);

That's not much entropy to seed a generator with lots
of internal state.

> /* Return random uniform 0 < x < 1 */
> double rand_prob(void);

That's peculiar: [0,1) is what I think most people expect.

> /* Return random uniform 0 <= x < range <= 65535 */
> unsigned int rand_index(unsigned short range);

The range restriction seems odd.

I'm not saying any of these design choices are "wrong,"
just that they're surprising.

--
Eric Sosman
http://www.velocityreviews.com/forums/(E-Mail Removed)d

James Dow Allen
Guest
Posts: n/a

 08-19-2013
On Monday, August 19, 2013 8:19:56 PM UTC+7, Eric Sosman wrote:
> On 8/19/2013 12:19 AM, James Dow Allen wrote:
> > void rand_seed(int seed);

>
> That's not much entropy to seed a generator with lots
> of internal state.

Answered in my previous post, more or less. Marsaglia's
or MTwister's numbers are very "random" just always one
of the same (up to 4 billion different) random sequences. ::whack::
If someone points me to a good write-up for seeding Marsaglia
*without worry that pathological case will degrade the PRNG*,
I will support it.

> > /* Return random uniform 0 < x < 1 */
> > double rand_prob(void);

> That's peculiar: [0,1) is what I think most people expect.

I want the mean value to be 0.5 exactly. (And I like symmetry.)
Anyway, the measure of (p==0 exactly) event would be small anyway
(2^-31).

> > /* Return random uniform 0 <= x < range <= 65535 */
> > unsigned int rand_index(unsigned short range);

> The range restriction seems odd.

Do you *ever* use random *non-power-of-two* indexes from
a larger range? A *non-factorable* range?
.... And rand_bits() is available if power of two.

I had larger range implemented but its value seemed outweighed
by added coding complexity. (I do like to think that if my
code has any virtues at all, *simplicity* is one of them.)

James

Eric Sosman
Guest
Posts: n/a

 08-19-2013
On 8/19/2013 9:38 AM, James Dow Allen wrote:
> On Monday, August 19, 2013 8:19:56 PM UTC+7, Eric Sosman wrote:
>> On 8/19/2013 12:19 AM, James Dow Allen wrote:
>>> [...]
>>> /* Return random uniform 0 <= x < range <= 65535 */
>>> unsigned int rand_index(unsigned short range);

>> The range restriction seems odd.

>
> Do you *ever* use random *non-power-of-two* indexes from
> a larger range? A *non-factorable* range?
> ... And rand_bits() is available if power of two.

"Indices," maybe not. "Values," certainly.

> I had larger range implemented but its value seemed outweighed
> by added coding complexity. (I do like to think that if my
> code has any virtues at all, *simplicity* is one of them.)

Tastes vary (that's all this discussion is about, really),
but a two-line rejection loop doesn't seem complex to me.

--
Eric Sosman
(E-Mail Removed)d

James Dow Allen
Guest
Posts: n/a

 08-19-2013
On Monday, August 19, 2013 8:51:49 PM UTC+7, Eric Sosman wrote:
> Tastes vary (that's all this discussion is about, really),
> but a two-line rejection loop doesn't seem complex to me.

Seven lines altogether, even without bit conservation,
since a completely new case is needed.
It was easier to just add them (as I've now done),
than to defend the indefensible.

James

Keith Thompson
Guest
Posts: n/a

 08-19-2013
James Dow Allen <(E-Mail Removed)> writes:
> I've polished up my own random library. If you're
> very masochistic^H^H^H^H^H^H^H^H^H^H^H eager to be a "beta tester",
> e-mail me. PRNG isn't novel; it's based on Marsaglia etc.
> Only "novelty" is avoiding unnecessary waste of random bits.
>
> Here's a summary of the external specification:
>
> /* Seed or reseed the random generator */
> void rand_seed(int seed);

If you want to have a specified number of bits in the seed, I suggest
using the typedefs from <stdint.h>, noting that the exact-width types
are not guaranteed to exist.

I note that the standard srand() takes an unsigned int seed argument,
which seems more convenient. Does yours work correctly with negative
seeds?

> /* Return cnt random bits. 0 < cnt < 33 */
> uint32_t rand_bits(unsigned int cnt);

I suggest spelling out "count".

[...]

--
Keith Thompson (The_Other_Keith) (E-Mail Removed) <http://www.ghoti.net/~kst>
Working, but not speaking, for JetHead Development, Inc.
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

Eric Sosman
Guest
Posts: n/a

 08-19-2013
On 8/19/2013 12:59 PM, James Dow Allen wrote:
> On Monday, August 19, 2013 8:51:49 PM UTC+7, Eric Sosman wrote:
>> Tastes vary (that's all this discussion is about, really),
>> but a two-line rejection loop doesn't seem complex to me.

>
> Seven lines altogether, even without bit conservation,
> since a completely new case is needed.
> It was easier to just add them (as I've now done),
> than to defend the indefensible.

Mine (based on a different PRNG, with a typedef and some
macros that shouldn't baffle anyone):

/**
* Returns a non-negative pseudo-random value strictly less
* than its argument, which must be in the range 1 through
* MSRANDMAX-MSRANDMIN+1, inclusive. Returned values are
* exactly equiprobable.
*/
msrand_t
msrandint(msrand_t limit)
{
msrand_t divisor, value;
divisor = (MSRANDMAX - MSRANDMIN + 1) / limit;
while ( (value = (msrand() - MSRANDMIN) / divisor) >= limit )
;
return value;
}

The loop could be micro-optimized a little, if anyone cared
(I never did, since even in the worst case the expected number
of divisions per call is slightly less than two).

--
Eric Sosman
(E-Mail Removed)d

Eric Sosman
Guest
Posts: n/a

 08-19-2013
On 8/19/2013 3:26 PM, Eric Sosman wrote:
> [...]
> msrand_t
> msrandint(msrand_t limit)
> {
> msrand_t divisor, value;
> divisor = (MSRANDMAX - MSRANDMIN + 1) / limit;
> while ( (value = (msrand() - MSRANDMIN) / divisor) >= limit )
> ;
> return value;
> }
>
> The loop could be micro-optimized a little, if anyone cared
> (I never did, since even in the worst case the expected number
> of divisions per call is slightly less than two).

Sorry: That should be "divisions per call *in the loop*,"
hence "slightly less than three all told."

--
Eric Sosman
(E-Mail Removed)d

James Dow Allen
Guest
Posts: n/a

 08-20-2013
On Tuesday, August 20, 2013 2:26:24 AM UTC+7, Eric Sosman wrote:

> > Seven lines altogether, even without bit conservation,

> ...
> while ( (value = (msrand() - MSRANDMIN) / divisor) >= limit )
> ;

I do NOT quibble with your fine code, but just seek to clarify
my "7 lines": In my older age, I'll often use {} when ;
could suffice. (Two of my 7 were together only "do {} ".
Similarly, I also spend two lines for the "if (range > OTHER) {}"

But much more importantly, your code will spend 31(?) random
bits to pick a value in 0...14 while (slightly less than) 4

James