[Question]
  [
Originally from [a CMC I proposed for the last BMG event](https://chat.stackexchange.com/transcript/message/58426955#58426955)
## Challenge
Given a non-negative integer \$n\$, create a 2D array of size \$2^n √ó 2^n\$ which is generated in the following manner:
1. Divide the matrix into four quadrants of size \$2^{n-1} √ó 2^{n-1}\$.
2. Visiting order of the four quadrants is defined to be the Z-shape (top-left, top-right, bottom-left, then bottom-right).
3. Recursively apply the ordering (steps 1-2) to each quadrant, until the ordering is defined for each cell in the matrix.
4. Visit each cell in the defined order, sequentially writing down 0, 1, 2, 3, ... to each cell.
You can output 1-based instead of 0-based (add 1 to all cells in the examples below).
Standard [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'") rules apply. The shortest code in bytes wins.
## Examples
```
n = 0:
[[0]]
n = 1:
[[0, 1],
 [2, 3]]
n = 2:
[[0,  1,  4,  5],
 [2,  3,  6,  7],
 [8,  9,  12, 13],
 [10, 11, 14, 15]]
n = 3:
[[0,  1,  4,  5,  16, 17, 20, 21],
 [2,  3,  6,  7,  18, 19, 22, 23],
 [8,  9,  12, 13, 24, 25, 28, 29],
 [10, 11, 14, 15, 26, 27, 30, 31],
 [32, 33, 36, 37, 48, 49, 52, 53],
 [34, 35, 38, 39, 50, 51, 54, 55],
 [40, 41, 44, 45, 56, 57, 60, 61],
 [42, 43, 46, 47, 58, 59, 62, 63]]
```
Brownie points for beating or tying with my ~~9~~ 6 bytes in Jelly or 19 bytes in J.
  
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 5 bytes
```
4·πó>2ƒ†
```
[Try it online!](https://tio.run/##y0rNyan8/9/k4c7pdkZHFvw/3O7@/78xAA "Jelly ‚Äì Try It Online") Neil saved a byte. Thanks!
```
4·πó       n-th Cartesian power of [1,2,3,4]
  >2     Replace 1,2,3,4 with 0,0,1,1
    ƒ†    Group indices by values
```
For n=2, `4·πó>2` generates a list of 16 pairs:
```
 1: [0, 0]
 2: [0, 0]
 3: [0, 1]
 4: [0, 1]
 5: [0, 0]
 6: [0, 0]
 7: [0, 1]
 8: [0, 1]
 9: [1, 0]
10: [1, 0]
11: [1, 1]
12: [1, 1]
13: [1, 0]
14: [1, 0]
15: [1, 1]
16: [1, 1]
```
Then `ƒ†` turns this into a list of lists:
* all indices where the value is `0 0`: namely `[ 1, 2, 5, 6]`,
* all indices where the value is `0 1`: namely `[ 3, 4, 7, 8]`,
* all indices where the value is `1 0`: namely `[ 9, 10, 13, 14]`,
* all indices where the value is `1 1`: namely `[11, 12, 15, 16]`.
The result is a Z-matrix.
## Why does this work?
Suppose we bump the input up to \$n=3\$. How does our output compare to that for \$n=2\$?
Let's call `4·πó>2` a function \$f\$. The list of 16 pairs above is the output of \$f(2)\$. By definition of the Cartesian power, \$f(3) = [0,0,1,1] \times f(2)\$ =
```
  [[0] + x for x in f(2)]     (indices  1..16)
+ [[0] + x for x in f(2)]     (indices 17..32)
+ [[1] + x for x in f(2)]     (indices 33..48)
+ [[1] + x for x in f(2)]     (indices 49..64)
```
There are now 8 values for `ƒ†` to group the indices by, ranging from `[0,0,0]` to `[1,1,1]`.
All indices in `1..32` end up in the first four groups, and because we repeated `[[0] + x for x in f(2)]` twice, `a[i+16]` is always in the same group as `a[i]`. Then, all indices in `33..64` end up in the last four groups, and their relative order is the same as the first half of the list (the only difference is they have `1` instead of `0` in front). Within each quadrant, the order is the same as it was in `f(2)`. So we get:
```
[0, 0, 0]:  1  2  5  6 17 18 21 22
[0, 0, 1]:  3  4  7  8 19 20 23 24
[0, 1, 0]:  9 10 13 14 25 26 29 30
[0, 1, 1]: 11 12 15 16 27 28 31 32
[1, 0, 0]: 33 34 37 38 49 50 53 54
[1, 0, 1]: 35 36 39 40 51 52 55 56
[1, 1, 0]: 41 42 45 46 57 58 61 62
[1, 1, 1]: 43 44 47 48 59 60 63 64
```
[Answer]
# [Charcoal](https://github.com/somebody1234/Charcoal), ~~27~~ ~~23~~ 18 bytes
```
‚âîÔº•Ôº∏¬≤ÔºÆ‚Ü®‚Å¥‚Ü®Œπ¬≤œÖÔº©Ôº•‚äóœÖ‚Å∫ŒπœÖ
```
[Try it online!](https://tio.run/##S85ILErOT8z5/9@xuDgzPU/DN7FAIyC/PLVIw0hHwTOvoLTErzQ3CcjV1NRRcEosTtUwgdKZOgpGmiDRUk1rroCizLwSDefE4hKwCS75pUk5qSkapUDpgJzSYpDiUqBiTev//03@65blAAA "Charcoal ‚Äì Try It Online") Link is to verbose version of code. Explanation: Now a port of @tsh's Python answer.
```
‚âîÔº•Ôº∏¬≤ÔºÆ‚Ü®‚Å¥‚Ü®Œπ¬≤œÖ
```
Map the numbers up to 2‚Åø by converting them to base 2 and then interpreting the result as base 4.
```
Ôº©Ôº•‚äóœÖ‚Å∫ŒπœÖ
```
Vectorised add the doubled array to itself as a matrix.
Previous ~~27~~ 23 byte answer uses a different approach that is interesting in its own right because some golfing languages have built-ins for things like group index by value or sort array by key.
```
‚âîÔº•Ôº∏‚Å¥ÔºÆ‚Ü®¬≤√∑‚Ü®Œπ‚Å¥¬¶¬≤œÖÔº©Ôº•‚äï‚åàœÖ‚åïÔº°œÖŒπ
```
[Try it online!](https://tio.run/##HYw7DsIwEER7TuFyLZkmcpcqgJBSgHIFY69gJXsT@RO4vbFTzps3Yz8m2tX4WqeU6M3wMBss6xcjaCVm3kp@lvBqUUolLiYhDJ3nG@3kEA5CSujWDrI7RY6nJRJnuJqUj7@ZbcSAnNG1/KNQApTu3ond5D0UJaiN5Virrufd/wE "Charcoal ‚Äì Try It Online") Link is to verbose version of code. Explanation:
```
‚âîÔº•Ôº∏‚Å¥ÔºÆ‚Ü®¬≤√∑‚Ü®Œπ‚Å¥¬¶¬≤œÖ
```
For each element in the final diagram, convert it to base 4, halve the digits, then convert from base 2. This gives the row number for that element.
```
Ôº©Ôº•‚äï‚åàœÖ‚åïÔº°œÖŒπ
```
For each row list the elements assigned to that row.
[Answer]
# [J](http://jsoftware.com/), 21 19 bytes
Found Bubbler's solution. :-)
```
2(*+/])4#.2#:@i.@^]
```
[Try it online!](https://tio.run/##y/qvpKeepmBrpaCuoKNgoGAFxLp6Cs5BPm7/jTS0tPVjNU2U9YyUrRwy9RziYv9rcmmkJmfkO1inaSppGGgqZOopmP4HAA "J ‚Äì Try It Online")
The first row is [A000695](http://oeis.org/A000695), the first column is the same but times two. In the comments Marc LeBrun notes that this is rebasing `n` from base 2 into base 4. So that's helpful with J:
```
2(*+/])4#.2#:@i.@^]
          2   i.@^]  0 1 ‚Ä¶ 2^n
           #:        digits in base 2
                      0 0 0
                      0 0 1
                          ‚Ä¶
                      1 1 0
                      1 1 1
       4#.           interpret as digits in base 4
                      0 1 4 5 16 17 20 21
2(*+/])              addition table of list*2 with list
```
[Answer]
# [R](https://www.r-project.org/), ~~65~~ ~~63~~ ~~59~~ ~~57~~ 55 bytes
```
for(i in seq(l=scan()))T=T^4%x%2^rbind(0:1,2:3);log2(T)
```
[Try it online!](https://tio.run/##Vc1BCsIwEEDRvacYkMIMKLTRVbW3yLrQ1kQGyiROIgri2dO6Ebr98PhaPFyP4J8yZQ6CTPCBNA2yqcSwty5leOkQo9PigyIDCyT3wLn7ASQi29n@XL0r0@vIcsO6bQ6mPdFlDneDlsp395f12qOyZPTrlcoC "R ‚Äì Try It Online")
Relies on the observation that each next iteration `N` relates to the previous iteration `P` by Kronecker product-like operation with the matrix `M` for `n=1`, but using summation instead of multiplication:
```
N = kronecker(4*P, M, FUN = "+")
```
Here, we transform the values with powers and logarithms in order to make use of the actual Kronecker product function, which is conveniently abbreviated as `%x%` operator.
[Answer]
# [Python 2](https://docs.python.org/2/), 81 bytes
```
R=[int(bin(x)[2:],4)for x in range(2**input())]
for y in R:print[x+y*2for x in R]
```
[Try it online!](https://tio.run/##K6gsycjPM/r/P8g2OjOvRCMpM0@jQjPayCpWx0QzLb9IoUIhM0@hKDEvPVXDSEsrM6@gtERDUzOWCyRXCZILsiooAuqMrtCu1DKC6wiK/f/fGAA "Python 2 ‚Äì Try It Online")
[Answer]
# [J](http://jsoftware.com/), 33 bytes
```
,~@^~&2$4(/:[:<.@-:4#.inv])@i.@^]
```
[Try it online!](https://tio.run/##y/qvpKeepmBrpaCuoKNgoGAFxLp6Cs5BPm7/deoc4urUjFRMNPStoq1s9Bx0rUyU9TLzymI1HTL1HOJi/2typSZn5CukKRhAGOrqMAFDdAEjdAHj/wA "J ‚Äì Try It Online")
The true insight of this idea is due to [Neil's excellent answer](https://codegolf.stackexchange.com/a/231027/15469).
I'll mention a couple of other insights I had which didn't help with bytes but might be useful to others:
1. The deltas along each row follow a repeating pattern whose unique elements are <http://oeis.org/A007583>: 1, 3, 11, 43, 171...
2. Similarly the deltas going down each column follow <http://oeis.org/A047849>: 1, 2, 6, 22, 86, 342...
## how
* `4...@i.@^]` Produces `0, 1, ... 4^n`, where n is the input.
* `/:[:<.@-:4#.inv]` Sort those according to row number, which we get by converting to base 4, halving, and flooring.
* `,~@^~&2$` Now shape that into a square matrix of side length `2^n`.
[Answer]
# [R](https://www.r-project.org/), ~~57~~ ~~52~~ 49 bytes
*Edit: -8 bytes by changing to iterative approach with `scan()` for input (instead of recursive function)*
```
for(i in seq(l=2*scan())-1)F=t(cbind(F,F+2^i));+F
```
[Try it online!](https://tio.run/##K/r/Py2/SCNTITNPoTi1UCPH1kirODkxT0NTU9dQ0822RCM5KTMvRcNNx03bKC5TU9Na2@2/8X8A "R ‚Äì Try It Online")
**Explanation**
~~Recursive function that~~ Each loop `i`creates two side-by-side copies of the columns of a matrix, adding ~~the maximum value of the first +1~~ `2^i` to the second copy, and then transposes it all.
This creates the top left & top right copies of the last matrix (`F`), transposed: so, every second loop creates a full 'Z' of four copies of the last matrix, untransposed.
So we need to loop twice the value of `n` times.
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 9 bytes
```
4*·∏∂b4:2·∏Ñƒ†
```
[Try it online!](https://tio.run/##y0rNyan8/99E6@GObUkmVkYPd7QcWfD/cLv7///GAA "Jelly ‚Äì Try It Online")
Pretty much a port of [Neil's Charcoal answer](https://codegolf.stackexchange.com/a/231027/68942), so go upvote them too.
```
4*·∏∂b4:2·∏Ñƒ†    Main Link; take x
4*           4 ** x
  ·∏∂          [0, 1, 2, ..., 4 ** x - 1]
   b4        Convert to base 4
     :2      Floor-divide by 2
       ·∏Ñ     Convert from binary
        ƒ†    Group the indices by their corresponding values
```
This answer can be ported directly into yuno.
# [yuno](https://github.com/hyper-neutrino/yuno), 9 bytes
```
4* Äb4:2…ÆG
```
[Try it online!](https://yuno.hyper-neutrino.xyz/#WyIiLCI0KsqAYjQ6MsmuRyIsIiIsIiIsIiIsWyIzIl1d "yuno online interpreter")
---
This approach is probably what Bubbler's 6-byter is going for, but a) I need to find how to get `0,0,1,1` in 3 bytes, and b) I need to figure out how to make the format valid without the `s2*$` being that long:
# [Jelly](https://github.com/DennisMitchell/jelly), 11 bytes
```
√ò.x2¬§·πó·ª§s2*$
```
[Try it online!](https://tio.run/##ARwA4/9qZWxsef//w5gueDLCpOG5l@G7pHMyKiT///8z "Jelly ‚Äì Try It Online")
# [Jelly](https://github.com/DennisMitchell/jelly), 16 bytes
```
,Uz0$FU
2*·∏∂B√ß√æ`·∏Ñ
```
[Try it online!](https://tio.run/##y0rNyan8/18ntMpAxS2Uy0jr4Y5tToeXH96X8HBHy////40B "Jelly ‚Äì Try It Online")
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), ~~49~~ 32 bytes
```
PositionIndex@Tuples[0|0|1|1,#]&
```
[Try it online!](https://tio.run/##y00syUjNTSzJTE78n2b7PyC/OLMkMz/PMy8ltcIhpLQgJ7U42qDGoMawxlBHOVbtf0BRZl6JQ5qDa3JGvoOyWl1wcmJeXVBiXnpqtIGOcex/AA "Wolfram Language (Mathematica) ‚Äì Try It Online")
Returns an `Association` where the values correspond with rows of the matrix (associations are ordered). For a traditional list-of-lists output, prepend `List@@` (+6 bytes).
Port of Lynn's [Jelly solution](https://codegolf.stackexchange.com/a/231062/81203).
---
Older solution:
```
Outer[#+##&,t=#~FromDigits~4&/@{0,1}~Tuples~#,t]&
```
[Try it online!](https://tio.run/##BcGxDkAwEADQj7mkiyYIq6QDVoJNDJemaKIldSZxv17vOaTdOCSrMa5V7B4yYYYEQEiqgNtwutpulm4uRareTOYfT891mJtB0iJiH6ynVK2q0fupQPCo0fOAfjNzJosl/g "Wolfram Language (Mathematica) ‚Äì Try It Online")
[Answer]
# ARM A32 machine code, 48 bytes
```
e3a0c000 e1a0215c e1b03152 112fff1e
e02c3112 ec432b10 f2800e00 f2810590
f2211110 f480180d e28cc001 eafffff4
```
Following the [AAPCS](https://github.com/ARM-software/abi-aa/blob/2bcab1e3b22d55170c563c3c7940134089176746/aapcs32/aapcs32.rst), this takes in r0 the address of an *array* of *array*s of 32-bit integers to be filled with the result, and takes n in r1.
Assembly:
```
.section .text
.code 32
.global zm
zm:
    mov r12, #0               @ Initialise counter to 0
loop:
    asr r2, r12, r1           @ Shift counter right by n, to obtain the top half
    asrs r3, r2, r1           @ Shift that right by n again, to check for ending (2^2n)
    bxne lr                   @ If so, return
    eor r3, r12, r2, lsl r1   @ Exclusive-or cancels out the top half, to obtain the bottom half
    vmov d0, r2, r3           @ Place the two halves into a 64-bit SIMD&FP register
    vmull.p8 q0, d0, d0       @ Multiply by itself as polynomials over GF(2);
               @ because this is a ring of characteristic 2, (a+b)^2 = a^2 + b^2 holds,
               @ thus squaring doubles each exponent in the polynomial, spreading out the bits.
               @ Results go into Q0, which is also D0 and D1.
    vshl.i64 d0, #1           @ Shift the squared top half left by 1 ...
    vorr d1, d0               @ ... and combine by OR with the bottom half
    vst1.32 {d1[0]}, [r0]!    @ Place the result into memory, advancing the pointer
    add r12, #1               @ Increment the counter
    b loop                    @ Repeat
```
(This should work on v7 or later from what I can tell, but in the testing setup I'm using (similar to [here](https://codegolf.stackexchange.com/questions/175236/sum-square-difference/231078#231078)), it gives 'instruction not supported' errors on versions lower than v8.4, which doesn't seem right.)
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), ~~¬†14 12¬†~~ 9 [bytes](https://github.com/DennisMitchell/jelly/wiki/Code-page)
### 9 byter
```
√ò.·πó·∏Ö4+√æ·∏§$
```
**[Try it online!](https://tio.run/##ARwA4/9qZWxsef//w5gu4bmX4biFNCvDvuG4pCT///8z "Jelly ‚Äì Try It Online")**
##### How?
Makes a table under \$x+2y\$ where both \$x\$ and \$y\$ are a prefix of the Moser‚ÄìDe Bruijn sequence.
```
√ò.·πó·∏Ö4+√æ·∏§$ - Link: n
√ò.        - bits -> [0,1]
  ·πó       - Cartesian power (n) -> binary representations of [0..2^n-1]
   ·∏Ö4     - convert from base four -> first 2^n terms of the Moser‚ÄìDe Bruijn sequence
        $ - last two links as a monad, f(S):
       ·∏§  -   double -> D=[2x for x in S]
      √æ   -   (s in S) table with (d in D):
     +    -     (s) add (d)
```
(Still no clue for the 6 byte code though!)
---
### 12 Byter
-1 thanks to [caird coinheringaahing](https://codegolf.stackexchange.com/users/66833/caird-coinheringaahing) (we can use 1-based values)
```
ƒ†;"FL+∆äZ∆ä‚Å∏¬°‚Å∫
```
A monadic Link that accepts a non-negative integer and yields the matrix (using the 1-based values option).
**[Try it online!](https://tio.run/##ASYA2f9qZWxsef//xKA7IkZMK8aKWsaK4oG4wqHigbr/w4fFkuG5mP//Mw "Jelly ‚Äì Try It Online")** (Footer formats as a Python list for clarity.)
##### How?
Starts with `[[1]]` and repeatedly applies a function that creates the required new entries to the right (and then below, via the transposed matrix) by adding the number of current elements to each element.
```
ƒ†;"FL+∆äZ∆ä‚Å∏¬°‚Å∫ - Link: non-negative integer, n
ƒ†            - group indices (of implicit [n]) by value -> [[1]] (our initial M)
          ¬°  - repeat...
         ‚Å∏   - ...number of times: chain's left argument -> n
        ∆ä    - ...what: last three links as a monad, f(M):
      ∆ä      -   last three links as a monad, g(M):
   F         -     flatten
    L        -     length
     +       -     add to elements of M
  "          -   zip (rows of M) with (rows of that) and apply:
 ;           -     concatenation
       Z     -   transpose
           ‚Å∫ - repeat last link (i.e. repeat f() n more times)
```
[Answer]
# Python 3 with numpy, ~~91~~ ~~87~~ ~~74~~ ~~70~~ ~~69~~ 67 bytes
```
lambda n:c_[s:=(c_[:1<<n]&(b:=1<<r_[:n]))@b]*2+s
from numpy import*
```
[Attempt This Online!](https://ato.pxeger.com/run?1=m72soLIkIz9vweI025ilpSVpuhY3nXMSc5NSEhXyrJLjo4utbDWAlJWhjU1erJpGkpUtkFUEFMiL1dR0SIrVMtIu5korys9VyCvNLahUyMwtyC8q0YIYtbqgKDOvRCNNw1hTEyKyYAGEBgA)
`c_[:1<<n]` produces an array of the numbers 0 (inclusive) to 2n (exclusive), vertically. `r_[:n]` produces 0 (inc.) to n (exc.), horizontally, and the left shift makes the corresponding powers of 2, saved in `b`. The bitwise AND, with broadcasting, produces a table of values, splitting the numbers into their constituent bits. The matrix product (`@`) is equivalent to taking the dot product of each broken-down number with `b`; this squares each nonzero bit, changing powers of 2 into powers of 4, and adds them back up, producing a horizontal sequence 0, 1, 4, 5, 16, 17, 20, 21, ..., which is saved in `s`. `c_` makes a vertical version of that sequence, which is doubled and added to the original sequence, which (by broadcasting) produces a table of sums, which is the required result.
---
Previously:
```
from numpy import*
lambda n:(s:=vander([4],8)@unpackbits(mat(r_[:1<<n],uint8),0))+s.T*2
```
`r_[:1<<n]` produces 0 (inc.) to 2n (exc.). It is passed to `mat` ("Interpret the input as a matrix"), which in particular makes it 2-dimensional, while also changing its type to `uint8`, which is necessary for the next function `unpackbits`, which breaks the numbers down into their binary representations, across the added dimension. `vander([4],8)` produces an array of powers of 4, which in the `@` (matrix multiplication) gets dot-producted with each binary vector, producing a similar `s` (transposed differently), and the rest is similar.
[Answer]
# [Haskell](https://www.haskell.org/), 61 bytes
```
f n|r<-iterate(>>= \x->[4*x,4*x+1])[0]!!n=[[x+2*y|x<-r]|y<-r]
```
[Try it online!](https://tio.run/##y0gszk7Nyfn/P00hr6bIRjezJLUosSRVw87OViGmQtcu2kSrQgeItQ1jNaMNYhUV82yjoyu0jbQqaypsdItiaypB5P/cxMw829zEAl@FgqLMvBIFFYU0BeP/AA "Haskell ‚Äì Try It Online")
Same idea as tsh's Python, but the list `r` (e.g. `[0,1,4,5,16,17,20,21]`) is created in a Haskell-y way.
[Answer]
# [Julia 1.0](http://julialang.org/), 42 bytes
```
!x=(x-=1)<0||(~n=n*2^2x.+!x;[~0 ~1;~2 ~3])
```
[Try it online!](https://tio.run/##yyrNyUw0rPj/X7HCVqNC19ZQ08agpkajLs82T8sozqhCT1uxwjq6zkChztC6zkihzjhW839KZnFBTmKlhqKx5n8A "Julia 1.0 ‚Äì Try It Online")
(one-indexed)
[Answer]
# JavaScript (ES6), ~~¬†89 86¬†~~ 79 bytes
*Saved 7 bytes thanks to @tsh*
```
n=>[...Array(1<<n)].map((_,y,a)=>a.map(g=(k=1,x)=>k>>n?0:x&k|(y&k|g(k*2,x))*2))
```
[Try it online!](https://tio.run/##y0osSyxOLsosKNHNy09J/Z9m@z/P1i5aT0/PsagosVLD0MYmTzNWLzexQEMjXqdSJ1HT1i4RzE231ci2NdSpAApk29nl2RtYVahl12hUAol0jWwtI6CMppaRpub/5Py84vycVL2c/HSNNA0DTU0uVBFDDBEjDBFjoDkA "JavaScript (Node.js) ‚Äì Try It Online")
[Answer]
# [Haskell](https://www.haskell.org/), 88 bytes
```
f 0=[[1]]
f n=let(#)=zipWith(++);m=n-1;a=f m;[b,c,d]=map(map(+4^m))<$>[a,b,c]in a#b++c#d
```
[Try it online!](https://tio.run/##FYq7CgMhEAD7fMWCKRQN5Eg6s/cHV6cQA3sPOYkrkljdzxtTTDHM7PR9bym1FuCKzg3enwJkTFuVQuERyzPWXWqtLGO@DJYwAFs3m8WsHpmK/KPvL1bqcR4dmZ58zEBi1noRa2PqhtC3Ccon5goywE21Hw "Haskell ‚Äì Try It Online")
[Answer]
# [Ruby](https://www.ruby-lang.org/), ~~90...~~ 64 bytes
```
->n{(k=0...2**n).map{|y|k.map{|x|k.sum{|z|[x,y][z%2][z/2]<<z}}}}
```
[Try it online!](https://tio.run/##KypNqvyfZvtf1y6vWiPb1kBPT89ISytPUy83saC6prImG8KoADKKS3Ora6pqoit0KmOjq1SNgIS@UayNTVUtEPwvUEiLNo79DwA "Ruby ‚Äì Try It Online")
### Quick explanation:
Every number can be expressed by interleaving the binary digits of the number of row and column.
[Answer]
# [Pip](https://github.com/dloscutoff/pip), 19 bytes
```
$+{TBg}FB4*^21MC Ea
```
[Try it here!](https://replit.com/@dloscutoff/pip) Or, here's a 20-byte equivalent in Pip Classic: [Try it online!](https://tio.run/##K8gs@P9fRbs6xCm91s3JRCvOyNDX2UhLK/H///@6wf@NAQ "Pip ‚Äì Try It Online")
### Explanation
Inspired by [Neil's comment](https://codegolf.stackexchange.com/questions/231022/recursive-z-matrix#comment529704_231022) about interleaving bit patterns, we can simply map the x and y coordinates to the right value (example x = 5, y = 9):
```
5    101    0 1 0 1
9   1001   1 0 0 1    10010011   147
```
We're not going to use Pip's `WV` operator to do the interleaving, since it works left-to-right and doesn't fill in with zeros; instead, observe that an equivalent formula is to convert to base 4, multiply the y value by 2, and sum:
```
5    101      1 0 1      10001 (2) =  101 (4) =  17 (10)
9   1001   1 0 0 1    10000010 (2) = 2002 (4) = 130 (10)   147
```
So:
```
$+{TBg}FB4*^21MC Ea
                 Ea   2^input
              MC      Map this function to a grid of coordinate pairs of that size:
  {TBg}                Convert each coordinate to binary
       FB4             Treat as base 4 and convert back to decimal
          *            Multiply by
           ^21         list containing 2 and 1
                       (multiplying the row by 2 and the column by 1)
$+                     Sum
```
[Answer]
# [Python 3](https://docs.python.org/3/) + numpy, 66 bytes
```
from numpy import*
f=lambda m:c_[m and(f(m-.5),4**m//2+f(m-.5))].T
```
[Try it online!](https://tio.run/##K6gsycjPM/7/P60oP1chrzS3oFIhM7cgv6hEiyvNNicxNyklUSHXKjk@OlchMS9FI00jV1fPVFPHREsrV1/fSBvK14zVC/lfUJSZVwJUYaCpyQVjGyKxjZDYxpqa/wE "Python 3 ‚Äì Try It Online")
Mostly straight-forward recursion. Splits the 2x2 update into 2 1x2 updates with intermittent transpose.
# [Python 3](https://docs.python.org/3/) + numpy, 68 bytes
```
from numpy import*
f=lambda m:m and log2(kron(4**f(m-.5),c_[1:3])).T
```
[Try it online!](https://tio.run/##TcaxCoMwEADQ3a@48S6oqKmL0L/oJiJpJa3o3YUQB78@nQTf9MKZfio2Zx@VQQ4OJ6wcNCZT@Ofu@L044IHByQK7fjvcogo@jPHIVd1T@ZnHdrATUf3KIa6S0GNDVFxvb@9ut0T5Dw "Python 3 ‚Äì Try It Online")
*Just noticed that @Kirill L.'s R answer uses the exact same idea.*
By componentwise exponentiating adding in the original problem becomes multiplying so friend Kronecker can do the heavy lifting for us.
Downside: will get numerically inaccurate pretty fast.
[Answer]
# [BQN](https://mlochbaum.github.io/BQN/), 18 [bytes](https://github.com/mlochbaum/BQN/blob/master/commentary/sbcs.bqn)
```
{(‚çâ‚àæ‚çâ+¬∑‚â†‚•ä)‚çü2‚çüùï©‚âç‚âç0}
```
[Try it here.](https://mlochbaum.github.io/BQN/try.html#code=RuKGkHso4o2J4oi+4o2JK8K34omg4qWKKeKNnzLijZ/wnZWp4omN4omNMH0KCkbCqDDigL8x4oC/MuKAvzM=)
How it works:
```
{(‚çâ‚àæ‚çâ+¬∑‚â†‚•ä)‚çü2‚çüùï©‚âç‚âç0}
              ‚âç‚âç0  # Starting from a rank 2 array containing 0
            ‚çüùï©     # apply the preceding function ùï© times where ùï© is the input
  ‚çâ‚àæ               # The transposed array joined to
    ‚çâ+¬∑‚â†‚•ä          # the transposed array plus the number of elements
 (       )‚çü2       # applied twice
```
[Answer]
# [Octave](https://www.gnu.org/software/octave/)/MATLAB, 67 bytes
```
function x=f(n);x=1;for i=1:n;k=4^i/4;x=[x,k+x;2*k+x,3*k+x];end;end
```
[Try it online!](https://tio.run/##y08uSSxL/f8/rTQvuSQzP0@hwjZNI0/TusLW0Dotv0gh09bQKs8629YkLlPfBCgaXaGTrV1hbaQFJHWMQWSsdWpeCgj/T8ksLtBI0zDR1PwPAA "Octave ‚Äì Try It Online")
[Answer]
# [Python 3](https://docs.python.org/3/) + [numpy](https://numpy.org/doc/stable/index.html), 113 bytes
A very literal (and long) implementation of the spec.
```
lambda n:sum(tile(R(R(c_[:2,2:4].T,1<<i,0),1<<i,1)<<2*i,(1<<n+~i,)*2)for i in r_[:n])
from numpy import*
R=repeat
```
[Try it online!](https://tio.run/##JYuxDsIgFAD3fsUbeUiMpU4N/YnGTRuDCvEl5UGQDl38dSQxN1xuuLSXd@Sh@ulWVxseLws8frYgCq1OzI3n/Tpqpcfzcryo3hhSJ/y7R2O0JCVa8eFLCqVGHzMQEENuHy/Y@RwD8BbSDhRSzEV285RdcrbUlImL8GJArD8 "Python 3 ‚Äì Try It Online")
[Answer]
# JavaScript (ES6), 72 bytes
```
f=(n,a=[0])=>n?f(n-1,a.flatMap(x=>[x*=4,x+1])):a.map(y=>a.map(x=>x+2*y))
```
[Try it online!](https://tio.run/##JYlBCoMwEEX3nmKWGRODFle1k56gJxAXgzUlxSaiUiKlZ08DhQ/v8f6T37yNq1v2yof7lCzlCa@Y@npAMv5qha8axdrOvN94EZFMH0tqVZTNgHhm/cr1IPOXfEd5Kg/EZMMKwgFB3YGDC0GbKSXCpwAYg9/CPOk5PIQVDrErvukH "JavaScript (Node.js) ‚Äì Try It Online")
Use the idea from [Lynn](https://codegolf.stackexchange.com/users/3852/lynn)'s [Haskell answer](https://codegolf.stackexchange.com/a/231061/44718) which saves many bytes on generating the list for iterate. Using base conversion as what my Python answer do will [make it much longer (79 bytes)](https://tio.run/##HYpBDoIwEADvvGJPZLdoVcINt4nvMMY0CGQNtKQYQ6O@vSqnmUzmbp92boJMj63ztzZ1nBwbtHzWWp9CsBFLpRxd9Ggn7BmvGyE2kueSH96V6v9hVxLROkQ2dpWFzVKUKhKlzgdAAYZ9DQJHhurHoiB4ZQCNd7MfWj34HjsUojr7pC8).
[Answer]
# Python 3.8 + numpy, ~~85~~ 84 bytes
```
import numpy;z=lambda n:n and numpy.block([[m:=z(n-1),m+(k:=4**~-n)],[m+2*k,m+3*k]])
```
Not as clever as some of the other python ones, but shorter than the other python 3 ones. Also 0 returns a scalar instead of a list, not sure if that's allowed.
[Answer]
# [Japt](https://github.com/ETHproductions/japt), 14 [bytes](https://en.wikipedia.org/wiki/ISO/IEC_8859-1)
Port of [tsh's Python solution](https://codegolf.stackexchange.com/a/231033/58974).
```
2pU √á¬§n4
¬£√ã+X√ë
```
[Try it](https://petershaggynoble.github.io/Japt-Interpreter/?v=1.4.6&code=MnBVIMekbjQKo8srWNE&footer=VmMg%2bSDyMnBOKW24&input=MwotUg)
[Answer]
# [Pari/GP](http://pari.math.u-bordeaux.fr/), 62 bytes
```
n->matrix(2^n,,i,j,(g=k->fromdigits(binary(k-1),4))(i)*2+g(j))
```
[Try it online!](https://tio.run/##DccxDoAgDADArzi2WgbRFZ9igjE0hYCkMujr0duuehXDtYfB9WK27JvKA3YvREKRgF0yW9Arn8LSbjikeH0hmRlpRQTB0U4MEbFXldIqBFj@fA "Pari/GP ‚Äì Try It Online")
]