Velocity Reviews > a gift for the mortensens

# a gift for the mortensens

Keith Thompson
Guest
Posts: n/a

 01-13-2010
frank <(E-Mail Removed)> writes:
[big snip]
> I've seen now several versions of
> i = (j / (RAND_MAX + 1.0)) * N;
> What is happening to 1222567701 / (RAND_MAX + 1.0) * 26?
> Of what type is (RAND_MAX + 1.0) ?

RAND_MAX is of type int. 1.0 is of type double. The "usual
arithmetic conversions" are applied to the operands (C99 6.5.6p4).
These are described in 6.3.1.8. RAND_MAX is converted from int to
double, and the addition yields a double result.

Grab a copy of <http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf>
if you haven't already. Or any decent C reference should explain it.

Something else I just thought of: An alternative to
(RAND_MAX + 1.0)
would be
(double)(RAND_MAX + 1)
The version with 1.0 has two advantages that I can think of:
it avoids the need for a cast, and it avoids integer overflow if
RAND_MAX == INT_MAX.

--
Keith Thompson (The_Other_Keith) http://www.velocityreviews.com/forums/(E-Mail Removed) <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

frank
Guest
Posts: n/a

 01-13-2010
Keith Thompson wrote:
> frank <(E-Mail Removed)> writes:
> [big snip]
>> I've seen now several versions of
>> i = (j / (RAND_MAX + 1.0)) * N;
>> What is happening to 1222567701 / (RAND_MAX + 1.0) * 26?
>> Of what type is (RAND_MAX + 1.0) ?

>
> RAND_MAX is of type int. 1.0 is of type double. The "usual
> arithmetic conversions" are applied to the operands (C99 6.5.6p4).
> These are described in 6.3.1.8. RAND_MAX is converted from int to
> double, and the addition yields a double result.
>
> Grab a copy of <http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf>
> if you haven't already. Or any decent C reference should explain it.
>
> Something else I just thought of: An alternative to
> (RAND_MAX + 1.0)
> would be
> (double)(RAND_MAX + 1)
> The version with 1.0 has two advantages that I can think of:
> it avoids the need for a cast, and it avoids integer overflow if
> RAND_MAX == INT_MAX.
>

I'm not sure that does what you think it does:

dan@dan-desktop:~/source\$ gcc -std=c99 -Wall -Wextra mort7.c -o out; ./out
mort7.c: In function main:
mort7.c:35: warning: integer overflow in expression
mort7.c:22: warning: unused variable g
mort7.c:21: warning: unused variable k
mort7.c:21: warning: unused variable i
3908406666 = 00010111010010111100000101011010
-2147483648.0000006 =
11000001111000000000000000000000000000000000000000 00000000000000
dan@dan-desktop:~/source\$ cat mort7.c
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <limits.h>

#define STRING "%7d6 = %s\n"
#define E_TYPE int
#define F_TYPE double
#define STRING2 "%7f6 = %s\n"
#define N 26

typedef E_TYPE e_type;
typedef F_TYPE f_type;

void bitstr (char *str, const void *obj, size_t n);

int
main (void)
{
int i, j, k;
double g;
e_type e;
f_type f;
char ebits[CHAR_BIT * sizeof e + 1];
char fbits[CHAR_BIT * sizeof f + 1];

srand (time (NULL));
j = rand ();
e = j;
bitstr (ebits, &e, sizeof e);
printf (STRING, e, ebits);

f = (double)(RAND_MAX + 1);
bitstr (fbits, &f, sizeof f);
printf (STRING2, f, fbits);

return 0;
}

void
bitstr (char *str, const void *obj, size_t n)
{
const unsigned char *const byte = obj;

while (n-- != 0)
{
mask = ((unsigned char) -1 >> 1) + 1;

do
{
*str++ = (char) (mask & byte[n] ? '1' : '0');
}
}
*str = '\0';
}

// gcc -std=c99 -Wall -Wextra mort7.c -o out; ./out
dan@dan-desktop:~/source\$

I think we all agree that none of these values are to be negative.
--
frank

frank
Guest
Posts: n/a

 01-13-2010
frank wrote:
> Keith Thompson wrote:

>> Something else I just thought of: An alternative to
>> (RAND_MAX + 1.0)
>> would be
>> (double)(RAND_MAX + 1)
>> The version with 1.0 has two advantages that I can think of:
>> it avoids the need for a cast, and it avoids integer overflow if
>> RAND_MAX == INT_MAX.
>>

>
> I'm not sure that does what you think it does:

Yes you do. You state it above, but I didn't understand that when I posted.
--
frank

frank
Guest
Posts: n/a

 01-13-2010
Phil Carmody wrote:
> Keith Thompson <(E-Mail Removed)> writes:
>> Richard Heathfield <(E-Mail Removed)> writes:
>>> Barry Schwarz wrote:
>>>> On Mon, 11 Jan 2010 15:14:23 -0800, Keith Thompson <(E-Mail Removed)>
>>>> wrote:
>>> <snip>
>>>
>>>>>
>>>>> i = (int) ((double) rand () / ((double) RAND_MAX + 1) * N);
>>>>>
>>>>> the cast to int is unnecessary, since the result is being assigned to
>>>>> an int object. The other two casts are necessary and appropriate,
>>>>> since in their absence the int values wouldn't be converted to double.
>>>> Only the second cast to double is necessary.
>>> Not even the second cast is necessary.
>>>
>>> i = (rand() / (RAND_MAX + 1.0)) * N;

>> Nicely done!
>>
>> In cases like this, though, I think I'd argue that using casts or not
>> is largely a matter of taste. I find that I'm a little uncomfortable
>> with the way the constant 1.0 imposes its type on the rest of the
>> expression, bubbling up through multiple levels of the tree.
>>
>> And I suppose that's inconsistent with most of what I've said about
>> using casts where implicit conversions would do the same job. Oh,
>> well.

>
> After years of dithering between different techniques, I now almost
> always use Richard's. I think nothing says 'use floating point' more
> concisely than actually sticking floating point numbers in the
> expression. Try it - you might like it
>
> Phil

I think it's important that the 1.0 bubbles up. You know a machine can
always give you one represented as a double. Adding RAND_MAX is then
small potatoes.

I've looked at a few different ways to do it, and the less good ones
involve building the integer, which isn't a good idea with RAND_MAX, and
then casting/converting it to a floating point type.

I still haven't figured out why we're creating a double type with a
value near RAND_MAX, but overflow is not the way to get there. Can
someone say a few words about
integer / double * integer?
--
frank

frank
Guest
Posts: n/a

 01-14-2010
Keith Thompson wrote:

> If you want random lowercase letters, you can declare
>
> const char letters[] = "abcdefghijklmnopqrstuvwxyz"
>
> and index into the array with a random number in the range 0..25.

I've got something now that looks alright, and it seems to test well for
the 200,000 times I put it through its paces:

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

#define N 26

char random_char (void);

int
main (void)
{
char a;
int i, counter;

srand (time (NULL));

counter = 0;

for (i = 0; i < RAND_MAX / 1000; ++i)
{
a = random_char ();
// printf("a is %c\n", a);
++counter;
}
printf ("counter is %d\n", counter);

return 0;
}

char
random_char (void)
{
char c;
int i, j;
const char letters[] = "abcdefghijklmnopqrstuvwxyz";
j = rand ();
i = (j / (RAND_MAX + 1.0)) * N;
assert (0 <= i && i < N);
c = letters[i];
return c;
}

// gcc -std=c99 -Wall -Wextra f1.c -o out; ./out

What I can't tell here is whether I've written something hugely
expensive in terms of computation. Would this code perform well?
--
frank

Phil Carmody
Guest
Posts: n/a

 01-14-2010
frank <(E-Mail Removed)> writes:
> I still haven't figured out why we're creating a double type with a
> value near RAND_MAX, but overflow is not the way to get there. Can
> someone say a few words about
> integer / double * integer?

Treated as ((integer/double)*integer). (integer/double) is
performed as (double/double) to yield a double. And then
(double*integer) performed as (double*double) to yield a
double.

It's integer/integer you normally want to be very wary of,
if you're embedding it within a tree of expressions that
will eventually become a double.

Phil
--
Any true emperor never needs to wear clothes. -- Devany on r.a.s.f1

Phil Carmody
Guest
Posts: n/a

 01-14-2010
Ben Bacarisse <(E-Mail Removed)> writes:
> Phil Carmody <(E-Mail Removed)> writes:
>> Ben Bacarisse <(E-Mail Removed)> writes:
>>> frank <(E-Mail Removed)> writes:

> <snip>
>>>> i = (int) ((double) rand () / ((double) RAND_MAX + 1) * N);
>>>
>>> Here be dragons. As 64 bit integers get more and more common we are
>>> edging towards a time when this will routinely fail because
>>> (double)RAND_MAX + 1 can be equal to RAND_MAX if RAND_MAX is big
>>> enough. This does not happen with a 32 RNG and normal IEEE
>>> double-precision numbers, but if RAND_MAX is big enough (and a signed
>>> 64-bit int is big enough) the +1 has no effect on (double)RAND_MAX.

>>
>> I'd like to see which real-world value of RAND_MAX and which
>> rounding mode could make that happen.

>
> I don't know of real world systems with 64-bit RAND_MAX. My warning
> was for the future of this common idiom. Of course, if you switch
> rand() for a 64-bit RNG (like KISS recently posted by George
> Marsaglia) you get the same problem if you use the idiom above.
>
>> Typically, RAND_MAX will
>> be odd, and adding 1 will cause a cascade of carries and leave
>> a result which requires fewer bits of precision to represent
>> accurately.

>
> That doesn't help, I don't think. Using my gcc, it is true that a
> large RAND_MAX get rounded up when converted to double, but so do a
> few of the returned values from rand(). Adding 1 does not have the
> effect of making the result of the division < 1.0 because the +1 has
> no effect.

AH, I didn't realise you were looking at that end of the issue.
Indeed, that is a problem.

> The probability of getting a large rand() from a 64-bit generator is
> very low, but it will happen one day! Here is the worse case example
> on my system. 512 returns from rand() can cause problems. You can
> fix it with long double but since that may be no wider than double,
> that option can just mask the problem until the code moves to another
> machine.
>
> #include <stdio.h>
> #include <limits.h>
>
> #define RAND_MAX LLONG_MAX
> #define N 10
>
> long long int rand(void) { return RAND_MAX-511; }
>
> int main(void)
> {
> int n = (rand() / (RAND_MAX + 1.0)) * N;
> if (n == N)
> puts("Don't use n as an array index!");

Yeah, but if you're mucking around with FP anyway, you shouldn't
expect results to not be inconveniently rounded occasionally
unless you've done some very thorough numerical analysis. 754 may

It does appear that implementers are very slow to migrate to
larger RAND_MAX, or in fact do anything to make rand() less of
the almost-useless toy that it currently is, so hopefully this
will never be a problem.

Phil
--
Any true emperor never needs to wear clothes. -- Devany on r.a.s.f1

frank
Guest
Posts: n/a

 01-16-2010
On Jan 13, 3:13*pm, Keith Thompson <(E-Mail Removed)> wrote:
> frank <(E-Mail Removed)> writes:
>
> [big snip]
>
> > I've seen now several versions *of
> > i = (j / (RAND_MAX + 1.0)) * N;
> > What is happening to 1222567701 / (RAND_MAX + 1.0) * 26?
> > Of what type is (RAND_MAX + 1.0) ?

>
> RAND_MAX is of type int. *1.0 is of type double. *The "usual
> arithmetic conversions" are applied to the operands (C99 6.5.6p4).
> These are described in 6.3.1.8. *RAND_MAX is converted from int to
> double, and the addition yields a double result.
>
> Grab a copy of <http://www.open-std.org/JTC1/SC22/WG14/www/docs/n1256.pdf>
> if you haven't already. *Or any decent C reference should explain it.
>
> Something else I just thought of: An alternative to
> * * (RAND_MAX + 1.0)
> would be
> * * (double)(RAND_MAX + 1)
> The version with 1.0 has two advantages that I can think of:
> it avoids the need for a cast, and it avoids integer overflow if
> RAND_MAX == INT_MAX.
>

Otherwise, if both operands have signed integer types or both have
unsigned
integer types, the operand with the type of lesser integer conversion
rank is
converted to the type of the operand with greater rank.
Otherwise, if the operand that has unsigned integer type has rank
greater or
equal to the rank of the type of the other operand, then the operand
with
signed integer type is converted to the type of the operand with
unsigned
integer type.
Otherwise, if the type of the operand with signed integer type can
represent
all of the values of the type of the operand with unsigned integer
type, then
the operand with unsigned integer type is converted to the type of the
operand with signed integer type.
Otherwise, both operands are converted to the unsigned integer type
corresponding to the type of the operand with signed integer type.

Can someone show an example for these integer promotions?