magic.jl

function magic!(M::Array{Int}, n::Int)
    # Fast magic square generation
    #
    #
    if n % 2 == 1
        # Odd-order magic square
        # Algorithm: http://en.wikipedia.org/wiki/Magic_square#Method_for_constructing_a_magic_square_of_odd_order
        #
        ndiv2 = div(n, 2);
        for I = 1:n
            intermediate1 = (I-1+ndiv2) % n
            intermediate2 = (2I - 2) % n
            @simd for J = 1:n
                # Using these instead of modulo operations gives a huge speed improvement
                # by Simon Kornblith, https://gist.github.com/simonster/6195af68c6df33ca965d
                x = ifelse(J + intermediate1 < n, J + intermediate1, J + intermediate1 - n)
                y = ifelse(J + intermediate2 < n, J + intermediate2, J + intermediate2 - n)
                @inbounds M[J, I] = n*x + y + 1
            end
        end
    elseif n % 4 == 0
        # Doubly-even order magic square
        # Algorithm: http://en.wikipedia.org/wiki/Magic_square#A_method_of_constructing_a_magic_square_of_doubly_even_order
        #            http://hal.archives-ouvertes.fr/docs/00/94/51/29/PDF/magicsquarev2.pdf
        # by Iain Dunning, https://gist.github.com/IainNZ/9b5f1eb1bcf923ed02d9
        for I = 1:n
            In = (I-1)*n
            @simd for J = 1:n
                @inbounds M[I,J] = In+J
            end
        end
        for outer_col = 1:4:n, outer_row = 1:4:n
            # Sides
            for i in (1,2), j in (0,3)
                @inbounds M[outer_row+i,outer_col+j] = n^2 - M[outer_row+i,outer_col+j]
            end
            # Top and bottom
            for i in (0,3), j in (1,2)
                @inbounds M[outer_row+i,outer_col+j] = n^2 - M[outer_row+i,outer_col+j]
            end
        end
    else
        # Even order magic square
        # Translated from https://pypi.python.org/pypi/magic_square/0.2
        # and devectorized by suggestion of the Julia-users group
        # by Phillip Berndt
        m = n >> 1
        k = m >> 1
        b = m*m
        b3 = 3b

        magic!(M, m)
        for J = 1:m
            Jm = J+m
            @simd for I = 1:m
                @inbounds M[I+m, J]   = M[I, J]
                @inbounds M[I,   Jm]  = M[I, J]
                @inbounds M[I+m, Jm] = M[I, J]
            end
        end

        for J = 1:n
            @simd for I = 1:n
                if I <= k
                    if J <= m
                        @inbounds M[I, J] += b3
                    end
                elseif I <= m
                    if 1+m <= J
                        @inbounds M[I, J] += b3
                   end
                elseif I <= n-k + 1
                    if J <= m
                        @inbounds M[I, J] += b + b
                    else
                        @inbounds M[I, J] += b
                    end
                else
                    if J <= m
                        @inbounds M[I, J] += b
                    else
                        @inbounds M[I, J] += b + b
                    end
                end
            end
        end

        M[1+k, 1+k]   += b3
        M[1, 1+k]     -= b3
        M[1, 1+m+k]   += b3
        M[1+k, 1+m+k] -= b3
    end
end

# Speed test:
@time begin
    for i = 3:1:1000
        M = Array(Int, i, i)
        magic!(M, i);
    end
end

xxcppversion.cpp

// Direct translation of the current version.
// Not nice code, sorry ;) Only for testing..

#include <time.h>
#include <stdio.h>

struct mat {
    int n;
    int *mem;

    mat(int n) {
        this->n = n;
        this->mem = new int[n*n];
    }

    ~mat() {
        delete this->mem;
    }

    inline int& operator()(const int row, const int col) {
        return this->mem[this->n * (col - 1) + row - 1];
    }
};

void iain_magic(mat &M, int n) {
    if(n % 2 == 1) {
        int ndiv2 = n / 2;
        for(int I = 1; I<=n; I++) {
            int intermediate1 = (I-1+ndiv2 % n);
            int intermediate2 = (2*I - 2) % n;
            for(int J = 1; J<=n; J++) {
                int x = (J + intermediate1 < n? J + intermediate1: J + intermediate1 - n);
                int y = (J + intermediate2 < n? J + intermediate2: J + intermediate2 - n);
                M(I, J) = n*x + y + 1;
				// M(I, J) = n*((J + intermediate1)%n) + ((J + intermediate2)%n) + 1;
            }
        }
    }
    else if(n % 4 == 0) {
        for(int I = 1; I <= n; I++) {
            int In = (I-1)*n;
            for(int J = 1; J<=n; J++) {
                M(I,J) = In+J;
            }
        }
        for(int outer_col = 1; outer_col <= n; outer_col += 4) {
            for(int outer_row = 1; outer_row <= n; outer_row += 4) {
                for(int i= 1; i<= 2; i++) {
                    for(int j=0; j<= 3; j++) {
                        M(outer_row+i,outer_col+j) = n*n - M(outer_row+i,outer_col+j);
                    }
                }
                for(int i= 0; i<= 3; i++) {
                    for(int j=1; j<= 2; j++) {
                        M(outer_row+i,outer_col+j) = n*n - M(outer_row+i,outer_col+j);
                    }
                }
            }
        }
    }
    else {
        int m = n >> 1;
        int k = m >> 1;
        int b = m*m;
        int b3 = 3*b;

        iain_magic(M, m);
        for(int J = 1; J<=m; J++) {
            int Jm = J+m;
            for(int I = 1; I <=m ; I++) {
                M(I+m, J)   = M(I, J);
                M(I,   Jm)  = M(I, J);
                M(I+m, Jm) = M(I, J);
            }
       }

       for(int J=1; J<=n; J++) {
            for(int I=1; I<=n; I++) {
                if(I <= k) {
                    if(J <= m)
                        M(I, J) += b3;
                    else if(I <= m)
                        if(1+m <= J) M(I, J) += b3;
                } else if(I <= n-k + 1) {
                    if(J <= m)
                        M(I, J) += b + b;
                    else
                        M(I, J) += b;
                } else {
                    if(J <= m)
                        M(I, J) += b;
                    else
                        M(I, J) += b + b;
                }
            }
       }

        M(1+k, 1+k)   += b3;
        M(1, 1+k)     -= b3;
        M(1, 1+m+k)   += b3;
        M(1+k, 1+m+k) -= b3;
    }
}

int main(int argc, char *argv[]) {
    clock_t start = clock();


    for(int i=3; i<1000; i++) {
        mat A(i);
        iain_magic(A, i);
    }
    printf("Needed %f s\n", (clock() - start) * 1. / CLOCKS_PER_SEC);
}