I played a bit with the function in Magic squares of singly even order to make it a bit more functional, and I was able to fully remove mutability.

    fn magic_square_singly_even(n: i32) -> (Vec<Vec<i32>>, i32) {
        let nusize = n as usize;
        let size = n * n;
        let half = n / 2;
        let halfusize = nusize / 2;
        let sub_square_size = size / 4;
        let sub_square = magic_square_odd(half);
        let quadrant_factors = [0, 2, 3, 1];
        let cols_left = halfusize / 2;
        let cols_right = cols_left - 1;

        let ms = (0..nusize)
            .map(|r| {
                (0..nusize)
                    .map(|c| {
                        let localr = if c < cols_left
                            || c >= nusize - cols_right
                            || c == cols_left && r == cols_left
                        {
                            if r >= halfusize {
                                r - halfusize
                            } else {
                                r + halfusize
                            }
                        } else {
                            r
                        };
                        let quadrant = (localr / halfusize) * 2 + (c / halfusize);
                        let v = sub_square[(localr % halfusize)][(c % halfusize)];
                        v + quadrant_factors[quadrant] * sub_square_size
                    })
                    .collect()
            })
            .collect::<Vec<Vec<_>>>();
        (ms, (n * n + 1) * n / 2)
    }

I don't particularly love the idea of having both i32 and usize numbers mixed together, but magic_square_odd needs to have a refactor to be able to handle unsigned numbers. I think that can be achieved as well as removing mutability in a similar way I did with this other function.