diff --git a/NEWS.md b/NEWS.md
index 748dda8738e25..23f92e53f5351 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -182,9 +182,12 @@ Standard library changes
 * `randn!(::MersenneTwister, ::Array{Float64})` is faster, and as a result, for a given state of the RNG,
   the corresponding generated numbers have changed ([#35078]).
 
+* `rand!(::MersenneTwister, ::Array{Bool})` is faster, and as a result, for a given state of the RNG,
+  the corresponding generated numbers have changed ([#33721]).
+
 * A new faster algorithm ("nearly division less") is used for generating random numbers
-  within a range ([#29240]).
-  As a result, the streams of generated numbers is changed.
+  within a range ([#29240]). As a result, the streams of generated numbers are changed
+  (for ranges, like in `rand(1:9)`, and for collections in general, like in `rand([1, 2, 3])`).
   Also, for performance, the undocumented property that, given a seed and `a, b` of type `Int`,
   `rand(a:b)` produces the same stream on 32 and 64 bits architectures, is dropped.
 
diff --git a/stdlib/Random/src/RNGs.jl b/stdlib/Random/src/RNGs.jl
index e95e5cc20145c..3ac2c3394a51c 100644
--- a/stdlib/Random/src/RNGs.jl
+++ b/stdlib/Random/src/RNGs.jl
@@ -586,8 +586,10 @@ function rand!(r::MersenneTwister, A::UnsafeView{UInt128}, ::SamplerType{UInt128
 end
 
 for T in BitInteger_types
-    @eval rand!(r::MersenneTwister, A::Array{$T}, sp::SamplerType{$T}) =
-        (GC.@preserve A rand!(r, UnsafeView(pointer(A), length(A)), sp); A)
+    @eval function rand!(r::MersenneTwister, A::Array{$T}, sp::SamplerType{$T})
+        GC.@preserve A rand!(r, UnsafeView(pointer(A), length(A)), sp)
+        A
+    end
 
     T == UInt128 && continue
 
@@ -602,6 +604,52 @@ for T in BitInteger_types
     end
 end
 
+
+#### arrays of Bool
+
+# similar to Array{UInt8}, but we need to mask the result so that only the LSB
+# in each byte can be non-zero
+
+function rand!(r::MersenneTwister, A1::Array{Bool}, sp::SamplerType{Bool})
+    n1 = length(A1)
+    n128 = n1 ÷ 16
+
+    if n128 == 0
+        bits = rand(r, UInt52Raw())
+    else
+        GC.@preserve A1 begin
+            A = UnsafeView{UInt128}(pointer(A1), n128)
+            rand!(r, UnsafeView{Float64}(A.ptr, 2*n128), CloseOpen12())
+            # without masking, non-zero bits could be observed in other
+            # positions than the LSB of each byte
+            mask = 0x01010101010101010101010101010101
+            # we need up to 15 bits of entropy in `bits` for the final loop,
+            # which we will extract from x = A[1] % UInt64;
+            # let y = x % UInt32; y contains 32 bits of entropy, but 4
+            # of them will be used for A[1] itself (the first of
+            # each byte). To compensate, we xor with (y >> 17),
+            # which gets the entropy from the second bit of each byte
+            # of the upper-half of y, and sets it in the first bit
+            # of each byte of the lower half; the first two bytes
+            # now contain 16 usable random bits
+            x = A[1] % UInt64
+            bits = x ⊻ x >> 17
+            for i = 1:n128
+                # << 5 to randomize the first bit of the 8th & 16th byte
+                # (i.e. we move bit 52 (resp. 52 + 64), which is unused,
+                # to position 57 (resp. 57 + 64))
+                A[i] = (A[i] ⊻ A[i] << 5) & mask
+            end
+        end
+    end
+    for i = 16*n128+1:n1
+        @inbounds A1[i] = bits % Bool
+        bits >>= 1
+    end
+    A1
+end
+
+
 ### randjump
 
 # Old randjump methods are deprecated, the scalar version is in the Future module.