
MATHEMATICA

fx[n_] := fx[n] = If[n == 1, 0, fx[n  1] + Sin[#*Pi/2]& @ Mod[Floor[Sqrt[ 4*(n  2) + 1]], 4]];
fy[n_] := fy[n] = If[n == 1, 0, fy[n  1]  Cos[k*Pi/2]& @ Mod[Floor[Sqrt[ 4*(n  2) + 1]], 4]];
b[_, _] = 0;
a[n_] := Module[{x, y, s, i, t, m}, {x, y} = {fx[n + 1], fy[n + 1]}; If[b[x, y] > 0, b[x, y], s = {};
For[i=1, True, i++, t = b[x+i, y+i]; If[t>0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[xi, yi]; If[t>0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[x+i, yi]; If[t>0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[xi, y+i]; If[t>0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[x+i, y]; If[t > 0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[xi, y]; If[t > 0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[x, y+i]; If[t > 0, s = Union[s, {t}], Break[]]];
For[i=1, True, i++, t = b[x, yi]; If[t > 0, s = Union[s, {t}], Break[]]];
m = 1; While[MemberQ[s, m], m++]; b[x, y] = m]];
Flatten[Position[a /@ Range[0, 10^4], 1]]  1 (* JeanFrançois Alcover, Feb 25 2020, after Alois P. Heinz *)
