Return to Answer

Here's a simple solution using Nearest. I've written it assuming the sites are in a list ps, rather than embedded inside a function p; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

Here's a simple solution using a single precomputed NearestFunction; much faster than $O(n^2)$. I've written it assuming the sites are in a list ps, rather than embedded inside a function p, because I think this way is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]] (* so it returns the index of the site *)

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

Here's a simple solution using Nearest. I've written it assuming the sites are in a list ps, rather than embedded inside a function p; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

P.S.

1. You could just define nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]], but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the NearestFunction.

2. I'm not too happy about having to put the N around EuclideanDistance. Unfortunately, NearestFunction doesn't accept something like $\sqrt2$ as the search radius.

Here's a simple solution using Nearest. I've written it assuming the sites are in a list ps, rather than embedded inside a function p; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

Edit: I forgot that you want the neighbours of all the sites collected in a big list. Well, then you just do

nearestSites = nthNearestSites[#, 1] & /@ Range[Length[ps]];
nextNearestSites = nthNearestSites[#, 2] & /@ Range[Length[ps]];

On a hundred thousand sites, these take 2.9 and 6.8 seconds on my machine respectively. On a million, they will probably take a couple of minutes.

P.S.

1. You could just define nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]], but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the NearestFunction.

2. I'm not too happy about having to put the N around EuclideanDistance. Unfortunately, NearestFunction doesn't accept something like $\sqrt2$ as the search radius.

Here's a simple solution using Nearest. I've written it assuming the sites are in a list ps, rather than embedded inside a function p; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes just a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

Here's a simple solution using Nearest. I've written it assuming the sites are in a list ps, rather than embedded inside a function p; I think this is easier to generate and manipulate. You may want to modify the code as appropriate.

num = 5000;
ps = RandomInteger[{1, 100}, {num, 2}];

nf = Nearest[Table[ps[[i]] -> i, {i, Length[ps]}]]

upToNthNearestSites[k_, 0] := {k} (* the "zeroth" nearest neighbour, i.e. itself *)
upToNthNearestSites[k_, n_] := Module[{pk, near, d},
  pk = ps[[k]];
  near = upToNthNearestSites[k, n - 1]; (* get the nearest neighbours up to order n-1 *)
  near = nf[pk, Length[near] + 1]; (* get one more; this is one of the nth nearest *)
  d = N@EuclideanDistance[pk, ps[[Last@near]]]; (* distance to the nth nearest *)
  nf[pk, {Infinity, d}] (* the solution is all the sites up to that distance *)
  ]
nthNearestSites[k_, n_] := Module[{pk, near0, near, d},
  pk = ps[[k]];
  near0 = upToNthNearestSites[k, n - 1];
  near = nf[pk, Length[near0] + 1];
  d = N@EuclideanDistance[pk, ps[[Last@near]]];
  near = nf[pk, {Infinity, d}];
  Complement[near, near0] (* same as above except remove neighbours closer than n *)
  ]

The nearest neighbours to the $k$th site are given by nthNearestSites[k, 1], the second nearest by nthNearestSites[k, 2], and so on. On my machine, even with a million points, the initial construction of nf takes a little over a second, and after that nthNearestSites[1, 2] takes negligible time.

P.S.

1. You could just define nthNearestSites[k_, n_] := Complement[upToNthNearestSites[k, n], upToNthNearestSites[k, n - 1]], but that would end up evaluating the $(n-1)$th neighbours twice (as well as the $(n-2)$th, the $(n-3)$th, and so on). In the implementation above, it makes exactly $2n$ calls to the NearestFunction.

2. I'm not too happy about having to put the N around EuclideanDistance. Unfortunately, NearestFunction doesn't accept something like $\sqrt2$ as the search radius.