K (ngn/k), 12 bytes

-2/+2 2\0,4\

Attempt This Online! (ATO because the original ngn/k page isn't working at the time of posting)

Uses bit halving followed by base -2. 0, is necessary to get 0 0 for input 0.

-2/+2 2\0,4\
          4\    convert to base 4
        0,      prepend a 0
    2 2\        convert each to base 2, forcing two bits per a base-4 digit
                (base conversion happens column-wise, so the result here is
                 a two-row matrix)
   +            transpose
-2/             evaluate each column in base -2

K (ngn/k), 72 68 66 62 53 bytes

{(-|:)/[0<d;0,d:m*x-1+s*s-1]+(m:*/s#-1)*2#-2!s:-_-%x}

Try it online!

-4 : Squeezed out a few bytes
-2 : A couple more
-4 : Oh yeah, we’re golfing…
-9 : Forgot ngn/k has sqrt

JavaScript (Node.js), 48 bytes

f=(x,y)=>x&1?f(x/2,-~y):[x/2,y].map(t=>t/2^-t%2)

Attempt This Online!

Scala, 79 bytes

Port of @ბიმო's Haskell answer in Scala.

Golfed version. Try it online!

n=>for{i<-0 to n;x<-0 to i}yield Seq((x,i-x),(x,x-1-i),(-1-x,x-1-i),(-1-x,i-x))

Ungolfed version. Try it online!

object Main {
  def main(args: Array[String]): Unit = {
    val result = for {
      i <- 0 to 20
      x <- 0 to i
    } yield List((x, i - x), (x, x - 1 - i), (-1 - x, x - 1 - i), (-1 - x, i - x))

    result.flatten.zipWithIndex.take(21).foreach {
      case (tuple, index) => println(s"$index: $tuple")
    }
  }
}

JavaScript, 77 70 bytes

Z=x=>(x&1n?x:-x)/2n
f=(x,y=1n)=>++x&1n?[Z(++x/2n),Z(y)]:f(x/2n-1n,++y)

Takes a BigInt as Input.

Attempt This Online!

First few values:

f(0n)=[0n,0n]
f(1n)=[0n,-1n]
f(2n)=[-1n,0n]
f(3n)=[0n,1n]
f(4n)=[1n,0n]
f(5n)=[-1n,-1n]
...

Explanation:

Uses more or less the same algorithm as Wislaw's solution:

Z is the standard bijection from the positive integers to the integers.

f computes positive integers k and m with x+1=2^(k-1)*(2m-1) (which is a bijection \$\mathbb{N} \to \mathbb{N}_{>0} \times \mathbb{N}_{>0}\$ by the uniqueness of prime factorizations) and then outputs the pair Z(m),Z(k)

JavaScript, 65 bytes

bijection \$\mathbb{N}_{>0} \to \mathbb{Z} \times \mathbb{Z}\$

Z=x=>(x&1n?x:-x)/2n
f=(x,y=1n)=>x&1n?[Z(++x/2n),Z(y)]:f(x/2n,++y)

Attempt This Online!

Vyxal, 5 bytes

Þn:Ẋi

Try it Online!

Þn    # All integers
  :Ẋ  # Cartesian product with self
    i # Index input into that

J, 11 bytes

_1j1+.@#.#:

Attempt This Online!

Uses complex base -1 + i to generate a Gaussian integer from the binary representation of the given nonnegative integer.

_1j1+.@#.#:
         #:  Convert to binary
_1j1   #.    Evaluate the binary sequence as base -1 + i
    +.@      Take apart its real and imaginary parts

If outputting a complex number were allowed, it would be 8 bytes (delete +.@).

I couldn't find a complete proof that this conversion is bijective, but this challenge mentions it, and I think the following argument is enough to prove it:

Given a Gaussian integer, the lowest bit of base -1 + i representation is determined by (real+imag)%2. Subtract it and divide by -1 + i, and you always get another Gaussian integer with smaller magnitude (at least until it becomes small enough, at which point the proof can be done individually). This conversion algorithm is deterministic and finite, and defines an inverse of conversion from binary to Gaussian using base -1 + i.

05AB1E, 35 bytes

>©DÝʒo®sÖ}àsÅÉʒ®sÖ}à<2÷‚εDÈi2÷ë>2÷(

Try it online! or as Test suite

Explanation

Consider

$$\begin{align*} f: \mathbb{N} &\longrightarrow \mathbb{N} \times \mathbb{N} \\ n&\longmapsto (x,y)\, , \end{align*} $$ where \$x\$ is the largest number so that \$2^x\$ divides \$n+1\$, and where \$2y+1\$ is the largest odd number which divides \$n+1\$. The inverse of \$f\$ is the well-known bijection \$f^{-1}(x,y) = 2^x(2y+1)-1\$.

Then consider $$\begin{align*} g: \mathbb{N} \times \mathbb{N} &\longrightarrow \mathbb{Z} \times \mathbb{Z} \\ (m,n)&\longmapsto (h(m),h(n))\, , \end{align*} $$ where $$\begin{align*} h: \mathbb{N} &\longrightarrow \mathbb{Z} \\ n&\longmapsto \begin{cases}\frac n2, & n \text{ even}\\ -\frac{n+1}2, & n \text{ odd}\end{cases}\, . \end{align*} $$ Since \$f\$, \$g\$ and \$h\$ all are bijections, the composition \$g\circ f:\mathbb{N} \longrightarrow \mathbb{Z} \times \mathbb{Z}\$ is a bijection.

The program simply computes \$g\circ f\$.

>©DÝʒo®sÖ}àsÅÉʒ®sÖ}à<2÷‚εDÈi2÷ë>2÷( # Full program

                                    # Implicit input: Integer n
>©                                  # Compute n+1 and save it to the register
  DÝ                                # Duplicate n+1 and push the list [0,...,n+1]
    ʒo®sÖ}                          # Only keep those numbers x so that 2^x divides n+1
          à                         # Get maximum element in the list.
           sÅÉ                      # Swap so that n+1 is on top and push [1,3,5,...,n+1]
              ʒ®sÖ}                 # Only keep those numbers z which divides n+1
                   à<2÷             # Compute y = (z-1)/2
                       ‚            # Push the pair [x,y]
                        ε           # Apply the function h to x (and y):
                           i        # if...
                         DÈ         # x is even
                            2÷      # then compute x/2
                              ë>2÷( # else compute -(x+1)/2
                                    # Implicit output: [h(x),h(y)]

JavaScript, 166 168 bytes/characters

New approach using a rectangular spiral like others are using.

function f(n){return b=Math,k=b.ceil((b.sqrt(n)-1)/2),t=2*k+1,m=b.pow(t,2),t+=4,m-t>n?(m-=t,m-t>n?(m-=t,m-t>n?[k,k-(m-n-t)]:[-k+(m-n),k]):[-k,-k+(m-n)]):[k-(m-n),-k]}

I used this answer on Math.SE, translated it to JS and compressed it using UglifyJS.

This approach does not use any loops nor does it create the spiral in any way.

Because the spiral's coordinates cover all integers the function is bijective in the sense of \$f:\mathbb{N}^0\rightarrow \mathbb{Z}^2\$.

Update: Saved 2 characters by storing Math in b.

Update 2: Replaced t-=1 with t+=4 to fix the issue that caused \$f(0)=f(8)\$. This does not generate a spiral anymore but works for all non-negative integers \$\mathbb{N}^0\$ (all natural numbers including \$0\$).

Matlab, 54 bytes

n=input('')+1;[i,j]=find(spiral(2*n)==n);disp([i,j]-n)

The key here is spiral, this creates a spiral matrix of an arbitrary size.

spiral(3)

returns

ans =

 7     8     9
 6     1     2
 5     4     3

In Matlab, this works for integers, but you will get memory issues way before that, (well, I am assuming you as the command spiral creates a full matrix first, that has size of about \$4n^2\$. At around \$n \simeq 10^4\$ this matrix takes up more than 1GB of space. So with 1TB of RAM you will get to around \$n \simeq 10^5\$, and about \$2.9\cdot 10^{11}\$ GB of RAM will get you to \$n=2^{32}\$.

Pyth, 15

u,-HyeGhGjQ2,ZZ

Try it online.

u             reduce
                lambda G,H:    [implicit]
  ,-HyeGhG         (H-2*G[-1],G[0])
  jQ2           base(input(),2)
  ,ZZ           (0,0)
              print result     [implicit]

A Python translation:

g=lambda Z,n:(n-2*Z[1],Z[0])
print reduce(g,binlist(input()),(0,0))

or iteratively:

(x,y)=(0,0)
for b in binlist(input()):
    (x,y)=(b-2*y,x)
print (x,y)

where binlist converts a number to a list of digits like binlist(4) = [1,0,0].

So, how does this work? It interprets the binary digits of the number \$n\$ as two interleaved numbers in base negative two like in my Python solution.

The binary number $$n=\dots b_5 b_4 b_3 b_2 b_1 b_0 $$ corresponds to the pair $$ (x,y) = (b_0 - 2 b_2 + 4 b_4 - 8 b_6 + \cdots, b_1 - 2 b_3 + 4 b_5 - 8 b_7 + \cdots).$$

If we hadn't yet processed the last digit \$b_0\$ of \$n\$, we'd have all indices higher by $1$, $$n'=\dots b_5 b_4 b_3 b_2 b_1 $$ corresponding to the pair $$ (x',y') = (b_1 - 2 b_3 + 4 b_5 - 8 b_7 + \cdots, b_2 - 2 b_4 + 4 b_6 - 8 b_8 + \cdots).$$

We can then express the new values once \$b_0\$ is read in terms of the old values

$$(x,y) = (b_0 - 2y', x').$$

So, by repeatedly performing the transformation \$(x,y) \to (b - 2y, x)\$ on each successive bit \$b\$ of \$n\$, we build the point \$(x,y)\$ that comes from extracting its digits.

Haskell, 50 bytes

(0!).succ
l!n=(last$(!).succ:[(,)|odd n])l$div n 2

Try it online or try it with its inverse!

Ungolfed

ntoN2 n = 0 ! (n + 1)

xCounter ! remainingNum
  | odd remainingNum = (xCounter, div remainingNum 2)
  | otherwise        = (xCounter + 1) ! div remainingNum 2

Explanation

This uses the fact that each \$(x,y) \in \mathbb{N}^2\$ can be 1-to-1 mapped to \$2^x \cdot (2\cdot y + 1) - 1 \in \mathbb{N}\$. The above operator (!) computes \$x\$ by dividing the input as long as it's even, keeping track with the zero-initialized variable l (xCounter). Once we reached the even number an integer division computes y.

Note that the actual function f (ntoN2) increments the input before starting with the procedure.

Mathematica, 46

SortBy[Tuples[Range[2#]-#,2],Norm][[#]]&[#+1]&

Sort the vectors their norm, then take the nth one.

Python 2, 49

Edit: Improved to 49 by using better one-step recursion for base -2.

def f(n):x,y=n and f(n/2)or(0,0);return n%2-2*y,x

Here's a Pyth version using reduce.

Edit: Improved to 52 by switching to base -2 from balanced ternary.

Python 2, 52

h=lambda n:n and n%2-2*h(n/4)
lambda n:(h(n),h(n/2))

Python 2, 54

h=lambda n:n and-~n%3-1+3*h(n/9)
lambda n:(h(n),h(n/3))

This uses digit interleaving like Runer112's solution, but with balanced ternary rather than signed binary. Python doesn't have built-in base conversion, so the code implements it recursively.

The helper function h (with 3 in place of the 9) takes a natural number and converts it from ternary to balanced ternary with the digit substitutions

0 -> 0 
1 -> +1
2 -> -1

So, for example, 19, which is 201 in base, becomes (-1)(0)(+1) in balanced ternary, which is (-1)*3^2 +(0)*3^1+(+1)*3^0 = -8.

Balanced ternary suffices to encode every integer, and so gives a mapping from natural numbers to integers.

To map to pairs of integers, we interleave the digits in n. To do so, we have h look at every other digit by doing n/9 as the recursive step rather than n/3. Then, for one coordinate, we shift n by floor-dividing by 3.

Here are the first 81 outputs, which cover the region [-4,4]^2.

0 (0, 0)
1 (1, 0)
2 (-1, 0)
3 (0, 1)
4 (1, 1)
5 (-1, 1)
6 (0, -1)
7 (1, -1)
8 (-1, -1)
9 (3, 0)
10 (4, 0)
11 (2, 0)
12 (3, 1)
13 (4, 1)
14 (2, 1)
15 (3, -1)
16 (4, -1)
17 (2, -1)
18 (-3, 0)
19 (-2, 0)
20 (-4, 0)
21 (-3, 1)
22 (-2, 1)
23 (-4, 1)
24 (-3, -1)
25 (-2, -1)
26 (-4, -1)
27 (0, 3)
28 (1, 3)
29 (-1, 3)
30 (0, 4)
31 (1, 4)
32 (-1, 4)
33 (0, 2)
34 (1, 2)
35 (-1, 2)
36 (3, 3)
37 (4, 3)
38 (2, 3)
39 (3, 4)
40 (4, 4)
41 (2, 4)
42 (3, 2)
43 (4, 2)
44 (2, 2)
45 (-3, 3)
46 (-2, 3)
47 (-4, 3)
48 (-3, 4)
49 (-2, 4)
50 (-4, 4)
51 (-3, 2)
52 (-2, 2)
53 (-4, 2)
54 (0, -3)
55 (1, -3)
56 (-1, -3)
57 (0, -2)
58 (1, -2)
59 (-1, -2)
60 (0, -4)
61 (1, -4)
62 (-1, -4)
63 (3, -3)
64 (4, -3)
65 (2, -3)
66 (3, -2)
67 (4, -2)
68 (2, -2)
69 (3, -4)
70 (4, -4)
71 (2, -4)
72 (-3, -3)
73 (-2, -3)
74 (-4, -3)
75 (-3, -2)
76 (-2, -2)
77 (-4, -2)
78 (-3, -4)
79 (-2, -4)
80 (-4, -4)

An alternative coding with quarter-imaginary turned out longer, though it is very pretty.

Python 2, 63

h=lambda n:n and n%4+2j*h(n/4)
lambda n:(h(n).real,h(n).imag/2)

In a language with less clunky handling of complex conversion, this would probably be a better approach. If we could output complex numbers, we could do:

Python 2, 38

f=lambda n:n and n%2+n/2%2*1j-2*f(n/4)

Python 2, 98 bytes

Let's start out with a simple approach:

def f(N):
 x=a=0;b=2
 while N:x+=1j**b;b+=a<1;a=a or b/2;N-=1;a-=1
 return int(x.real),int(x.imag)

It just forms a rectangular spiral N units long on a 2d grid, starting from the origin, and returns the coordinates of the last point.

The function is bijective since:

• Each point can be covered, given a long enough spiral
• Each point will only be intersected by the spiral once

The spiral looks something like this (except starting at 0 rather than 1):

CJam, 24 22 21 bytes

My brain has trouble understanding the math other solutions are using. But my brain definitely understands binary, so here's a soultion based on bit manipulation!

li4b2fmd2/z{)(\2b^}%p

Try it online.

Explanation

This approach treats the input as two interleaved binary values, one for each output number. All but the least significant bit of each encode a magnitude, and the least significant bit signals whether or not to take the bitwise complement of this magnitude. In this implementation, the odd-positioned bits correspond to the first output number (and the even-positioned bits correspond to the second) and an LSB of 0 signals to take the complement.

For example, given an input of 73, uninterleaving its binary representation of 1001001b produces 0 1|0 (odd-positioned bits) and 1 0 0|1 (even-positioned bits). The first value has a magnitude of 01b = 1 and should be complemented for a final value of ~1 = -2, and the second value has a magnitude of 100b = 4 and should not be complemented.

Informal demonstration of correctness

I made a testing program that places each input from zero to a user-specified number minus one in its output location on a 2D grid. You can try it online as well. Here's an output of this program showing how the algorithm maps 0-99:

      -8 -7 -6 -5 -4 -3 -2 -1  0  1  2  3  4  5  6  7  8

-8                      92 84 86 94                     
-7                      88 80 82 90                     
-6                      76 68 70 78                     
-5                   96 72 64 66 74 98                  
-4                60 52 28 20 22 30 54 62               
-3                56 48 24 16 18 26 50 58               
-2                44 36 12  4  6 14 38 46               
-1                40 32  8  0  2 10 34 42               
 0                41 33  9  1  3 11 35 43               
 1                45 37 13  5  7 15 39 47               
 2                57 49 25 17 19 27 51 59               
 3                61 53 29 21 23 31 55 63               
 4                   97 73 65 67 75 99                  
 5                      77 69 71 79                     
 6                      89 81 83 91                     
 7                      93 85 87 95                     
 8                                                      

The fill pattern looks a bit weird, but it is in fact bijective! With each successive power of 4, it fills a square with double the previous side length. For example, here's how the algorithm maps 0-15:

      -2 -1  0  1  2

-2    12  4  6 14   
-1     8  0  2 10   
 0     9  1  3 11   
 1    13  5  7 15   
 2                  

This makes up the 4x4 square in the middle of the 8x8 square of 0-63:

      -4 -3 -2 -1  0  1  2  3  4

-4    60 52 28 20 22 30 54 62   
-3    56 48 24 16 18 26 50 58   
-2    44 36 12  4  6 14 38 46   
-1    40 32  8  0  2 10 34 42   
 0    41 33  9  1  3 11 35 43   
 1    45 37 13  5  7 15 39 47   
 2    57 49 25 17 19 27 51 59   
 3    61 53 29 21 23 31 55 63   
 4                              

Which makes up the 8x8 square in the middle of the 16x16 square of 0-255:

         -8  -7  -6  -5  -4  -3  -2  -1   0   1   2   3   4   5   6   7   8

 -8     252 244 220 212 124 116  92  84  86  94 118 126 214 222 246 254    
 -7     248 240 216 208 120 112  88  80  82  90 114 122 210 218 242 250    
 -6     236 228 204 196 108 100  76  68  70  78 102 110 198 206 230 238    
 -5     232 224 200 192 104  96  72  64  66  74  98 106 194 202 226 234    
 -4     188 180 156 148  60  52  28  20  22  30  54  62 150 158 182 190    
 -3     184 176 152 144  56  48  24  16  18  26  50  58 146 154 178 186    
 -2     172 164 140 132  44  36  12   4   6  14  38  46 134 142 166 174    
 -1     168 160 136 128  40  32   8   0   2  10  34  42 130 138 162 170    
  0     169 161 137 129  41  33   9   1   3  11  35  43 131 139 163 171    
  1     173 165 141 133  45  37  13   5   7  15  39  47 135 143 167 175    
  2     185 177 153 145  57  49  25  17  19  27  51  59 147 155 179 187    
  3     189 181 157 149  61  53  29  21  23  31  55  63 151 159 183 191    
  4     233 225 201 193 105  97  73  65  67  75  99 107 195 203 227 235    
  5     237 229 205 197 109 101  77  69  71  79 103 111 199 207 231 239    
  6     249 241 217 209 121 113  89  81  83  91 115 123 211 219 243 251    
  7     253 245 221 213 125 117  93  85  87  95 119 127 215 223 247 255    
  8                                                                        

Haskell, 78 74 bytes

(concat[[(x,i-x),(x,x-1-i),(-1-x,x-1-i),(-1-x,i-x)]|i<-[0..],x<-[0..i]]!!)

Test run:

*Main> mapM_ (print . (concat[[(x,i-x),(x,x-1-i),(-1-x,x-1-i),(-1-x,i-x)]|i<-[0..],x<-[0..i]]!!) ) [0..20]
(0,0)
(0,-1)
(-1,-1)
(-1,0)
(0,1)
(0,-2)
(-1,-2)
(-1,1)
(1,0)
(1,-1)
(-2,-1)
(-2,0)
(0,2)
(0,-3)
(-1,-3)
(-1,2)
(1,1)
(1,-2)
(-2,-2)
(-2,1)
(2,0)

How it works: list all pairs in the first quadrant in the following order

  |
 2| #4
  |
 1| #2  #5
  | 
 0| #1  #3  #6
  +---------------
     0   1   2   3 

mirror each point into the other quadrants to make a list of 4 element lists. Concatenate all into a single list and take the nth element.

Edit: function doesn't need a name, reordered math. expressions.

dc, 49

[1+2~2*1-*n]sm?dsa8*1+v1-2/dd1+*2/lar-dlmx32P-lmx

This starts off by arranging the non-negative integers on a grid thus:

..| 
4 | 14
3 |  9 13
2 |  5  8 12
1 |  2  4  7 11
0 |  0  1  3  6 10
Y +-----------------
  X  0  1  2  3  4 ...

Note that how the grid positions are filled diagonally with increasing N. Notice the Y=0 line contains the triangular number sequence, given by N = X(X+1)/2. This is a quadratic equation which is solved using the normal formula, using just the +ve root, so that we can determine X from N when Y=0. Next is some simple arithmetic shuffling to give unique {X,Y} for every N.

This provides the required bijective quality, but X and Y are only non-negative, but the question requires all possible integers. So X and Y are mapped using ((t+1)/2)*((t+1)~2*2-1) to give all possible integers.

dc has arbitrary precision numbers, so input range to 2^31-1 is no problem. Note also that default precision is 0 decimal digits, and the sqrt() and / round down which is the behavior needed here.

Output:

$ for i in {0..10}; do dc biject.dc <<< $i; echo; done
0 0
0 -1
-1 0
0 1
-1 -1
1 0
0 -2
-1 1
1 -1
-2 0
0 2
$