\\ 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)