fix #9037: segfault caused by incorrect alignment assumptions

The rand!(::MersenneTwister, A::Array{Float64}) assumed that the
array A was 16-byte aligned, which can be false (e.g. after a call
to resize!) and lead to a segfault in libdSFMT.

base/dSFMT.jl

@@ -39,19 +39,18 @@ function dsfmt_init_by_array(s::DSFMT_state, seed::Vector{UInt32})
           s.val, seed, length(seed))
 end
 
-
-# precondition for dsfmt_fill_array_*:
-# the underlying C array must be 16-byte aligned, which is the case for "Array"
-function dsfmt_fill_array_close1_open2!(s::DSFMT_state, A::Array{Float64}, n::Int)
-    @assert dsfmt_min_array_size <= n <= length(A)  && iseven(n)
+function dsfmt_fill_array_close1_open2!(s::DSFMT_state, A::Ptr{Float64}, n::Int)
+    @assert Int(A) % 16 == 0 # the underlying C array must be 16-byte aligned
+    @assert dsfmt_min_array_size <= n && iseven(n)
     ccall((:dsfmt_fill_array_close1_open2,:libdSFMT),
           Void,
           (Ptr{Void}, Ptr{Float64}, Int),
           s.val, A, n)
 end
 
-function dsfmt_fill_array_close_open!(s::DSFMT_state, A::Array{Float64}, n::Int)
-    @assert dsfmt_min_array_size <= n <= length(A)  && iseven(n)
+function dsfmt_fill_array_close_open!(s::DSFMT_state, A::Ptr{Float64}, n::Int)
+    @assert Int(A) % 16 == 0 # the underlying C array must be 16-byte aligned
+    @assert dsfmt_min_array_size <= n && iseven(n)
     ccall((:dsfmt_fill_array_close_open,:libdSFMT),
           Void,
           (Ptr{Void}, Ptr{Float64}, Int),

base/random.jl

@@ -33,7 +33,7 @@ end
 @inline mt_pop!(r::MersenneTwister) = @inbounds return r.vals[r.idx+=1]
 
 function gen_rand(r::MersenneTwister)
-    dsfmt_fill_array_close1_open2!(r.state, r.vals, length(r.vals))
+    dsfmt_fill_array_close1_open2!(r.state, pointer(r.vals), length(r.vals))
     mt_setfull!(r)
 end
 
@@ -224,15 +224,32 @@ end
 
 rand!(r::MersenneTwister, A::AbstractArray{Float64}) = rand_AbstractArray_Float64!(r, A)
 
-fill_array!(s::DSFMT_state, A::Array{Float64}, n::Int, ::Type{CloseOpen}) = dsfmt_fill_array_close_open!(s, A, n)
-fill_array!(s::DSFMT_state, A::Array{Float64}, n::Int, ::Type{Close1Open2}) = dsfmt_fill_array_close1_open2!(s, A, n)
+fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{CloseOpen}) = dsfmt_fill_array_close_open!(s, A, n)
+fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{Close1Open2}) = dsfmt_fill_array_close1_open2!(s, A, n)
 
 function rand!{I<:FloatInterval}(r::MersenneTwister, A::Array{Float64}, n=length(A), ::Type{I}=CloseOpen)
-    if n < dsfmt_get_min_array_size()
+    # depending on the alignment of A, the data written by fill_array! may have
+    # to be left-shifted by up to 15 bytes (cf. unsafe_copy! below) for
+    # reproducibility purposes;
+    # so, even for well aligned arrays, fill_array! is used to generate only
+    # the n-2 first values (or n-3 if n is odd), and the remaining values are
+    # generated by the scalar version of rand
+    n2 = (n-2) ÷ 2 * 2
+    if n2 < dsfmt_get_min_array_size()
         rand_AbstractArray_Float64!(r, A, n, I)
     else
-        fill_array!(r.state, A, 2*(n ÷ 2), I)
-        isodd(n) && (A[n] = rand(r, I))
+        pA = pointer(A)
+        align = Int(pA) % 16
+        if align > 0
+            pA2 = pA + 16 - align
+            fill_array!(r.state, pA2, n2, I)
+            unsafe_copy!(pA, pA2, n2)
+        else
+            fill_array!(r.state, pA, n2, I)
+        end
+        for i=n2+1:n
+            A[i] = rand(r, I)
+        end
     end
     A
 end

test/random.jl

@@ -238,3 +238,19 @@ let mt = MersenneTwister(0)
     resize!(AF64, 2*length(mt.vals))
     @test Base.Random.rand_AbstractArray_Float64!(mt, AF64)[end]  == 0.432757268470779
 end
+
+# Issue #9037
+let mt = MersenneTwister()
+    a = Array(Float64, 0)
+    resize!(a, 1000) # could be 8-byte aligned
+    b = Array(Float64, 1000) # should be 16-byte aligned
+    c8 = Array(UInt8, 8001)
+    c = pointer_to_array(Ptr{Float64}(pointer(c8, 2)), 1000) # Int(pointer(c)) % 16 == 1
+
+    for A in (a, b, c)
+        srand(mt, 0)
+        rand(mt) # this is to fill mt.vals, cf. #9040
+        rand!(mt, A) # must not segfault even if Int(pointer(A)) % 16 != 0
+        @test A[end-4:end] == [0.49508297796349776,0.3408340446375888,0.3211229457075784,0.9103565379264364,0.16456579813368521]
+    end
+end