WIP: Use single call to dsfmt_gv_genrand_uint32() for small ranges

[mschauer]

This is just an idea regarding #5549, this way one looses the type stability in RandIntGen, but gains faster number generation, when length(r) <= typemax(Uint32) and only one call to dsfmt_gv_genrand_uint32 to obtain entropy is needed.

[JeffBezanson]

This is probably not going to be a good place to give up type stability.

[mschauer]

I could not see any time difference in repeated calls to rand(int64(1:10)), rand(int64(1:10), 10000) on the other hand is two times faster than before here. But still worth to make everything type stable if this idea is interesting at all.

[mschauer]

I have changed the procedure a bit and rebased. The overall performance of seperate calls is slightly better than before due to less calls to rand(Uint32) for small ranges:

Master:

julia> @timen 1000000 rand(uint(1:10^15));
elapsed time: 0.071746714 seconds (0 bytes allocated)

julia> @timen 1000000 rand(uint(1:10^5));
elapsed time: 0.067608864 seconds (0 bytes allocated)

PR:

julia> @timen 1000000 rand(uint(1:10^15));
elapsed time: 0.07122678 seconds (0 bytes allocated)

julia> @timen 1000000 rand(uint(1:10^5));
elapsed time: 0.057847697 seconds (0 bytes allocated)

Note that these are hard to compare, because of the rejection part of the algorithm.

[lindahua]

The indentation of this function seems to be a little bit messy. Would you please clean up following julia's four-space convention?

[mschauer]

Sorry, some hidden tabs removed.

[lindahua]

What's the performance of rand(1:1000, 10^6) before & after?

[mschauer]

I get

before> @time rand(1:1000, 10^6);
elapsed time: 0.038024456 seconds (8000160 bytes allocated)
after> @time rand(1:1000, 10^6);
elapsed time: 0.029658467 seconds (8000160 bytes allocated)

[lindahua]

Good.

+1 for merging this.

[mschauer]

I'll test on 32 bit.

[StefanKarpinski]

This is probably not going to be a good place to give up type stability.

You don't return a different type. The point is that you only need 32 bits of entropy to return an Int64 within the range of Int32. If the amount of entropy used by rand(1:n) only depends on n and not on the type of n then you can get the same sequence of random values on both 32- and 64-bit systems, even though the types of the values will be different. Being able to generate consistent sequences across platforms with different word sizes is important.

[StefanKarpinski]

These convert calls are unnecessary since x is a typed variable.

[StefanKarpinski]

@mschauer, can you clarify about the type stability? Does this actually introduce any?

[mschauer]

No, I was maybe too generous with rebasing, this does not apply anymore.

[JeffBezanson]

Only a very early version of this was type-unstable. No longer relevant. This seems like a good change.

[mschauer]

I found a bug in the constructors and need a bit more time testing this. A problem is not straight forward to test because the results only can differ if rejections happen, which is unfrequent. I will add some real tests too.

[mschauer]

I am puzzled how the new function be faster than the old one for big ranges...
it only happens for certain ranges, for example for 0x1:0x0000040000000000 (the red dot is below the black one right in the figure.)

After some trial I am thinking that either this happens because the closing modulo operation takes less time or or I am missing something crucial.

Timings for some different sized ranges of rand(1:t, 10^6). PR is red dots, master black circles. Timings really depend on the range size, as g.a + rem(x, g.k) takes different time in each case. Red and black lines are PR and master without the division operation (replaced by g.a + x).

PS: The small peak at 0x1:0x0000000080000000 is expected, as this is the range which just doesn't fit two times into typemax(Uint32).

[mschauer]

To test I replaced the random number generator with a pseudo"random" number generator that just counts 0x00000000, 0x00000001, etc and everything is fine. Once #5585 converges to something things will likely change and #5550 can be fixed, should one wait or do this in a separate pull request?

[mschauer]

I will take 5585 into account and submit a new pull request.
@JeffBezanson : On jb/range1change I noted

int32(2:3)
> 2:4

[mschauer]

Closed in favor of 5741