Spiral neighbourhoods If we take the natural numbers and roll them up counter clock-wise into a spiral we end up with the following infinite spiral:

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


Given some number in that spiral your task is to determine its neighbours - meaning the element above, left, right and below it.

Example

If we have a look at 27 we can see that it has the following neighbours:


above: 28
left: 10
right: 52
below: 26


So the output would be: [28,10,52,26]

Rules


Input will be a number \$n \geq 0\$ in any default I/O format
Output will be a list/matrix/.. of that numbers' 4 neighbours in any (consistent!) order
You may work with a spiral that starts with 1 instead of 0, however you should specify that in your answer


Examples

The output is in the format [above,left,right,below] and uses a 0-based spiral:

0  -&gt;  [3,5,1,7]
1  -&gt;  [2,0,10,8]
2  -&gt;  [13,3,11,1]
3  -&gt;  [14,4,2,0]
6  -&gt;  [5,19,7,21]
16  -&gt;  [35,37,15,17]
25  -&gt;  [26,24,50,48]
27  -&gt;  [28,10,52,26]
73  -&gt;  [42,72,74,112]
101  -&gt;  [100,146,64,102]
2000  -&gt;  [1825,1999,2001,2183]
1000000  -&gt;  [1004003,1004005,999999,1000001]
R, 156 bytes



function(n){g=function(h)c(0,cumsum(h((4*(0:(n+2)^2)+1)^.5%%4%/%1/2)))
x=g(sinpi)
y=g(cospi)
a=x[n]
b=y[n]
which(x==a&amp;(y==b+1|y==b-1)|y==b&amp;(x==a+1|x==a-1))}


Try it online!


posted another R answer since it's a slightly different approach than @ngn
1-indexed
neighbours are always sorted by ascending value
saved 6 bytes removing round and using cospi(x)/sinpi(x) which are more precise than cos(x*pi)/sin(x*pi) in case of half numbers (0.5, 1.5 etc...)
saved another byte removing the minus on y coordinates since the result is the same (just up/down neighbours are reversed)




Explanation :

If we look at the matrix coordinates of the values, considering the first value 0 placed at x=0, y=0, they are :

x = [0,  1,  1,  0, -1, -1, -1,  0,  1,  2,  2,  2,  2,  1,  0, ...] 
y = [0,  0,  1,  1,  1,  0, -1, -1, -1, -1,  0,  1,  2,  2,  2, ...]


The x coordinates follow the A174344 OEIS sequence with the recursive formula :  

a(1) = 0, a(n) = a(n-1) + sin(mod(floor(sqrt(4*(n-2)+1)),4)*pi/2)


The same formula holds for y matrix coordinates, but with cos instead of sin and negated :

a(1) = 0, a(n) = a(n-1) - cos(mod(floor(sqrt(4*(n-2)+1)),4)*pi/2)


So, in R we can translate the formula to this function, taking sinpi/cospi as parameter :

g=function(h)c(0,cumsum(h((4*(0:(n+2)^2)+1)^.5%%4%/%1/2)))


and we generate the two coordinates vectors (we don't negate the y coords since we'll get the same result, just with up/down neighbours reversed) :

x=g(sinpi)
y=g(cospi)


Note that we have generated (n+2)^2 coordinates, which are more than the minimum necessary coordinates containing both n and their neighbours (a tighter bound would be (floor(sqrt(n))+2)^2 but unfortunately is less "golfy").

Therefore, now that we have all the coordinates, we first search the coordinates a,b corresponding to our n : 

a=x[n]
b=y[n]


finally we select the positions of their neighbours, i.e. : 


the up/down neighbours where x == a and y == b+1 or b-1
the right/left neighbours where y == b and x == a+1 or a-1


using :

which(x==a&amp;(y==b+1|y==b-1)|y==b&amp;(x==a+1|x==a-1))