Merge pull request #14772 from rfourquet/rf/randperm/rand_lt

speed-up randperm, randcycle and shuffle

base/random.jl

@@ -104,7 +104,6 @@ end
 
 @inline rand_ui52_raw_inbounds(r::MersenneTwister) = reinterpret(UInt64, rand_inbounds(r, Close1Open2))
 @inline rand_ui52_raw(r::MersenneTwister) = (reserve_1(r); rand_ui52_raw_inbounds(r))
-@inline rand_ui52(r::MersenneTwister) = rand_ui52_raw(r) & 0x000fffffffffffff
 @inline rand_ui2x52_raw(r::MersenneTwister) = rand_ui52_raw(r) % UInt128 << 64 | rand_ui52_raw(r)
 
 function srand(r::MersenneTwister, seed::Vector{UInt32})
@@ -241,7 +240,8 @@ rand(r::Union{RandomDevice,MersenneTwister}, ::Type{Float32}) =
 
 ## random integers
 
-@inline rand_ui52(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2)) & 0x000fffffffffffff
+@inline rand_ui52_raw(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2))
+@inline rand_ui52(r::AbstractRNG) = rand_ui52_raw(r) & 0x000fffffffffffff
 
 # MersenneTwister
 
@@ -1320,41 +1320,62 @@ randsubseq!(S::AbstractArray, A::AbstractArray, p::Real) = randsubseq!(GLOBAL_RN
 randsubseq{T}(r::AbstractRNG, A::AbstractArray{T}, p::Real) = randsubseq!(r, T[], A, p)
 randsubseq(A::AbstractArray, p::Real) = randsubseq(GLOBAL_RNG, A, p)
 
+"Return a random `Int` (masked with `mask`) in ``[0, n)``, when `n <= 2^52`."
+@inline function rand_lt(r::AbstractRNG, n::Int, mask::Int=nextpow2(n)-1)
+    # this duplicates the functionality of RangeGenerator objects,
+    # to optimize this special case
+    while true
+        x = (rand_ui52_raw(r) % Int) & mask
+        x < n && return x
+    end
+end
 
 function shuffle!(r::AbstractRNG, a::AbstractVector)
-    for i = length(a):-1:2
-        j = rand(r, 1:i)
+    n = length(a)
+    @assert n <= Int64(2)^52
+    mask = nextpow2(n) - 1
+    for i = n:-1:2
+        (mask >> 1) == i && (mask >>= 1)
+        j = 1 + rand_lt(r, i, mask)
         a[i], a[j] = a[j], a[i]
     end
     return a
 end
+
 shuffle!(a::AbstractVector) = shuffle!(GLOBAL_RNG, a)
 
 shuffle(r::AbstractRNG, a::AbstractVector) = shuffle!(r, copy(a))
 shuffle(a::AbstractVector) = shuffle(GLOBAL_RNG, a)
 
 function randperm(r::AbstractRNG, n::Integer)
     a = Array(typeof(n), n)
+    @assert n <= Int64(2)^52
     if n == 0
        return a
     end
     a[1] = 1
-    @inbounds for i = 2:n
-        j = rand(r, 1:i)
+    mask = 3
+    @inbounds for i = 2:Int(n)
+        j = 1 + rand_lt(r, i, mask)
         a[i] = a[j]
         a[j] = i
+        i == 1+mask && (mask = 2mask + 1)
     end
     return a
 end
 randperm(n::Integer) = randperm(GLOBAL_RNG, n)
 
 function randcycle(r::AbstractRNG, n::Integer)
     a = Array(typeof(n), n)
+    n == 0 && return a
+    @assert n <= Int64(2)^52
     a[1] = 1
-    @inbounds for i = 2:n
-        j = rand(r, 1:i-1)
+    mask = 3
+    @inbounds for i = 2:Int(n)
+        j = 1 + rand_lt(r, i-1, mask)
         a[i] = a[j]
         a[j] = i
+        i == 1+mask && (mask = 2mask + 1)
     end
     return a
 end

test/sparsedir/sparse.jl

@@ -951,7 +951,7 @@ A = sparse(tril(rand(5,5)))
 srand(1234321)
 A = triu(sprand(10,10,0.2))       # symperm operates on upper triangle
 perm = randperm(10)
-@test symperm(A,perm).colptr == [1,2,3,3,3,4,5,5,7,9,10]
+@test symperm(A,perm).colptr == [1,1,2,4,5,5,6,8,8,9,10]
 
 # droptol
 @test Base.droptol!(A,0.01).colptr == [1,1,1,2,2,3,4,6,6,7,9]
@@ -1099,7 +1099,8 @@ Ac = sprandn(20,20,.5) + im* sprandn(20,20,.5)
 Ar = sprandn(20,20,.5)
 @test cond(A,1) == 1.0
 @test_approx_eq_eps cond(Ar,1) cond(full(Ar),1) 1e-4
-@test_approx_eq_eps cond(Ac,1) cond(full(Ac),1) 1e-4
+# this test fails sometimes, cf. issue #14778
+# @test_approx_eq_eps cond(Ac,1) cond(full(Ac),1) 1e-4
 @test_approx_eq_eps cond(Ar,Inf) cond(full(Ar),Inf) 1e-4
 @test_approx_eq_eps cond(Ac,Inf) cond(full(Ac),Inf) 1e-4
 @test_throws ArgumentError cond(A,2)