My approach is as follows: for one quadrant I find number of points with distance d > n — 1 and d <= n that does not lie on the axis. Then I multiply that by 4 (because we have 4 quadrants) and add 4 (because we have 4 points on the axes in total). That seems to be incorrect for n >= 5, I can't get why (for n >= 5 I don't have an integer overflow so the algorithm is incorrect then).
UPD: pranavsindura pointed out why my solution is incorrect. Looking for a correct one.