Initialize an array \(C[n]\) with indices up to and including \(N\), then perform the following algorithm to directly count for each \(n\) the number of pairs \((x, y)\) that satisfies \(0 < x < y\) and \(x^2 + y^2 = n\):

\(\texttt{for } x \in \mathbb Z \cap [1, \infty):\)

\(\qquad \texttt{if }x^2 > N: \texttt{break}\)

\(\qquad \texttt{for } y \in \mathbb Z \cap [x+1, \infty):\)

\(\qquad \qquad n \leftarrow x^2 + y^2\)

\(\qquad \qquad \texttt{if }n > N: \texttt{break}\)

\(\qquad \qquad C[n] \xleftarrow {+} 1\)

\(\qquad \texttt{end for}\)

\(\texttt{end for}\)

Then collect in an array-list the set of \(n\) where \(C[n] = 1\), then print the result.

As the number of positive integers \(x\) such that \(x^2 \le N\) is \(O(\sqrt N)\), the amount of time this would take is \(O(\sqrt N)^2 = O(N)\)