\\ Function a363763 for efficient calculation of single terms of A363763.
\\ Author: Hugo Pfoertner http://www.pfoertner.org/
\\ Change history:
\\ 2023-07-19 initial version
\\
\\ Auxiliary functions B3 and f3
\\ Expansion of Landau's function S(x) for the number of nonnegative integers <= x
\\ which can be expressed as a sum of two squares.
\\ Parameters fitted to 200000 terms of A077773
\\ instead of the Landau-Ramanujan and the Shanks constants
\\ (valid only asymptotically for x -> oo) are used,
\\ giving a better match of the observed distribution,
\\ also when compared with A164775 in a range x <10^16.
B3(x) = x/sqrt(log(x)) * (0.766747 + 0.30901/log(x) + 3.00549/(log(x))^2);
\\ expected number of hits in the interval [ n^2 + 1, n^2+2*n ]
f3(n) = B3((n+1/2)^2) - B3((n-1/2)^2) - 1;
\\ The limits of the scatter band of +-6*sqrt(n)
\\ are derived from the calculation of 200000 terms of A077773.
a363763(n) = {if (n<4,return([0,1,2,4][n+1]));
                  my (a77773mid = round (solve (x=n/2, 4*n, f3(x)-n)),
                  halfrange=ceil(6*sqrt(n)),
                  kfound=-1);
\\ diagnostic output
\\ print([n,a77773mid,halfrange]);
for (k = a77773mid-halfrange, a77773mid+halfrange,
     my (k1=k^2+1, k2=(k+1)^2-1, m);
\\ Consider using m=parsum(...) instead of m=sum(...) when multithreaded execution is available
     m = sum(j=k1,k2, sumdiv(j,d, (d%4==1)-(d%4==3))>0);
     if (m==n, kfound=k; break));
  kfound};
\\ Expected results: n=1000: 2570 (5s single core Intel I5 @ 3.3 GHz)
\\                   n=10000: 29356 (5 min real time 4-threaded on a Raspberry pi 4)
\\                   n=500000: 1747191 (~= 3 days single threaded on a Xeon W 2145 @ 4.3 GHz)