Mergesort implementation in Julia
I'm trying to implement the merge sort algorithm in Julia, but I cannot seem to understand the recursion step needed for the algorithm. My code is the following:
mₐ = [1, 10, 7, 4, 3, 6, 8, 2, 9]
b₁(t, z, half₁, half₂)= ((t<=length(half₁)) && (z<=length(half₂))) && (half₁[t]<half₂[z])
b₂(t, z, half₁, half₂)= ((z<=length(half₂)) && (t<=length(half₁)) ) && (half₁[t]>half₂[z])
function Merge(m₁, m₂)
N = length(m₁) + length(m₂)
B = zeros(N)
i = 1
j = 1
for k in 1:N
if b₁(i, j, m₁, m₂)
B[k] = m₁[i]
i += 1
elseif b₂(i, j, m₁, m₂)
B[k] = m₂[j]
j += 1
elseif j >= length(m₂)
B[k] = m₁[i]
i += 1
elseif i >= length(m₁)
B[k] = m₂[j]
j += 1
end
end
return B
end
function MergeSort(M)
if length(M) == 1
return M
elseif length(M) == 0
return nothing
end
n = length(M)
i₁ = n ÷ 2
i₂ = n - i₁
h₁ = M[1:i₁]
h₂ = M[i₂:end]
C = MergeSort(h₁)
D = MergeSort(h₂)
return Merge(C, D)
end
MergeSort(mₐ)
It always gets stuck when C becomes a single element because it returns it and then splits it again, the only solution is to make it a loop once it is a single element. However, this would not be a recursive approach.
Solution
Taking @Sundar R answer and suggestions. This is a working implementation
#implementation of MergeSort in julia
# merge function, it joins two ordered arrays and returning one single ordered array
function merge(m₁, m₂)
N = length(m₁) + length(m₂)
# create a zeros array of the same input type (int64)
B = zeros(eltype(m₁), N)
i = 1
j = 1
for k in 1:N
if !checkbounds(Bool, m₁, i)
B[k] = m₂[j]
j += 1
elseif !checkbounds(Bool, m₂, j)
B[k] = m₁[i]
i += 1
elseif m₁[i]<m₂[j]
B[k] = m₁[i]
i += 1
else
B[k] = m₂[j]
j += 1
end
end
return B
end
# merge mergesort, this function recursively orders m/2 sub array given an array M
function mergeSort(M)
# base cases
if length(M) == 1
return M
elseif length(M) == 0
return nothing
end
# dividing array in two
n = length(M)
i₁ = n ÷ 2
# be careful with the indexes, thank you @Sundar R
i₂ = i₁ + 1
h₁ = M[1:i₁]
h₂ = M[i₂:end]
# recursively sorting the array
C = mergeSort(h₁)
D = mergeSort(h₂)
return merge(C, D)
end
#test the function
mₐ = [1, 10, 7, 4, 3, 6, 8, 2, 9]
b = mergeSort(mₐ)
println(b)
The issue is with the indices used for splitting, specifically i₂. n - i₁ is the number of elements in the second half of the array, but not necessarily the index where the second half starts - for that you just want i₂ = i₁ + 1.
With i₂ = n - i₁, when n is 2 i.e. when you come down to [1, 10] as the array to sort, i₁ = n ÷ 2 is 1, and i₂ is (2 - 1) = 1 also. So instead of splitting it into [1], [10], you end up "splitting" it into [1], and [1, 10], hence the infinite looping.
Once you fix that, there's a BoundsError from Merge because of a minor mistake: the elseif conditions should check for >, not >= (since Julia uses 1-based indexing, j is still a valid index when j == length(m₂)).
Some other suggestions:
-
zeros(N)returns aFloat64array, so the result here will always be a float array. I'd suggestzeros(eltype(m₁), N)instead. - It feels like
b₁andb₂only complicate the code and make it less clear, I'd suggest a simple nestedifthere, an outer one to check the indices - look upcheckbounds, for eg.checkbounds(Bool, m₁, i)- and an inner one to see which is greater. - Julia convention is to use lowercase for functions, so
mergeandmergesortinstead ofMergeandMergeSort
To add to the previous answers, which deal with some of the problems in your existing code, here is for reference a relatively efficient and straightforward Julia implementation of mergesort:
# Top-level function will allocate temporary arrays for convenience
function mergesort(A)
S = similar(A)
return mergesort!(copy(A), S)
end
# Efficient in-place version
# S is a temporary working (scratch) array
function mergesort!(A, S, n=length(A))
width = 1
swapcount = 0
while width < n
# A is currently full of sorted runs of length `width` (starting with width=1)
for i = 1:2*width:n
# Merge two sorted lists, left and right:
# left = A[i:i+width-1], right = A[i+width:i+2*width-1]
merge!(A, i, min(i+width, n+1), min(i+2*width, n+1), S)
end
# Swap the pointers of `A` and `S` such that `A` now contains merged
# runs of length 2*width.
S,A = A,S
swapcount += 1
# Double the width and continue
width *= 2
end
# Optional, if it is important that `A` be sorted in-place:
if isodd(swapcount)
# If we've swapped A and S an odd number of times, copy `A` back to `S`
# since `S` will by now refer to the memory initially provided as input
# array `A`, which the user will expect to have been sorted in-place
copyto!(S,A)
end
return A
end
# Merge two sorted subarrays, left and right:
# left = A[iₗ:iᵣ-1], right = A[iᵣ:iₑ-1]
@inline function merge!(A, iₗ, iᵣ, iₑ, S)
left, right = iₗ, iᵣ
@inbounds for n = iₗ:(iₑ-1)
if (left < iᵣ) && (right >= iₑ || A[left] <= A[right])
S[n] = A[left]
left += 1
else
S[n] = A[right]
right += 1
end
end
end
This is enough to get us in the same ballpark as Base's implementation of the same algorithm
julia> using BenchmarkTools
julia> @benchmark mergesort!(A,B) setup = (A = rand(50); B = similar(A))
BenchmarkTools.Trial: 10000 samples with 194 evaluations.
Range (min … max): 497.062 ns … 1.294 μs ┊ GC (min … max): 0.00% … 0.00%
Time (median): 501.438 ns ┊ GC (median): 0.00%
Time (mean ± σ): 526.171 ns ± 49.011 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
█▅ ▁ ▁ ▃▇▄ ▁ ▂
█████▇▇▆▇█▇████▇▅▆▅▅▅▆█▆██▄▅▅▄▆██▆▆▄▄▆██▅▃▄██▄▅▅▃▃▃▃▄▅▁▄▄▃▁█ █
497 ns Histogram: log(frequency) by time 718 ns <
Memory estimate: 0 bytes, allocs estimate: 0.
julia> issorted(mergesort(rand(50)))
true
julia> issorted(mergesort(rand(10_000)))
true
julia> @benchmark Base.sort!(A, alg=MergeSort) setup=(A = rand(50))
BenchmarkTools.Trial: 10000 samples with 216 evaluations.
Range (min … max): 344.690 ns … 11.294 μs ┊ GC (min … max): 0.00% … 95.73%
Time (median): 352.917 ns ┊ GC (median): 0.00%
Time (mean ± σ): 401.700 ns ± 378.399 ns ┊ GC (mean ± σ): 3.57% ± 3.76%
█▇▄▄▄▂▁▂▁▂▃▁▁ ▃▂ ▁ ▁▁ ▁
████████████████▇██████▆▆▆▅▆▆▆▆▅▃▅▅▄▅▃▅▅▄▆▅▄▅▄▅▃▄▄██▇▅▆▆▇▆▄▅▅ █
345 ns Histogram: log(frequency) by time 741 ns <
Memory estimate: 336 bytes, allocs estimate: 3.
though both cost a good bit more in terms of both time and memory (the latter due to the need for the working array) in most numeric cases than a similarly efficient pure-Julia implementation of quicksort!:
julia> @benchmark VectorizedStatistics.quicksort!(A) setup = (A = rand(50))
BenchmarkTools.Trial: 10000 samples with 993 evaluations.
Range (min … max): 28.854 ns … 175.821 ns ┊ GC (min … max): 0.00% … 0.00%
Time (median): 35.268 ns ┊ GC (median): 0.00%
Time (mean ± σ): 38.703 ns ± 7.478 ns ┊ GC (mean ± σ): 0.00% ± 0.00%
▂ ▃█▁ ▃▃ ▃▆▂ ▂ ▃ ▂ ▁ ▂ ▂
█▆▃▅▁▁▄▅███▆███▆▆███▁▇█▇▅▇█▆▇█▁▆▅▃▅▄▄██▅▆▅▇▅▄▃▁▄▃▁▄▁▃▃▃▁▄▄▇█ █
28.9 ns Histogram: log(frequency) by time 68.7 ns <
Memory estimate: 0 bytes, allocs estimate: 0.