Re: Suggestion: Change the specification for grand in glib

[Sebastian Wilhelmi <wilhelmi ira uka de>]

Hi Sverre,

> > now we have it and I would like to keep it. I hope the rational for this is
> > now more clear.
> 
> Yes it is, but I do not like it. If the specification is not changed,
> the parameter max should be renamed.

We could call it 'begin' and 'end' just like the STL parameters, if that would
help. but calling them lower_closed_bound and upper_open_bound seems a bit
odd. 

> > That must be due to the bad multiplicator. If choosen carefully it should
> > work.
> 
> Lets name G_RAND_DOUBLE_TRANSFORM T.  The correct multiplicator to use
> if only one call to g_rand_int is performed is:
> 
> T'= T + T*T

This of course depends on, what G_RAND_DOUBLE_TRANSFORM is. I would like it to
be 1/2^32, such that we get the open intervall and in that case you're
obviously right.

> but with this multiplicator g_rand_double may return 1.0.

Yes, but only with bad rounding luck, I suppose.

> If we round down or to zero we will get the range [0..1), otherwise [0..1].
> Lets look at what is happening.
> 
> T is the following constant written in hexadecimal notation in the
> style 0.dddddddd:
> 
> T: 0.00000001
> 
> When we multiply T with 2^32-1 we get:
> 
> 0.FFFFFFFF
> 
> If we instead use T':
> 
> T': 0.0000000100000001
> 
> and multiply this with 2^32-1 we get:
> 
> 0.FFFFFFFFFFFFFFFF
> 
> Only the first 53 ones will be represented unless it is rounded up to
> 1.0.
> 
> It should be obvious that T' is the perfect multiplicator if we only
> want to make one call to g_rand_int.

Agreed.

> > Actually that seems like a good idea, even if it requires
> > g_rand_double to fetch to 32 bit random numbers. But your implementation
> > really looks weird to me (That can be my fault, though). Can you please
> > explain, what your are doing there.
> 
> Lets start the explanation by forgetting the if test in the program.
> I continue with explaining the for loop.
> 
> (Please see the attached table at the bottom of this mail.)
> 
> The normal case is that initially r1 != 0 and transform is assigned
> 0.00000001. r1 == 0 initially means that g_rand_double will return a
> number below 0.00000001. If we don't perform an extra call to
> g_rand_int g_rand_double will return only 32 bits. The loop will
> perform an extra call to g_rand_int and assign transform
> 0.0000000000000001.
> 
> r1 == 0 initially and r1 == 0 after the first extra call to g_rand_int
> in the for loop means that g_rand_double will return a number below
> 0.0000000000000001 . If we don't perform an extra call to g_rand_int
> also in this case g_rand_double will again return only 32 bits.
> 
> After the loop if r1 < 1<<20 lesser than 53 bits will be returned, and
> we therfore perform an extra call to g_rand_int.
>
> If we at all should test for the case r1 == 0 there is no extra cost
> making this into a loop. The overhead concerned the for-loop and the
> if-test is small compared to the overhead in g_rand_int.
> 
> But for the use needed in g_rand_int_range the following simplified
> code is perfectly good. I don't think that this code will perform
> significantly faster so I suggest that my original suggestion should
> be used:
> 
> gdouble
> g_rand_double (GRand* rand)
> {
>   return (g_rand_int (rand) + g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM)
>          * G_RAND_DOUBLE_TRANSFORM;
> }

This is actually, what I would suggest and what I will use, see below.

> Original suggestion:
> 
> gdouble
> g_rand_double (GRand* rand)
> {
>   gdouble transform;
>   guint32 r1 = g_rand_int (rand);
>   gdouble r2 = g_rand_int (rand);
> 
>   for (transform = G_RAND_DOUBLE_TRANSFORM; r1 == 0;
>        transform *= G_RAND_DOUBLE_TRANSFORM)
>     r1 = g_rand_int (rand);
> 
>   if (r1 < 1<<20)
>     r2 += g_rand_int (rand) * G_RAND_DOUBLE_TRANSFORM;
> 
>   return (r1 + r2 * G_RAND_DOUBLE_TRANSFORM) * transform;
> }
> 
> Some expressions:
> -----------------
> 
> r1: aaaaaaaa
> 
> r2: bbbbbbbb
> 
> r1 + r2 * G_RAND_DOUBLE_TRANSFORM:
> 
>                    aaaaaaaa.bbbbbbbbcccccccc
>                        ^        ^      ^
>                        |        |      |
>                       r1       r2    r1<1<<20
> 
> +--------------+----------------+-----------------------------------+
> | Number of    |                |                                   |
> | iterations   |    transform   |(r1 + r2 * G_RAND_DOUBLE_TRANSFORM)|
> | in the loop  |                |     * transform                   |
> +--------------+----------------+-----------------------------------+
> |              |                |      0.aaaaaaaabbbbbbbb           |
> |              |                |             ^       ^             |
> |      0       |   0.00000001   |             |       |             |
> |              |                |            r1      r2             |
> +--------------+----------------------------------------------------+
> |              |                |   0.00000000aaaaaaaabbbbbbbb      |
> |              |   0.00000000   |                 ^       ^         |
> |      1       |     00000001   |                 |       |         |
> |              |                |                r1      r2         |
> +--------------+----------------+-----------------------------------+
> |              |   0.00000000   |0.0000000000000000aaaaaaaabbbbbbbb |
> |              |     00000000   |                      ^       ^    |
> |      2       |     00000001   |                      |       |    |
> |              |                |                     r1      r2    |
> +--------------+----------------+-----------------------------------+
> |              |                |                                   |
> |              |                |                                   |
> .              .                .                                   .
> .              .                .                                   .

Come on, this is silly. Firstly case 0 will happen in 4294967295 of 4294967296
cases not to talk about the cases 1 and 2 (I think, it is more probable, that
the world collapses, than to see the case 2). And if it does, nobody needs the
higher accuracy of the number. Actually it could harm, because it provides a
better granularity, than the user might expect. 

I mean, if you return an integer in [0,32^2-1] you are not returning less
information, if the first bits are zero. Even if double allows you to return
more bits, if the first ones after the point are zero, it doesn't seem to make
any sense to really do that. 

Unless you provide some good examples, where that is necessary or helps, I'm
very unlikely to add this. 

But nontheless thanks for bringing up the issue with the bad multiplicator and
the idea to set all 52 bits after the point for double.

I'm going to check in some changes. Please have a look.

Bye,
Sebastian
-- 
Sebastian Wilhelmi
mailto:wilhelmi ira uka de
http://goethe.ira.uka.de/~wilhelmi