Magic squares in Julia

magic.jl

# These are possibly naïve implementations of magic squares generation
# for a micro-benchmark to compare Julia to Python and Matlab
#

function magic_matlab(n::Int64)
	# Works exactly as Matlab's magic.m

	if n % 2 == 1
		p = (1:n)
		M = n * mod(broadcast(+, p', p - div(n+3, 2)), n) + mod(broadcast(+, p', 2p - 2), n) + 1
		return M
	elseif n % 4 == 0
		J = div([1:n] % 4, 2)
		K = J' .== J
		M = broadcast(+, [1:n:(n*n)]', [0:n-1])
		M[K] = n^2 + 1 - M[K]
		return M
	else
		p = div(n, 2)
		M = magic_matlab(p)
		M = [M M+2p^2; M+3p^2 M+p^2]
		if n == 2
			return M
		end
		i = (1:p)
		k = (n-2)/4
		j = convert(Array{Int}, [(1:k); ((n-k+2):n)])
		M[[i; i+p],j] = M[[i+p; i],j]
		i = k+1
		j = [1; i]
		M[[i; i+p],j] = M[[i+p; i],j]
		return M
	end
end
@vectorize_1arg Int magic_matlab

function magic_python(n::Int64)
	# Works exactly as magic_square.py (from pypy)

	if n % 2 == 1
		m = (n >> 1) + 1
		b = n^2 + 1

		M = reshape(repmat(1:n:b-n, 1, n+2)[m:end-m], n+1, n)[2:end, :] +
		    reshape(repmat(0:(n-1), 1, n+2), n+2, n)[2:end-1, :]'
		return M
	elseif n % 4 == 0
		b = n^2 + 1
		d = reshape(1:b-1, n, n)

		d[1:4:n, 1:4:n] = b - d[1:4:n, 1:4:n]
		d[1:4:n, 4:4:n] = b - d[1:4:n, 4:4:n]
		d[4:4:n, 1:4:n] = b - d[4:4:n, 1:4:n]
		d[4:4:n, 4:4:n] = b - d[4:4:n, 4:4:n]
		d[2:4:n, 2:4:n] = b - d[2:4:n, 2:4:n]
		d[2:4:n, 3:4:n] = b - d[2:4:n, 3:4:n]
		d[3:4:n, 2:4:n] = b - d[3:4:n, 2:4:n]
		d[3:4:n, 3:4:n] = b - d[3:4:n, 3:4:n]

		return d
	else
		m = n >> 1
		k = m >> 1
		b = m^2

		d = repmat(magic_python(m), 2, 2)

		d[1:m, 1:k] += 3*b
		d[1+m:end, 1+k:m] += 3*b
		d[1+k, 1+k] += 3*b
		d[1+k, 1] -= 3*b
		d[1+m+k, 1] += 3*b
		d[1+m+k, 1+k] -= 3*b
		d[1:m,1+m:n-k+1] += b+b
		d[1+m:end, 1+m:n-k+1] += b
		d[1:m, 1+n-k+1:end] += b
		d[1+m:end, 1+n-k+1:end] += b+b

		return d
   end
end
@vectorize_1arg Int magic_python

print("Matlab version: ")
@time magic_matlab(3:1000)

print("Python version: ")
@time magic_python(3:1000)