New random_unit_interval implementation

Re-wrote the random_unit_interval function to generate doubles with 48 bits of random mantissa (3.5e-15), as opposed to the current 32-bit version (2.3e-10).

This implementation uses rand_r (to blend the initial state bytes), and erand48 to generate the result over the initial state bytes. The implementation of both functions is inlined to avoid relatively expensive function calls.

This change passes both the testRandom and testRandomSpacing tests.

@matthieu, @folkert, testRandomSpacing hasn't finished quite yet, but the results look decent:

[02973.3] main: Categories           |  1e-05  1e-06  1e-07  1e-08  1e-09  1e-10
[02973.3] main: total =  52344913920 | 522609  52521   5165    517     50      7
[02973.3] main: expected             | 523449  52344   5234    523     52      5
[02973.3] main: Difference           |    840    177     69      6      2      2
[02973.3] main: 5 sigma difference   |   3615   1140    360    110     35     10

Speed-wise, the new function is ~1.5x slower than the original (configured with --enable-optimization).

What do you think?

Cheers, Pedro

This looks good! It indeed seems to pass testRandomSpacing and it is less files than the one using Random123 and speed-wise it is not much different from the method using Random123.

added 1 commit

• ddc2a6fd - remove Random123 sub-directory, not needed in this branch.

Compare with previous version

@folkert, cool! If we're going with this branch, can you update the descriptions in doc/RTD/? Cheers!

Thanks both!

@nnrw56 do we worry about the re-implementation of the rand_r function that may not work on other architectures? e.g. ARM?

@folkert could you try this on one of your isolated galaxy examples to make sure we haven't broken the higher level things?

mentioned in merge request !862 (closed)

@matthieu, the code was taken from the glibc implementation, so it's essentially the same code we get when we call the function directly. The implementations also don't have anything platform-specific about them. For safety/sanity, I could replace the short int and unsigned int types with int16_t and uint32_t respectively, just to be sure we don't break anywhere.

@nnrw56, yes I can update the RTD.

@matthieu, sure I will test it on a isolated galaxy example.

@nnrw56 I was just wondering whether glibc would use a different implementation on a different architecture.

IIRC each function may have an architecture-specific implementation for performance reasons, but they should produce identical results to the generic implementation, which is what we're using here.

@nnrw56, is it easy in this RNG to extend the number of different random number types? Or should we test a few different numbers using testRandom and if they work add them?

@folkert, I don't understand the question :) What do you mean by different number types?

We want to have different random numbers for the same particle and the same time step, currently we do that by changing the number type like for example random_number_stellar_feedback and random_number_star_formation, we know that we need more types than are currently available, so we also need to extend this by a few more types and potentially add a few spare ones such that we do not need to think about it in the future.

We have multiple physical processes in the simulations that require some level of stochasticity. At a given time-step (ti_current) and for a given particle (a given ID), we do need different uncorrelated numbers. Essentially, using additional salt, which is what the last argument (the enum) to the function call is.

@folkert beat me to it with a better answer...

Got it :) So you don't test changes in the last argument in testRandom?

We do but that is only for the currently available numbers (or maybe only the first 4).

Document that this is taken straight from glibc with a link to the original?

Done.

changed this line in version 3 of the diff

Basically, Folkert's question is whether we can simply extend the enum with any value or whether we should be careful about what we pick.