Re: The answer is "All of it." of course.
Jonathan de Boyne Pollard
>
> for (uint32 N = 0; N < Max; N++)
> num = Data[num];
>
Your program is buggy. It's not actually calculating the metric that it
is purporting to print out.
> uint32 Random = rand() * rand();
>
Not random enough the first time, was it? (-:
Peter Olcott
"Jonathan de Boyne Pollard"
<J.deBoynePoll...@NTLWorld.COM> wrote in message
news:IU.D20100327.T...@J.de.Boyne.Pollard.localhost...
A single invocation of rand() would not by itself provide
random enough coverage over size when size is much greater
than 32767.
Dr Malcolm McLean
> "Jonathan de Boyne Pollard"<J.deBoynePollard-newsgro...@NTLWorld.COM> wrote in message
is unlikely to be what you want.
rand() * (RAND_MAX + 1) + rand();
is the uniform distribution.
pete
What makes you think that (RAND_MAX + 1) is defined?
--
pete
Jasen Betts
I assumed you were doing that because you didn't want a uniform distribution.
(no primes above 32767, denser near the middle, etc...)
but it's worse than that rand()*rand() has (approximately) a 1 in 16384.25
chance of producing the value 0
so your attempt to improve rand() has actually made it worse.
this means your test program is likely to find a loop starting at 0
and going for a few tens of thousands of steps before repeating.
you could do this
uint32 Random = rand() * (RAND_MAX+1) + rand();
but /dev/urandom is likely to serve you better if you want more randomness.
(or dev random if you need truly random and are prepared to wait for it
to arrive).
long longrand()
{
static int fd = -1;
long rv;
if( fd < 0)
fd=open("/dev/urandom",O_RDONLY);
if( fd<0 || read(fd,&rv,sizeof rv) < sizeof rv)
{
if( fd>=0)close(fd);
fd=-1;
return -1;
}
return rv & 0x7fffffff;
}
--- news://freenews.netfront.net/ - complaints: ne...@netfront.net ---
Eric Sosman
> Dr Malcolm McLean wrote:
>> [...]
>> However rand() * rand() will give you a binomial distribution, which
>> is unlikely to be what you want.
>> rand() * (RAND_MAX + 1) + rand();
>> is the uniform distribution.
>
> What makes you think that (RAND_MAX + 1) is defined?
And what makes him think that rand()*rand() has a binomial
distribution, even if it doesn't overflow?
--
Eric Sosman
eso...@ieee-dot-org.invalid
Peter Olcott
"Jasen Betts" <ja...@xnet.co.nz> wrote in message
news:hokglb$ked$1...@reversiblemaps.ath.cx...
> On 2010-03-27, Peter Olcott <NoS...@OCR4Screen.com> wrote:
>>
>> "Jonathan de Boyne Pollard"
>><J.deBoynePoll...@NTLWorld.COM> wrote in
>>message
>> news:IU.D20100327.T...@J.de.Boyne.Pollard.localhost...
>>> >
>>>>
>>>> for (uint32 N = 0; N < Max; N++)
>>>> num = Data[num];
>>>>
>>> Your program is buggy. It's not actually calculating
>>> the
>>> metric that it is purporting to print out.
>>>
>>>> uint32 Random = rand() * rand();
>>>>
>>> Not random enough the first time, was it? (-:
>>>
>> RAND_MAX = 32767
>> A single invocation of rand() would not by itself provide
>> random enough coverage over size when size is much
>> greater
>> than 32767.
>
> I assumed you were doing that because you didn't want a
> uniform distribution.
> (no primes above 32767, denser near the middle, etc...)
No I did this so that spatial locality could not as easily
improve the cache hit ratio.
>
> but it's worse than that rand()*rand() has
> (approximately) a 1 in 16384.25
> chance of producing the value 0
>
> so your attempt to improve rand() has actually made it
> worse.
>
> this means your test program is likely to find a loop
> starting at 0
> and going for a few tens of thousands of steps before
> repeating.
>
> you could do this
>
> uint32 Random = rand() * (RAND_MAX+1) + rand();
OK, I may try this. the results don't look random enough
right now.
>
> but /dev/urandom is likely to serve you better if you want
> more randomness.
> (or dev random if you need truly random and are prepared
> to wait for it
> to arrive).
That may slow things down too much, if I want true random I
an get a hardware random number generator dongle. For my
purposes I only need the numbers to be very widely
dispersed. I could generated my own pseudo random numbers
based on a lookup table.
Peter Olcott
"Dr Malcolm McLean" <malcolm...@btinternet.com> wrote
in message
news:4f098cf5-f4b8-4fdf...@z3g2000yqz.googlegroups.com...
On 27 Mar, 06:44, "Peter Olcott" <NoS...@OCR4Screen.com>
wrote:
> "Jonathan de Boyne
> Pollard"<J.deBoynePollard-newsgro...@NTLWorld.COM> wrote
> in message
>
Dr Malcolm McLean
> What makes you think that (RAND_MAX + 1) is defined?
>
> And what makes him think that rand()*rand() has a
> binomial distribution, even if it doesn't overflow?
>
A web search turns up that there's no neat way of obtaining the
probability density function of a the product of several uniformly-
distributed random variables. Basically you have to go to log space
and treat is as a sum (which approximates the normal distribution as N
increases).
Jasen Betts
ints are long (32bit) on linux, and I (and possibly several others)
didn't notice the crossposting to comp.lang.c (where ints can be short)
Phred Phungus
[broadening x-post to sci.math]
I'm reading this in comp.lang.c, and I don't think we have the firepower
to determine this distribution, hence the crosspost.
Hello sci.math. If rand has a flat pdf on the interval [0, 32768], what
is the distribution of rand squared?
Thanks for your comment, and cheers,
--
fred
Ray Vickson
If U is uniformly distributed on interval [0,M], then U^2 will be
distributed on interval [0,M^2]. Find the distribution as follows: for
any x in (0,M^2) the cumulative distribution of U^2 is F(x) = Pr{U^2
<= x} = Pr{U <= sqrt(x)}; we don't have negative values of U, so we
don't consider the interval from -sqrt(x) to 0. Anyway, for uniform U
on [0,M], Pr{U <= y} = y/M for y in [0,M]. For y = sqrt(x) this gives
Pr{U^2 <= x} = sqrt(x)/M [note: the right-hand-side = 1 when x = M^2,
as it should]. Thus, the density function of U^2 is f(x) = dF(x)/dx =
(1/2)/[M*sqrt(x)].
Note: this assumes continuously-distributed uniform U. For an *actual*
rand, U will be discrete. Furthermore, for a typical congruential
generator, the value 0 is never attained, and the distribution will be
slightly non-uniform over the other values 1, 2, ..., 32768.
R.G. Vickson
Phred Phungus
> On Mar 27, 10:46 pm, Phred Phungus <Ph...@example.invalid> wrote:
>> Dr Malcolm McLean wrote:
>>
>> [broadening x-post to sci.math]
>> Hello sci.math. If rand has a flat pdf on the interval [0, 32768], what
>> is the distribution of rand squared?
> If U is uniformly distributed on interval [0,M], then U^2 will be
> distributed on interval [0,M^2]. Find the distribution as follows: for
> any x in (0,M^2) the cumulative distribution of U^2 is F(x) = Pr{U^2
> <= x} = Pr{U <= sqrt(x)}; we don't have negative values of U, so we
> don't consider the interval from -sqrt(x) to 0. Anyway, for uniform U
> on [0,M], Pr{U <= y} = y/M for y in [0,M]. For y = sqrt(x) this gives
> Pr{U^2 <= x} = sqrt(x)/M [note: the right-hand-side = 1 when x = M^2,
> as it should]. Thus, the density function of U^2 is f(x) = dF(x)/dx =
> (1/2)/[M*sqrt(x)].
Ok. So, if one is to verify, then you're gonna heat up the
pseudorandoms and have M^2 bins. If the product of r = rand() and s =
rand() equals q, then bin number q, initially set to zero, increments.
q2) Does anyone have a statistical model for this event being
pseudorandom as opposed to being clunky in the lower bits?
>
> Note: this assumes continuously-distributed uniform U. For an *actual*
> rand, U will be discrete. Furthermore, for a typical congruential
> generator, the value 0 is never attained, and the distribution will be
> slightly non-uniform over the other values 1, 2, ..., 32768.
Thanks so much for your comment and your minor bounds error above like
the one I made. For clarity's sake, I should stipulate that the
interval is closed about zero and a value stipulated by the C
implementation. I'd like to choose this as [0, 32767] because then I
know I won't be in any trouble with M^2.
q3) How different is the analysis when it goes from continuous to discrete?
--
fred
Phil Carmody
> Dr Malcolm McLean wrote:
> [broadening x-post to sci.math]
>
>> On 27 Mar, 13:35, Eric Sosman <esos...@ieee-dot-org.invalid> wrote:
>>> What makes you think that (RAND_MAX + 1) is defined?
>>>
>>> And what makes him think that rand()*rand() has a
>>> binomial distribution, even if it doesn't overflow?
>>
>> I nodded.
>> A web search turns up that there's no neat way of obtaining the
>> probability density function of a the product of several uniformly-
>> distributed random variables.
What's not neat about the continuous case?
>> Basically you have to go to log space
>> and treat is as a sum (which approximates the normal distribution as N
>> increases).
And why do you introduce an N? There is no N.
> I'm reading this in comp.lang.c, and I don't think we have the
> firepower to determine this distribution, hence the crosspost.
>
> Hello sci.math. If rand has a flat pdf on the interval [0, 32768],
> what is the distribution of rand squared?
We're not squaring anything, we're multiplying 2 (purportedly)
independent values.
Phil
--
I find the easiest thing to do is to k/f myself and just troll away
-- David Melville on r.a.s.f1
Dr Malcolm McLean
wrote:
> Phred Phungus <Ph...@example.invalid> writes:
> > Dr Malcolm McLean wrote:
> > [broadening x-post to sci.math]
>
> >> On 27 Mar, 13:35, Eric Sosman <esos...@ieee-dot-org.invalid> wrote:
> >>> What makes you think that (RAND_MAX + 1) is defined?
>
> >>> And what makes him think that rand()*rand() has a
> >>> binomial distribution, even if it doesn't overflow?
>
> >> I nodded.
> >> A web search turns up that there's no neat way of obtaining the
> >> probability density function of a the product of several uniformly-
> >> distributed random variables.
>
> What's not neat about the continuous case?
>
>
Ray Vickson
wrote:
This is ambiguous in the original question. Having the line y =
rand()*rand() in a program will (probably) call rand twice and give
the product of (supposedly) independent uniform random variables. On
the other hand, the instruction y = rand()^2 will give the square of a
supposedly uniform random variable. So, both questions are valid.
Maybe the OP did not realize, or just forgot, the difference between
X*X and X^2 when X = rand() instruction. Interestingly, getting the
distribution of the square is easy, but not of the product; the latter
can be done but is not easy and is the subject of research papers.
R.G. Vickson
Liviu
On Mar 27, 10:46 pm, Phred Phungus <Ph...@example.invalid> wrote:
>> Dr Malcolm McLean wrote:
>>
>> [broadening x-post to sci.math]
>>
>> > On 27 Mar, 13:35, Eric Sosman <esos...@ieee-dot-org.invalid> wrote:
>> >> What makes you think that (RAND_MAX + 1) is defined?
>>
>> >> And what makes him think that rand()*rand() has a
>> >> binomial distribution, even if it doesn't overflow?
>>
>>> I nodded.
>>> A web search turns up that there's no neat way of obtaining the
>>> probability density function of a the product of several uniformly-
>>> distributed random variables. Basically you have to go to log space
>>> and treat is as a sum (which approximates the normal distribution
>>> as N increases).
>>
>> I'm reading this in comp.lang.c, and I don't think we have the
>> firepower to determine this distribution, hence the crosspost.
>>
>> Hello sci.math. If rand has a flat pdf on the interval [0, 32768],
>> what is the distribution of rand squared?
>
> distributed on interval [0,M^2]. Find the distribution as follows: for
> any x in (0,M^2) the cumulative distribution of U^2 is F(x) = Pr{U^2
> <= x} = Pr{U <= sqrt(x)}; we don't have negative values of U, so we
> don't consider the interval from -sqrt(x) to 0. Anyway, for uniform U
> on [0,M], Pr{U <= y} = y/M for y in [0,M]. For y = sqrt(x) this gives
> Pr{U^2 <= x} = sqrt(x)/M [note: the right-hand-side = 1 when x = M^2,
> as it should]. Thus, the density function of U^2 is f(x) = dF(x)/dx =
> (1/2)/[M*sqrt(x)].
>
> Note: this assumes continuously-distributed uniform U. For an *actual*
> rand, U will be discrete. Furthermore, for a typical congruential
> generator, the value 0 is never attained, and the distribution will be
> slightly non-uniform over the other values 1, 2, ..., 32768.
Following the background, I think the real underlying question was
more like "given independent discrete variables U1, U2
uniformly distributed on [0, M] what is the distribution of U1*U2".
Sort of the discrete counterpart to the continuous case of
http://mathworld.wolfram.com/UniformProductDistribution.html.
Except that the discrete case does not seem nearly as amenable and
smooth. For example, 0 has a probability of 2/M, primes in [2, M]
2/M^2, primes >M 0, the entire range [M(M-1)+1, M^2] also 0.
Don't know that there is a "nice" closed form for this distribution.
Liviu
David Schwartz
> That may slow things down too much, if I want true random I
> an get a hardware random number generator dongle. For my
> purposes I only need the numbers to be very widely
> dispersed. I could generated my own pseudo random numbers
> based on a lookup table.
On Linux, you can open /dev/urandom and read an unlimited supply of
random bits. So long as your distribution sets things up properly,
they should pass most statistical tests for randomness and should be
cryptographically strong.
DS
Phred Phungus
> This is ambiguous in the original question. Having the line y =
> rand()*rand() in a program will (probably) call rand twice and give
> the product of (supposedly) independent uniform random variables. On
> the other hand, the instruction y = rand()^2 will give the square of a
> supposedly uniform random variable. So, both questions are valid.
> Maybe the OP did not realize, or just forgot, the difference between
> X*X and X^2 when X = rand() instruction. Interestingly, getting the
> distribution of the square is easy, but not of the product; the latter
> can be done but is not easy and is the subject of research papers.
Ok, Ray, I was unclear about this myself. From my perspective, I'm not
OP, but from the perspective of sci.math, I introduced the original
question.
The interesting question I was thinking about is what is the
distribution of q:
r = rand();
s = rand();
q = r * s;
That I lapsed into a notation of talking about rand squared shows that I
lack many of the tools of writing about mathematics on the net. The
Germans are very good at it in de.sci.mathematik.
If there is an OP still paying attention, maybe he/she can specify what
he is looking for.
Cheers,
--
fred
Ersek, Laszlo
> David Schwartz <dav...@webmaster.com> writes:
> The problem here isn't so much the true randomness of the sequence as
> the distribution and the avoiding of accesses turning into a short
> loop.
(The local nntp server I used to use seems to have a problem with
forwarding my messages outwards. At least I was unable to find my messages
with google, and they appear absent also from news.eternal-september.org.
I went kind of crazy seeing this question discussed in multiple newsgroups
and nobody reflecting on my pertinent(?) message -- reposted below. I'm
switching to a different news client and the aforementioned news server,
in the hope that I'll be able to repost Sattolo's name to all relevant
newsgroups. Excuse the massive cross-posting. I set Followup-To to c.u.p.,
because that's where I saw the topic first.)
lacos
X-News: ludens comp.unix.programmer:167386
From: la...@ludens.elte.hu (Ersek, Laszlo)
Subject: Re: Can extra processing threads help in this case?
Date: 24 Mar 2010 02:49:32 +0100
Message-ID: <RfITPGAL3QUd@ludens>
In article <pe9k77-...@wilbur.25thandClement.com>, William Ahern <wil...@wilbur.25thandClement.com> writes:
> Peter Olcott <NoS...@ocr4screen.com> wrote:
>> "Eric Sosman" <eso...@ieee-dot-org.invalid> wrote in
>> message news:ho5tof$lon$1...@news.eternal-september.org...
> <snip>
>> > But if there's another CPU/core/strand/pipeline, it's possible that one
>> > processor's stall time could be put to productive use by another if
>> > there were multiple execution threads.
> <snip>
>> is there any possible way that this app is not memory bound that you can
>> provide a specific concrete example of?
>
> Your question was answered.
>
> You're hung up on your numbers and preconceived ideas. Your application
> could be BOTH memory bound AND able to benefit from multiple CPUs. But it's
> nearly impossible to guess without knowing at least the algorithm; more
> specifically, the code.
This thread has inspired me to write some code. I wrote it last night in
a sleep deprivation induced delirium. (I almost modified the upper limit
of LOG2_NEL from 31 to 32 before posting, but then decided to post the
version with 31 so that it works where "size_t" is a 32-bit unsigned
integer type with conversion rank at least that of "int".) I was fairly
sure about knowing what it measures. Now I'm fairly sure I was stupid.
So please tell me, I ask you kindly, what does the following code
measure? Thank you.
First some test results (inserting whitespace generously):
$ time -p ./jump_around 1
starting threads
14.5209 s/thread
real 16.43
user 16.37
sys 0.02
$ time -p ./jump_around 2
starting threads
15.571 s/thread
real 17.64
user 32.95
sys 0.05
$ time -p ./jump_around 4
starting threads
15.587 s/thread
real 33.25
user 64.16
sys 0.08
CPU (from <http://mysite.verizon.net/pchardwarelinks/current_cpus.htm>):
Athlon 64 X2 6000+ MMX 3DNow! SSE SSE2 SSE3 (Windsor)
(128-bit on-Die unbuffered DDR2 PC6400 mem controller, AMD-V)
(dual core)
940 pins, 3000MHz (200x15), (64-bit dual-pumped bus), 1.4v
2x 64KB data (2-way)
2x 64KB instruction (2-way)
2x 1MB on-Die unified L2 (16-way exclusive)
Here cometh the code. Please excuse the bugs, both the blatant and the
insidious.
Thanks,
lacos
/*
$Id: jump_around.c,v 1.4 2010/03/23 00:02:22 lacos Exp $
This program creates a single-cycle permutation of [0, 2**LOG2_NEL - 1] in a
static array with Sattolo's algorithm. Then it starts the specified number of
threads. Each sub-thread will traverse the permutation NROUND times. Each
thread starts at a different position in the loop (at offsets increasing
simply one by one in the buffer). Since the shuffle is driven by a uniform
distribution PRNG, the sub-threads should divide the loop more or less evenly
among themselves. The program measures and prints the user CPU time consumed
by a single sub-thread on average. Waiting for memory should be included in
this value.
Compile and link with eg.
$ gcc -ansi -pedantic -Wall -Wextra -o jump_around jump_around.c -l pthread
uint32_t and uint64_t are used only when their sizes are considered
significant for some reason.
*/
#define _XOPEN_SOURCE 500
#include <stdlib.h> /* size_t, strtol(), EXIT_FAILURE, jrand48() */
#include <inttypes.h> /* uint32_t, uint64_t */
#include <assert.h> /* assert() */
#include <stdio.h> /* fprintf(), fflush() */
#include <sys/resource.h> /* getrusage() */
#include <pthread.h> /* pthread_create(), pthread_join() */
/* Number of rounds a single sub-thread traverses the cycle. */
#define NROUND 10u
/* Highest number of sub-threads allowed on the command line. */
#define MAXTHR 4L
/*
Base two logarithm of the number of elements in the array (that is, log2 of
the cycle length). Make sure it is in [log2(MAXTHR), 31].
*/
#define LOG2_NEL 24
/* Number of elements; IOW, cycle length. */
#define NEL ((size_t)1 << LOG2_NEL)
/* 28 == LOG2_NEL will result in a 1GB array. */
static uint32_t buf[NEL];
/*
"*arg" points to a size_t object holding the start offset in the buffer.
*/
static void *
meander(void *arg)
{
unsigned round;
for (round = 0u; round < NROUND; ++round) {
size_t pos,
cnt;
pos = *(const size_t *)arg;
for (cnt = 0u; cnt < NEL; ++cnt) {
pos = buf[pos];
}
assert(pos == *(const size_t *)arg);
}
return 0;
}
int
main(int argc, char **argv)
{
long nthr;
/* Get number of sub-threads. */
{
char *endptr;
if (2 != argc || (nthr = strtol(argv[1], &endptr, 10)) < 1L
|| nthr > MAXTHR || '\0' != *endptr) {
(void)fprintf(stderr, "pass number of threads (1-%ld)\n", MAXTHR);
return EXIT_FAILURE;
}
}
/* Initialize array with Sattolo's algorithm. */
{
size_t cnt;
short unsigned xsubi[3] = { 1u, 2u, 3u };
/* Identity permutation. */
for (cnt = 0u; cnt < NEL; ++cnt) {
buf[cnt] = cnt;
}
/* Swaps. */
for (cnt = 0u; cnt < NEL - 1u; ++cnt) {
uint32_t save;
size_t rnd;
save = buf[cnt];
/*
Choose a pseudo-random number in [0, NEL - 2u - cnt] with uniform
distribution.
(uint32_t)jrand48(xsubi) covers [0, 2**32 - 1].
See also the constraint on LOG2_NEL.
*/
rnd = (uint64_t)(uint32_t)jrand48(xsubi) * ((NEL - 2u - cnt) + 1u) >> 32;
buf[cnt] = buf[cnt + 1u + rnd];
buf[cnt + 1u + rnd] = save;
}
}
(void)fprintf(stderr, "starting threads\n");
{
int ret;
struct rusage start, stop;
ret = getrusage(RUSAGE_SELF, &start);
assert(0 == ret);
{
long widx;
size_t startpos[MAXTHR];
pthread_t worker[MAXTHR];
for (widx = 0L; widx < nthr; ++widx) {
startpos[widx] = widx;
ret = pthread_create(worker + widx, 0, &meander, startpos + widx);
assert(0 == ret);
}
for (widx = 0L; widx < nthr; ++widx) {
ret = pthread_join(worker[widx], 0);
assert(0 == ret);
}
}
ret = getrusage(RUSAGE_SELF, &stop);
assert(0 == ret);
(void)fprintf(stdout, "%Lg s/thread\n",
(((long double)stop.ru_utime.tv_usec - start.ru_utime.tv_usec)
/ 1000000.0L + stop.ru_utime.tv_sec - start.ru_utime.tv_sec) / nthr);
if (EOF == fflush(stdout)) {
return EXIT_FAILURE;
}
}
return EXIT_SUCCESS;
}