Mergesort implementation in Julia

Question

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)

Answer 1

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.

Answer 2

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:

1. zeros(N) returns a Float64 array, so the result here will always be a float array. I'd suggest zeros(eltype(m₁), N) instead.
2. It feels like b₁ and b₂ only complicate the code and make it less clear, I'd suggest a simple nested if there, an outer one to check the indices - look up checkbounds, for eg. checkbounds(Bool, m₁, i) - and an inner one to see which is greater.
3. Julia convention is to use lowercase for functions, so merge and mergesort instead of Merge and MergeSort