juliarandom / randomnumbers.jl Goto Github PK

Hi,
is it possible to set one of these RNGs as the system default somehow? I'd like to use these myself, but much of what I use / write is library code, so I would like to not have to pass r to rand every time.
Thanks!

Tests fail for RandomNumbers on OSX (10.14.3) on Julia 1.10.

Procedure:

1. ] up
2. test RandomNumbers

Results:
Testing RandomNumbers
Resolving package versions...
Status /var/folders/26/dyl9dp2157z1g3sck948l6s40000gn/T/tmpK4hg52/Manifest.toml
[e6cf234a] RandomNumbers v1.1.0
[2a0f44e3] Base64 [@stdlib/Base64]
[ade2ca70] Dates [@stdlib/Dates]
[8ba89e20] Distributed [@stdlib/Distributed]
[b77e0a4c] InteractiveUtils [@stdlib/InteractiveUtils]
[76f85450] LibGit2 [@stdlib/LibGit2]
[56ddb016] Logging [@stdlib/Logging]
[d6f4376e] Markdown [@stdlib/Markdown]
[44cfe95a] Pkg [@stdlib/Pkg]
[de0858da] Printf [@stdlib/Printf]
[3fa0cd96] REPL [@stdlib/REPL]
[9a3f8284] Random [@stdlib/Random]
[ea8e919c] SHA [@stdlib/SHA]
[9e88b42a] Serialization [@stdlib/Serialization]
[6462fe0b] Sockets [@stdlib/Sockets]
[8dfed614] Test [@stdlib/Test]
[cf7118a7] UUIDs [@stdlib/UUIDs]
[4ec0a83e] Unicode [@stdlib/Unicode]
ERROR: LoadError: ArgumentError: Package RandomNumbers does not have Requires in its dependencies:

• If you have RandomNumbers checked out for development and have
  added Requires as a dependency but haven't updated your primary
  environment's manifest file, try Pkg.resolve().
• Otherwise you may need to report an issue with RandomNumbers
  Stacktrace:
  [1] require(::Module, ::Symbol) at ./loading.jl:836
  [2] include at ./boot.jl:326 [inlined]
  [3] include_relative(::Module, ::String) at ./loading.jl:1038
  [4] include(::Module, ::String) at ./sysimg.jl:29
  [5] top-level scope at none:2
  [6] eval at ./boot.jl:328 [inlined]
  [7] eval(::Expr) at ./client.jl:404
  [8] top-level scope at ./none:3
  in expression starting at /Users/jeffreyvarner/.julia/packages/RandomNumbers/aomjt/src/RandomNumbers.jl:32
  ERROR: LoadError: LoadError: Failed to precompile RandomNumbers [e6cf234a-135c-5ec9-84dd-332b85af5143] to /Users/jeffreyvarner/.julia/compiled/v1.1/RandomNumbers/pgCpR.ji.
  Stacktrace:
  [1] error(::String) at ./error.jl:33
  [2] compilecache(::Base.PkgId, ::String) at ./loading.jl:1197
  [3] _require(::Base.PkgId) at ./loading.jl:960
  [4] require(::Base.PkgId) at ./loading.jl:858
  [5] require(::Module, ::Symbol) at ./loading.jl:853
  [6] include at ./boot.jl:326 [inlined]
  [7] include_relative(::Module, ::String) at ./loading.jl:1038
  [8] include at ./sysimg.jl:29 [inlined]
  [9] include(::String) at /Users/jeffreyvarner/.julia/packages/RandomNumbers/aomjt/test/runtests.jl:7
  [10] top-level scope at none:0
  [11] eval(::Module, ::Any) at ./boot.jl:328
  [12] top-level scope at /Users/jeffreyvarner/.julia/packages/RandomNumbers/aomjt/test/runtests.jl:7
  [13] include at ./boot.jl:326 [inlined]
  [14] include_relative(::Module, ::String) at ./loading.jl:1038
  [15] include(::Module, ::String) at ./sysimg.jl:29
  [16] include(::String) at ./client.jl:403
  [17] top-level scope at none:0
  in expression starting at /Users/jeffreyvarner/.julia/packages/RandomNumbers/aomjt/test/common.jl:1
  in expression starting at /Users/jeffreyvarner/.julia/packages/RandomNumbers/aomjt/test/runtests.jl:1
  ERROR: Package RandomNumbers errored during testing

I noticed a commit a few hours ago (on 2/14/19 around 2:30 EST USA) to fix this issue, was fine before the commit. Broken after?

I think I'll be working on this package with Julia 0.5, but as long as it can run well at 0.4 I won't change the REQUIRE file.

Any advice?

https://nullprogram.com/blog/2017/09/21/
"xoroshiro128+ fails PractRand very badly. [..] Since this article was published, its authors have supplanted it with xoshiro256**. It has essentially the same performance, but better statistical properties. xoshiro256** is now my preferred PRNG"

See also my JuliaLang/julia#27614 (comment) for others to possibly implement (mostly Google's Randen based on AES, the one I list here would however be my first choice) and updated top comment there.

RandomNumbers.jl fails to build on julia 0.6

`julia> Pkg.build("RandomNumbers")
INFO: Building RandomNumbers
====================================================================[ ERROR: RandomNumbers ]=====================================================================

LoadError: LoadError: this intrinsic must be compiled to be called
while loading /Users/snirgaz/.julia/v0.6/RandomNumbers/deps/Random123/build.jl, in expression starting on line 38
while loading /Users/snirgaz/.julia/v0.6/RandomNumbers/deps/build.jl, in expression starting on line 1

=================================================================================================================================================================

========================================================================[ BUILD ERRORS ]=========================================================================

WARNING: RandomNumbers had build errors.

• packages with build errors remain installed in /Users/snirgaz/.julia/v0.6
• build the package(s) and all dependencies with Pkg.build("RandomNumbers")
• build a single package by running its deps/build.jl script

=================================================================================================================================================================

`

I am getting this error while trying to build RandomNumbers in Julia 1.0.

┌ Error: Error building `RandomNumbers`: 
│ ERROR: LoadError: LoadError: UndefVarError: is_bsd not defined

Is there a plan to register this package? I would like to be able to start using it to some extent, but it cannot be a dependency until it's registered.

One important aspect of this project is that the RNGs will be compared together. I want to discuss about the methods of benchmark here.

I have written a page on testing PRNGs for high-quality randomness. My goal is to use the results of statistical randomness tests to inform decisions on which PRNGs are high-quality, more specifically, by testing whether the PRNG produces a sequence of numbers that behave like independent uniform random numbers. If a PRNG does well in statistical testing (and meets my other criteria of a high-quality PRNG), only then should we discuss its other aspects (such as performance and ease of implementation).

The following are some points on PRNG testing:

• There are many kinds of PRNGs, including traditional, counter-based, and splittable PRNGs, as well as PRNGs that combine two or more other PRNGs, and there are different approaches to testing all these kinds.

• Many PRNGs have a cycle length so large that no (empirical) statistical test can evaluate all the numbers in the cycle in a reasonable time. At most, the test can only look at a small fraction of the PRNG's behavior (e.g., PractRand has a limit of 32 TiB, or 2^46 bytes). The most demanding simulations may require still fewer random numbers, e.g., 10^11 numbers at most. I currently consider a PRNG statistically good if it doesn't fail PractRand at 1 TiB (2^40 bytes) (where I don't consider "unusual" or "suspicious" results, as opposed to "FAILs", to be failures). Is this threshold too low or too high?

• Different PRNGs have different ways to produce independent random number sequences ("streams"), such as by incrementing a seed, or discarding a huge number of outputs. A given stream strategy can be tested by interleaving the outputs of those streams. I am currently testing two or four PRNGs with consecutive seeds, as well as two or four PRNGs that are "jumped ahead" in the manner supported by the PRNG.

• A relatively low-quality PRNG can be combined with another number sequence to produce a high-quality PRNG, enough to show no PractRand failure at 1 TiB (see "Combined PRNGs" for details). A notable example is L'Ecuyer's combined multiple recursive generators.

• Counter-based PRNGs and splittable PRNGs (such as the PRNG used in JAX) can be based on hash functions and other mixing functions, so it's useful to test these functions to see if they are viable for use in these kinds of PRNGs. (I give some of these constructions in detail.) On the subject of mixing functions, a testing procedure by P. Evensen may be of interest.

It seems there are several extremely poor seeds for Xorshift128Plus which result in repeating sequences. For example, seed 512:

julia> rand(Xorshifts.Xorshift128Plus(512), 30)
30-element Vector{Float64}:
 0.12500214576741975
 0.1250021457673638
 0.015626089647719832
 0.15625003352785027
 0.2656260896476079
 0.12500214576747748
 0.015627162531226313
 0.15625217929502666
 0.14062608964765588
 1.1368683772161603e-13
 1.0728837338547237e-6
 2.1457675316582936e-6
 0.12500214576741975
 0.1250021457673638
 0.015626089647719832
 0.15625003352785027
 0.2656260896476079
 0.12500214576747748
 0.015627162531226313
 0.15625217929502666
 0.14062608964765588
 1.1368683772161603e-13
 1.0728837338547237e-6
 2.1457675316582936e-6
 0.12500214576741975
 0.1250021457673638
 0.015626089647719832
 0.15625003352785027
 0.2656260896476079
 0.12500214576747748

Examining the first million seeds, we can easily find several such examples which seem to fail in a similar way:

julia> for i = 1:1_000_000
           xs = rand(Xorshifts.Xorshift128Plus(i), 30)
           n = length(unique(xs))
           if n != length(xs)
               @info "Repeated numbers for $i (length $n)"
           end
       end
[ Info: Repeated numbers for 512 (length 12)
[ Info: Repeated numbers for 1024 (length 12)
[ Info: Repeated numbers for 1536 (length 12)
[ Info: Repeated numbers for 2048 (length 12)
[ Info: Repeated numbers for 2560 (length 12)
[ Info: Repeated numbers for 3072 (length 12)
[ Info: Repeated numbers for 3584 (length 12)

Is it possible to issue a new release (to include MT fix and speed benchmarks)?

I took the fastest version of fasta at the Debian Benchmark Game, I think it was that one (I can send you my full version):
https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/fasta-julia-4.html

I thought I could make it faster with your RNG (yes, not allowed by the rules, just testing).

I changed and added below but it got slower!

time ~/julia-1.4.0-DEV-8ebe5643ca/bin/julia -O3 fasta_xor.jl 250000 >/dev/null

I thought Xoroshiro (maybe not this variant?) was one of the fastest RNG. I'll check others, but thought I should let you know as it's even slower than Julia's default too, not just compared this 18-bit non-standard RNG in the program I trying out better ones for. And if I add the modulus(%) to other RNGs (needed for the same range, as it's used in the original code), it just gets even slower, that limits numbers to 18-bit (at first I thought it to 16-bit so I was first trying such a replacement).

const IM = UInt32(139968)
const IA = UInt32(3877)
const IC = UInt32(29573)
const last_rnd = Ref(UInt32(42))
#gen_random() = (last_rnd[] = (last_rnd[] * IA + IC) % IM)

using RandomNumbers.Xorshifts
r2 = Xoroshiro128Plus(0x1234567890abcdef)  # with a certain seed. Note that the seed must be non-zero.
# gen_random() = UInt32(rand(r2, UInt16)) #  % IM) # also slower with UInt16 there, and I think I reported for that, not the even slower:
gen_random() = UInt32(rand(r2, UInt32)) #  % IM)

Version 1.1.1 of Julia doesn't change much or -O2 vs -O3.

Possibly this is because of my very old laptop (would be nice if you checked):

julia> versioninfo()
Julia Version 1.4.0-DEV.17
Commit 8ebe5643ca (2019-08-20 08:32 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: Intel(R) Core(TM)2 Duo CPU P8400 @ 2.26GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-6.0.1 (ORCJIT, penryn)
Environment:
JULIA_NUM_THREADS = 1

If this isn't just a problem with my laptop, it could be inlining related? And if the RNG in that program is simply faster, then maybe (even with it probably then a flawed RNG, broken by RNGCrush?) could be an option for your library with a clear WARNING to users.

Note, using RandomNumbers.Xorshifts is noticeably slow and you might think I'm not taking it into account, but running for a minute+ (with larger number on the command-line when calling the program) doesn't help.

These two lines of tempering transform should use left shift:
https://github.com/sunoru/RandomNumbers.jl/blob/v1.2.0/src/MersenneTwisters/bases.jl#L77-L78.

As a verification, when seeded with 5489, MT19937 generates a UInt32 sequence of 0xD091BB5C, 0x22AE9EF6, 0xE7E1FAEE, ....

I see that gen_seed is a common function for all RNGs in RandomNumbers. It would be useful to able to call gen_seed(rng) and have it return the appropriately sized tuple for that rng. This could make code more generic.
I can't see whether this is currently possible. My workaround has been:

Base.length(::Type{NTuple{N,T}}) where {N,T} = N
seed_length(rng) = length(seed_type(rng))
gen_seed(rng::RandomNumbers.AbstractRNG) = gen_seed(output_type(rng), seed_length(rng))

Not sure if Base.length is the right function to extend here, or if there already exists a function for this.

When building, Pkg.build("RandomNumbers"), on Windows without gcc, get error due to this:

check_compiler() = success(`gcc --version`)

Could replace by the following which fails more gracefully:

function check_compiler()
    try
        run(`gcc --version`)
        return true
    end
    return false
end

However, I noticed that build() on Windows tries run(`mingw32-make`),
and if that fails then downloads a binary dll.
In this case should the check_compiler() be removed?

if have_aesni() #&& check_compiler()
    build()
else
    warn("AES-NI will not be compiled.")
end

I commented out the compiler check, and the build seemed to work.
Also Pkg.test("RandomNumbers") reported that all tests passed.

If this is intended then definition and calls to check_compiler() should be removed?

I tested one common PCG declared and initialized as:

using RNG.PCG
r = PermutedCongruentialGenerator(PCGStateSetseq, PCG_XSH_RR, UInt64)
srand(r, (42 % UInt64, 54 % UInt64))

It passed all the tests of the TestU01 big crush battery as expected, but the long CPU time frightened me. It took over 21 hours on my laptop.

========= Summary results of BigCrush =========

 Version:          TestU01 1.2.3
 Generator:
 Number of statistics:  160
 Total CPU time:   21:32:37.44

 All tests were passed

For comparison, the test results of Base.Random.rand are:

========= Summary results of BigCrush =========

 Version:          TestU01 1.2.3
 Generator:
 Number of statistics:  160
 Total CPU time:   03:49:04.43
 The following tests gave p-values outside [0.001, 0.9990]:
 (eps  means a value < 1.0e-300):
 (eps1 means a value < 1.0e-15):

       Test                          p-value
 ----------------------------------------------
 80  LinearComp, r = 0              1 - eps1
 81  LinearComp, r = 29             1 - eps1
 ----------------------------------------------
 All other tests were passed

Obviously it's unacceptable with such slow speed. I've tried to use @code_warntype to find where could be the bottleneck, and it shows some strange positions of codes are analyzed to return Any, which was thought to be optimized. For example,

julia> @code_warntype rand(r, UInt32)
Variables:
  #self#::Base.Random.#rand
  pcg::RNG.PCG.PermutedCongruentialGenerator{RNG.PCG.PCGStateSetseq{StateUIntType<:Union{UInt128,UInt16,UInt32,UInt64,UInt8},MethodType<:Union{Val{:RXS_M_XS},Val{:XSH_RR},Val{:XSH_RS},Val{:XSL_RR_RR},Val{:XSL_RR}}},Val{:XSH_RR},UInt64,UInt32}
  #unused#::Type{UInt32}

Body:
  begin 
      $(Expr(:meta, :inline)) # /home/sunoru/.julia/v0.5/RNG/src/./PCG/main.jl, line 40:
      return (top(typeassert))((RNG.PCG.pcg_random)((top(getfield))(pcg::RNG.PCG.PermutedCongruentialGenerator{RNG.PCG.PCGStateSetseq{StateUIntType<:Union{UInt128,UInt16,UInt32,UInt64,UInt8},MethodType<:Union{Val{:RXS_M_XS},Val{:XSH_RR},Val{:XSH_RS},Val{:XSL_RR_RR},Val{:XSL_RR}}},Val{:XSH_RR},UInt64,UInt32},:state)::RNG.PCG.PCGStateSetseq{StateUIntType<:Union{UInt128,UInt16,UInt32,UInt64,UInt8},MethodType<:Union{Val{:RXS_M_XS},Val{:XSH_RR},Val{:XSH_RS},Val{:XSL_RR_RR},Val{:XSL_RR}}})::Any,$(Expr(:static_parameter, 4)))::UInt32
  end::UInt32

How can I fix this problem?

Hi, I thought the Xoroshiro128Plus was super fast. The docs state:

a significant improvement in speed and in statistical quality. Therefore, Xoroshiro128Plus is the current best suggestion for replacing other low-quality generators

I have this package in my juliarc.jl file because of this. But

using BenchmarkTools, RandomNumbers
r = Xoroshiro128Plus()
myrand = rand(r)

@benchmark rand()
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     6.291 ns (0.00% GC)
  median time:      6.658 ns (0.00% GC)
  mean time:        7.337 ns (0.00% GC)
  maximum time:     80.114 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000

@benchmark myrand()
BenchmarkTools.Trial: 
  memory estimate:  16 bytes
  allocs estimate:  1
  --------------
  minimum time:     27.372 ns (0.00% GC)
  median time:      28.381 ns (0.00% GC)
  mean time:        36.034 ns (2.64% GC)
  maximum time:     2.122 μs (91.72% GC)
  --------------
  samples:          10000
  evals/sample:     995

It appears in fact to be several factors slower?

A PkgEval run for a Julia pull request which changes the generated numbers for rand(a:b) indicates that the tests of this package might fail in Julia 1.5 (and on Julia current master branch).

Apologies if this is a false positive. Cf. https://github.com/JuliaCI/NanosoldierReports/blob/ab6676206b210325500b4f4619fa711f2d7429d2/pkgeval/by_hash/52c2272_vs_47c55db/logs/RandomNumbers/1.5.0-DEV-87d2a04de3.log

Is anyone working on that? Thank you!

Previously I realized how bad the current implementation of using integers with small entropy to initialize a big-state RNG is. If we directly use small integers like 1, 2, 3 to be the initial states of xor-based RNGs (such as Xorshift), the starting sequences will probably poorly random. This is caused by the zeros in the state.
E.g.: https://sunoru.github.io/RNG.jl/latest/man/mersenne-twisters/#Examples-1

MT19937 uses an initialization function to avoid this situation. It performs a iterative transformation on one UInt32 seed to generate a 623-UInt32s state. (And I still don't understand how people come up with the constants 1812433253 and 6364136223846793005)

C++'s seed_seq may be a valuable reference, although this blog shows its flaws.

I have looked at how dSFMT and libstdc++ deal with it:
https://github.com/MersenneTwister-Lab/dSFMT/blob/c6bf8a8dab3710b7abd86c4b68d0e6b4aa5e6db1/dSFMT.c#L554
https://github.com/atgreen/gcc/blob/76cc869aa9710764c7732d9bca0740fcf6865a27/libstdc%2B%2B-v3/include/bits/random.tcc#L3398
and they seem to be using similar methods..

@simonbyrne

Julia Version 1.0.0
Commit 5d4eaca0c9 (2018-08-08 20:58 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-4500U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, haswell)

I get the following:

┌ Error: Error building `RandomNumbers`:
│ ERROR: LoadError: LoadError: UndefVarError: info not defined
│ Stacktrace:
│  [1] build() at C:\Users\nicol\.julia\packages\RandomNumbers\cpJea\deps\Random123\build.jl:13
│  [2] top-level scope at C:\Users\nicol\.julia\packages\RandomNumbers\cpJea\deps\Random123\build.jl:40
│  [3] include at .\boot.jl:317 [inlined]
│  [4] include_relative(::Module, ::String) at .\loading.jl:1038
│  [5] include(::Module, ::String) at .\sysimg.jl:29
│  [6] include(::String) at .\client.jl:388
│  [7] top-level scope at none:0
│  [8] include at .\boot.jl:317 [inlined]
│  [9] include_relative(::Module, ::String) at .\loading.jl:1038
│  [10] include(::Module, ::String) at .\sysimg.jl:29
│  [11] include(::String) at .\client.jl:388
│  [12] top-level scope at none:0
│ in expression starting at C:\Users\nicol\.julia\packages\RandomNumbers\cpJea\deps\Random123\build.jl:39
│ in expression starting at C:\Users\nicol\.julia\packages\RandomNumbers\cpJea\deps\build.jl:1
└ @ Pkg.Operations C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.0\Pkg\src\Operations.jl:1068

I think that info should be @info in https://github.com/sunoru/RandomNumbers.jl/blob/838e69e8199941b2331dd98836df669b75b49264/deps/Random123/build.jl#L13

This library seems to be pure Julia except for Random123. I would like to see it end up as pure Julia because that makes it much easier to install and transfer in some ways. Then on the page for them you can just say that using Random123 needs to be done. I don't think that would hurt because they don't seem to perform that well, so adding the binary dependency for what is not likely to be used is likely not a good decision now that this is past research and in use.

On one of my PCs (which runs Windows 7), Pkg.test() fails.

Files in the expected sub-directory contain CR-LF (\r\n) pairs, whereas files created in the actual sub-directory contain only LF characters.

On another PC (which runs Windows 10), files in both sub-directories contain LF character only, and Pkg.test() succeeds.

Perhaps the GitHub download of expected files is different for these PCs.
I'm not sure what is different between these two PCs, other than the Windows version.

I have a Julia fix, which simply replaces CR-LF with LF before testing, and could submit a PR.
However, I was wondering if someone else can reproduce this and confirm it is a general issue, rather than just something peculiar to my machine.

I see this library has been updated on master but it hasn't been tagged.

Just a short question because I have not found any hint for this in the docs. Why is the performance of the Base Mersenne Twister much faster than the performance of this one implemented here.

MWE:

using RandomNumbers
rng = Xorshifts.Xoroshiro128Plus(0)
randn(rng,Float64) # 0.0

Right link is: http://sunoru.github.io/RandomNumbers.jl/latest/ I think
But currently it's linking to: http://sunoru.github.io/RandomNumbers.jl/latest/man/bases/

Great work on this package - thank you!

This is a question rather than an issue. I hope that's OK to ask here.

I would like to use a RNG for parallel simulations, and so want statistically independent RNGs for each worker process.

For Base.Random.MersenneTwister, I use something like this:

rng = Base.Random.MersenneTwister(0)
rngseed = Base.Random.make_seed()
srand(rng, rngseed)

rngs = randjump(rng, numtrials)

f(trial, rng) = ...
pmap(f, 1:numtrials, rngs)

I would like to use a RNG from the Random123 family.
Is this a good choice for parallel sims?

My code would be something like:

rng = RandomNumbers.Random123.Threefry4x()
rngseed = RandomNumbers.gen_seed(rng)
srand(rng, rngseed)

rngs = Vector{typeof(rng)}(numtrials)
for trial = 1:numtrials
    rngs[trial] = copy(rng)
    RandomNumbers.Random123.set_counter!(rngs[trial], trial)
end

f(trial, rng) = ...
pmap(f, 1:numtrials, rngs)

Is using set_counter! this way the right approach to give independent RNGs?
Will it work for all RNGs in the Random123 family?
Why would I choose 2x or 4x version?
Which particular RNG would you recommend?

It seems almost all the implementations choose the (a,b,c) parameters arbitrarily. Some require input by the user. I'm wondering what our package should provide...

ref: http://www.jstatsoft.org/v08/i14/paper

MethodError is raised when trying to use some RNGs from this package.

Steps to reproduce (for example):
using RNG.Random123; randn(Philox4x())

Haven't yet investigated whether this is a package issue or upstream, but want to file an issue somewhere so I don't forget.

On Windows 10 Julia version 1.0.0 I install RandomNumbers. There is a failure.
When running the build tests:

WARNING: could not import Random.srand into T1
[ Info: Testing Generic functions
WARNING: could not import Random.srand into T2
[ Info: Testing WrappedRNG
WARNING: could not import Random.srand into T3
[ Info: Testing PCG
WARNING: could not import Random.srand into T4
[ Info: Testing MersenneTwisters
WARNING: could not import Random.srand into T5
[ Info: Testing Xorshifts
ERROR: LoadError: LoadError: UndefVarError: srand not defined
Stacktrace:
[1] top-level scope at C:\Users\hearn.julia\packages\RandomNumbers\5tfiq\test\Xorshifts\runtests.jl:39 [inlined]
[2] top-level scope at .\none:0
[3] include at .\boot.jl:317 [inlined]
[4] include_relative(::Module, ::String) at .\loading.jl:1038
[5] include at .\sysimg.jl:29 [inlined]
[6] include(::String) at C:\Users\hearn.julia\packages\RandomNumbers\5tfiq\test\runtests.jl:8
[7] top-level scope at none:0
[8] eval(::Module, ::Any) at .\boot.jl:319
[9] top-level scope at C:\Users\hearn.julia\packages\RandomNumbers\5tfiq\test\runtests.jl:8
[10] include at .\boot.jl:317 [inlined]
[11] include_relative(::Module, ::String) at .\loading.jl:1038
[12] include(::Module, ::String) at .\sysimg.jl:29
[13] include(::String) at .\client.jl:388
[14] top-level scope at none:0
in expression starting at C:\Users\hearn.julia\packages\RandomNumbers\5tfiq\test\Xorshifts\runtests.jl:19
in expression starting at C:\Users\hearn.julia\packages\RandomNumbers\5tfiq\test\runtests.jl:1
ERROR: Package RandomNumbers errored during testing

This is more of a question than an issue. Apologies if this is not the right place!

I'm porting some code to using RandomNumbers because I need parallel random number generation. I'm using an array of Random123.Threefry4x as an rng for each thread.

My current code does something like rand(Multinomial(1, multinomial_probabilities))). Is there a way I can pass in the specific thread rng to be used for the Multinomial draw rather than the default global rng?

Thanks,
Rohit

As discussed before, I'd like to implement various RNGs in this package. And temporarily they are listed in the following:

Four from Random123:

• Threefry
• Philox
• AESNI
• ARS

The xorshift families:

• xoroshiro128+
• xorshift128+
• xorshift1024*

Six from SPRNG

• Combined Multiple Recursive Generator
• 48 Bit Linear Congruential Generator with Prime Addend
• 64 Bit Linear Congruential Generator with Prime Addend
• Modified Lagged Fibonacci Generator
• Multiplicative Lagged Fibonacci Generator
• Prime Modulus Linear Congruential Generator

In fact if we put the parallel stuff aside, the work seems pretty easy to be done. I have roughly read all the docs of the RNGs mentioned in #1 and have some ideas.

TODOs is generally as the followings (in order):

• Complete the design of test suites and main framework of RNGs. #4
• Complete a method of benchmark. It will also used on the current Base.rand. #2
• Implement the various RNGs in pure Julia. Make sure that they can pass the tests.
• Dig into the codes for better performance, maybe write some C or llvm codes and even do some changes to Julia itself.
• Select the RNG of the best performance. Try to merge it to Base.Random.
• Take parallel into account.
• Documentation and result reports etc.

In my opinion it is suitable to work on the first three by the midterm evaluation.

Issues have been created in various aspects, and this is one for general discussion. I'm willing to see any advice.

Is it possible to store (to a string/file etc.) and later load the RNG state in order to resume the random stream?

Tests of random number generators are very important.

TestU01 is a great library to do this, and there is a ready-made binding. Its user guide has showed nice reliability.

However the standard crush test battery of TestU01 are too slow (which would exceed the limitation of travis CI), so I think RNGTest.smallcrushJulia may be a good choice in runtests.jl and manually run RNGTest.crushTestU01 on my own computer.

Moreover, the tests in SPRNG seems also helpful. It has physical model tests, and some useful tests for parallel algorithms. It may also be considered.

From http://xoroshiro.di.unimi.it

All generators, being based on linear recurrences, provide jump functions that make it possible to simulate any number of calls to the next-state function in constant time, once a suitable jump polynomial has been computed. We provide ready-made jump functions for a number of calls equal to the square root of the period, to make it easy generating non-overlapping sequences for parallel computations.

Would it be possible to implement jump functions in this package?
This seems like a nice addition for parallel usage.

I am using Julia version 5. If I try the following:

using RandomNumbers.Xorshifts

r = Xoroshiro128Star()
arr = [1,2,3,4]
shuffle!(r, arr)

It will explode with a MethodError.

I was wondering whether there are plans to add to this package of family of very cheap PRNG. I, personally, would like to compare Xoroshiro to something that I just came across:

// Cheap PRNG from "Numerical Recipes in C", page 284.
static inline uint32_t cheap_generator () {
	static unsigned long idum;
	idum = 1664525L*idum + 1013904223L;

	return idum;
}

Sure, that's so simple that one could just implement that directly, but it might be nice to have it as part of this package. Some applications (chaotic PDEs for me) probably don't need any high quality PRNG, therefore testing it against the (presumably) cheapest possible PRNG would be great.

As it's written in the documentation, the standard approach of for rand(Float32)/rand(Float64)

An UInt64 will firstly be converted to a Float64 that is perfectly uniformly distributed in [1.0, 2.0), and then minus one. This is a very fast approach, but not completely ideal, since the one bit is wasted. The current default RNG in Base.Random library does the same thing, so it also causes some tricky problems

is probably good enough in most cases, but not perfect. For Float32 it samples only from 2^23 = 8,388,608 floats in [+0,1) (as that’s the number of floats in [1,2)) whereas there are 1,065,353,216 (=2^30-2^23) in that range in total. So rand(Float32) samples from about 1% possible floats, and rand(Float64) from ~0.1% (2^52/(2^62-2^52)). I propose a function that solves that problem by drawing the exponent bits at correct chances first and then combine them with significant bits, similar to the following

using RandomNumbers.Xorshifts
const Xoroshiro128P = Xoroshiro128Plus()

"""Random number generator for floats in [0,1) that samples from *all* floats.""" 
function randfloat(::Type{Float32})
    # create exponent bits in 0000_0000 to 0111_1111
    # at following chances
    # e=01111110 at 50.0% for [0.5,1.0)
    # e=01111101 at 25.0% for [0.25,0.5)
    # e=01111100 at 12.5% for [0.125,0.25)
    # ... etc
    ui = rand(Xoroshiro128P,UInt128) >> 1 | 0x1 # first bit always 0, last always 1
    lz = leading_zeros(ui) % UInt8              # 2 leading zeros is 50%, 3 is 25%, etc.  
    e = ((0x7f - lz) % UInt32) << 23            # convert lz to exponent bits of float32

    # significant bits
    f = rand(Xoroshiro128P,UInt32) >> 9 
    
    # combine exponent and significand (sign always 0)
    return reinterpret(Float32,e | f)
end

I don't know whether something like this exists already in Julia somewhere, I don't even know whether this is an established algorithm. I thought it would be nice to have this function and different to implemented here, it could take the actual RNG as an argument. I've just used Xoroshiro128+ as it's fast. I'm happy to implement a Float64 version too, this requires technically a 1024bit UInt, but in practice one could first draw a UInt64 and only if 0x0 is drawn (does that actually happen?) one draws another UInt64 and count the leading zeros again. The above implementation seems to be reasonably fast, about 3x slower than Xoroshiro128+ with the "[1,2)-1"-technique but still faster than Julia's Base rand(Float32):

julia> @btime randfloat($Float32)   # this approach
  4.603 ns (0 allocations: 0 bytes)
julia> @btime rand($Float32)          # Julia's Base MT19337
  6.675 ns (0 allocations: 0 bytes)
julia> @btime rand($Xoroshiro128P,$Float32)    # Xoroshiro128+ with [1,2)-1 technique
  1.544 ns (0 allocations: 0 bytes)

Again, there might be ways to speed this up further by first drawing a UInt64 and only if that's 0x0 draw another one.

Dear @sunoru

Very nice job on RandomNumbers. I'm trying to run the benchmarks locally, but I get an error:

julia> include("speeds.jl")
ERROR: LoadError: LoadError: UndefVarError: RNG not defined
Stacktrace:
 [1] include_from_node1(::String) at ./loading.jl:576
 [2] include(::String) at ./sysimg.jl:14
 [3] include_from_node1(::String) at ./loading.jl:576
 [4] include(::String) at ./sysimg.jl:14
while loading common.jl, in expression starting on line 4
while loading speeds.jl, in expression starting on line 7

I think you should add to the basic section of the documentation how to use this library for generating numbers from some basic distributions. A mention of how to use it with randn is first, but then also how to (if it's possible) use it with Distributions.jl for random numbers from a larger class of distributions. I believe everything to do these should already be implemented (they should have an argument spot for the RNG?), so it's likely just a documentation thing to make it more accessible. If not, appropriate changes around Julia should be made to allow using these faster RNGs in place of the Base RNG for calculating random numbers from distributions.

AES encryption is an important part of the RNGs for secure use. There are some methods to use AES:

1. Implement the encryption algorithm in pure Julia. There are also existing packages such as AES.jl (but this is not registered), and Nettle.jl (though this also provides many other functions).
2. Use the mature OpenSSL. This may also have two different ways: to wrap it by ourselves or use the existing ones (for example, Crypto.jl, but it is not completed so doesn't guarantee the security)...

It would be better if we can directly use AES-NI on the most kinds of CPU, so I'm also wondering how will llvm do with this? Anyway, OpenSSL should have been already optimized about such things.

You had a comment about moving circular shifts to assembly. One alternative is to use llvmcall:

function llvm_rotr(x::UInt64, y::UInt64)
    Base.llvmcall("""
                %3 = lshr i64 %0, %1
                %4 = sub i64 64, %1
                %5 = shl i64 %0, %4
                %6 = or i64 %3, %5
                ret i64 %6
              """,UInt64, Tuple{UInt64,UInt64},x,y)
end

then using 0.4, or 0.5 with -O3 (which is how packages are precompiled), gives the appropriate rorq instruction.

The "obvious" approach of converting to a float, then multiplying by a scale factor is slow, and has the potential to give an answer of 1.0.

The easiest option is:

import Base: significand_mask, exponent_one
f1(u::UInt64) = reinterpret(Float64, exponent_one(Float64) | significand_mask(Float64) & u) - 1.0

unfortunately this has the downside that the last bit will always be zero, so you only get 52 bits of randomness per float: see JuliaLang/julia#16344.

A slightly more advanced option (based on this proposal) is:

import Base: significand_mask, significand_bits, exponent_half
function f2(u::UInt64)
    f = reinterpret(Float64, exponent_half(Float64) | 0x001f_ffff_ffff_ffff & u)
    if (u >> significand_bits(Float64)) &1 == 1
        f-= 0.5
    end
    f
end

This gets us 53 bits, but introduces a branch. There might be some clever stuff we can do here though.

I think "re-using" the name AbstractRNG is not ideal:
abstract type AbstractRNG{T<:Number} <: Base.Random.AbstractRNG end

AbstractRNG is defined and exported by Base and RandomNumbers

using RandomNumbers
AbstractRNG

results in following:

WARNING: both RandomNumbers and Base export "AbstractRNG"; uses of it in module Main must be qualified
ERROR: UndefVarError: AbstractRNG not defined

Thus code that uses AbstractRNG without RandomNumbers, will no longer work when using RandomNumbers, without explicitly qualifying. i.e. that code will need modifying.

This can be problematic when one doesn't have access to the modules containing references to AbstractRNG.

Maybe this is some form of type piracy?

It would be nice to have a few chaotic PRNGs, which are highly praised by some practitioners.

E.g --

1. gjrand -- http://gjrand.sourceforge.net/ and
2. jsf -- http://burtleburtle.net/bob/rand/smallprng.html
3. sfc -- https://gist.github.com/imneme/f1f7821f07cf76504a97f6537c818083

Thanks!

From what I've read, Xoroshiro is often suggested as the fastest RNG with decent random properties. Using the implementation here, I do see it beating Base's MT for small samples, but getting relatively slower as the sample size passes ~100:

julia> r0 = Random.default_rng();

julia> r1 = RandomNumbers.Xorshifts.Xoroshiro128Plus(1);

julia> @btime rand($r0);
  2.458 ns (0 allocations: 0 bytes)

julia> @btime rand($r1);
  1.925 ns (0 allocations: 0 bytes)

julia> @btime rand($r0, 8);
  41.771 ns (1 allocation: 144 bytes)

julia> @btime rand($r1, 8);
  34.454 ns (1 allocation: 144 bytes)

julia> @btime rand($r0, 128);
  171.495 ns (1 allocation: 1.14 KiB)

julia> @btime rand($r1, 128);
  180.105 ns (1 allocation: 1.14 KiB)

julia> @btime rand($r0, 10^7);
  6.606 ms (2 allocations: 76.29 MiB)

julia> @btime rand($r1, 10^7);
  25.872 ms (2 allocations: 76.29 MiB)

Is this to be expected? Perhaps this information could be added to the docs.

julia> versioninfo()
Julia Version 1.5.0-DEV.803
Commit 8ab87d2205 (2020-05-03 08:07 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
  JULIA_NUM_THREADS = 8

julia> versioninfo()
Julia Version 0.7.0-DEV.1758
Commit 8757abf (2017-09-11 20:45 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, ivybridge)
Environment:

julia> using RandomNumbers.Random123

julia> key = 0x88966ab1262b56a53efae875dab0d616
0x88966ab1262b56a53efae875dab0d616

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0x5318a25ddf82379b96f0f594b2a7f096, RandomNumbers.Random123.AESNIKey(0x00007fa41bc721250000000000000011, 0x52a43d4352a442e7496363c2496363c2, 0x1a00368148a40bc21a00492553632ae7, 0x17653d800d650b0145c100c35fc149e6, 0xcdf0328bda950f0bd7f0040a923104c9, 0x6f19b170a2e983fb787c8cf0af8c88fa, 0x4ba8e26924b153198658d0e2fe245c12, 0xeed6ff58a57e1d3181cf4e2807979eca, 0xa7d8c41d490e3b45ec7026746dbf685c, 0xcb45d0776c9d146a25932f2fc9e3095b, 0xbeb78c2f75f25c58196f48323cfc671d), 0x00000000000000000000000000000000)

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0xe067f92982571746f3ed5fb6dd484ea9, RandomNumbers.Random123.AESNIKey(0x00007fa406ec001000007fa4088551e0, 0x470a3223470a4d8741e64d9741e63233, 0x26a0672161aa550226a0188567465512, 0xfb1b9f35ddbbf814bc11ad169ab1b593, 0x960fd0776d144f42b0afb7560cbe1a40, 0xb29a4443249594344981db76f92e6c20, 0x3c97df1a8e0d9b59aa980f6de319d41b, 0x59f017eb6567c8f1eb6a53a841f25cc5, 0x7fc45c0726344bec4353831da839d0b5, 0x77485812088c04152eb84ff96debcce4, 0xf5628d46822ad5548aa6d141a41e9eb8), 0x00000000000000000000000000000000)

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0xe067f92982571746f3ed5fb6dd484ea9, RandomNumbers.Random123.AESNIKey(0x00007fa406ec001000007fa4088551e0, 0x470a3223470a4d8741e64d9741e63233, 0x26a0672161aa550226a0188567465512, 0xfb1b9f35ddbbf814bc11ad169ab1b593, 0x960fd0776d144f42b0afb7560cbe1a40, 0xb29a4443249594344981db76f92e6c20, 0x3c97df1a8e0d9b59aa980f6de319d41b, 0x59f017eb6567c8f1eb6a53a841f25cc5, 0x7fc45c0726344bec4353831da839d0b5, 0x77485812088c04152eb84ff96debcce4, 0xf5628d46822ad5548aa6d141a41e9eb8), 0x00000000000000000000000000000000)

It gets a different result when do ccall for initkey the first time, but the subsequent results are wrong as well. It worked fine before 0.7 with the correct results:

julia> versioninfo()
Julia Version 0.6.0-rc1.0
Commit 6bdb395 (2017-05-07 00:00 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, ivybridge)

julia> using RandomNumbers.Random123
INFO: Recompiling stale cache file /home/sunoru/.julia/lib/v0.6/RandomNumbers.ji for module RandomNumbers.
key = 
julia> key = 0x88966ab1262b56a53efae875dab0d616
0x88966ab1262b56a53efae875dab0d616

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0x7b69c3e1f0c2a2bb3fb27e95b37ae3d3, RandomNumbers.Random123.AESNIKey(0x88966ab1262b56a53efae875dab0d616, 0x823392740aa5f8c52c8eae6012744615, 0x247f4189a64cd3fdace92b3880678558, 0x098bee932df4af1a8bb87ce7275157df, 0x549757915d1cb90270e81618fb506aff, 0x03131a3f57844dae0a98f4ac7a70e2b4, 0x51043c0b5217263405936b9a0f0b9f36, 0x225a1c38735e20332149060724da6d9d, 0x5304e98d715ef5b50200d5862349d381, 0x5efee83a0dfa01b77ca4f4027ea42184, 0xd15c87a68fa26f9c82586e2bfefc9a29), 0x00000000000000000000000000000000)

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0x7b69c3e1f0c2a2bb3fb27e95b37ae3d3, RandomNumbers.Random123.AESNIKey(0x88966ab1262b56a53efae875dab0d616, 0x823392740aa5f8c52c8eae6012744615, 0x247f4189a64cd3fdace92b3880678558, 0x098bee932df4af1a8bb87ce7275157df, 0x549757915d1cb90270e81618fb506aff, 0x03131a3f57844dae0a98f4ac7a70e2b4, 0x51043c0b5217263405936b9a0f0b9f36, 0x225a1c38735e20332149060724da6d9d, 0x5304e98d715ef5b50200d5862349d381, 0x5efee83a0dfa01b77ca4f4027ea42184, 0xd15c87a68fa26f9c82586e2bfefc9a29), 0x00000000000000000000000000000000)

julia> AESNI1x(key)
RandomNumbers.Random123.AESNI1x(0x7b69c3e1f0c2a2bb3fb27e95b37ae3d3, RandomNumbers.Random123.AESNIKey(0x88966ab1262b56a53efae875dab0d616, 0x823392740aa5f8c52c8eae6012744615, 0x247f4189a64cd3fdace92b3880678558, 0x098bee932df4af1a8bb87ce7275157df, 0x549757915d1cb90270e81618fb506aff, 0x03131a3f57844dae0a98f4ac7a70e2b4, 0x51043c0b5217263405936b9a0f0b9f36, 0x225a1c38735e20332149060724da6d9d, 0x5304e98d715ef5b50200d5862349d381, 0x5efee83a0dfa01b77ca4f4027ea42184, 0xd15c87a68fa26f9c82586e2bfefc9a29), 0x00000000000000000000000000000000)