Re: Suggestion: Change the specification for grand in glib

[Sverre Johansen <sj ifi uio no>]

Sebastian Wilhelmi <wilhelmi ira uka de> writes:

> 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.

> 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

but with this multiplicator g_rand_double may return 1.0. 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.


> Ah, so what your basically saying is, that your implementation of
> g_rand_double is setting all 52 bits of gdouble, not only 32 as mine and it
> can thus be used to circumvent the 'if (dist > 2^16)' part of
> g_rand_int_range.

Correct.

> 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;
}

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    |
+--------------+----------------+-----------------------------------+
|              |       	       	|      	       	       	       	    |
|	       |		|				    |
.	       .		.				    .
.      	       .       	       	.				    .


Sverre Hvammen Johansen