A Closed Formula and Fast Summation for the Number of Pairs (x, y) with x^2 + y^2 Dividing n (OEIS A239680)

Let a(n) be the number of pairs (x,y) in {0, ..., n-1}^2 such that x^2+y^2 divides n (OEIS A239680). The original OEIS entry contains only a brute-force program of cost O(n^2). We prove an exact multiplicative formula a(n) = h(n) + s(n) - 2*[n=1], where h = 1*1*chi_4 is a Dirichlet convolution with h(2^e) = e+1, h(p^e) = binomial(e+2, 2) for p = 1 mod 4, h(p^e) = floor(e/2)+1 for p = 3 mod 4, and s(n) is the number of square divisors of n. Equivalently, the Dirichlet generating series of a(n) is established as zeta(z)^2 * L(z, chi_4) + zeta(z)*zeta(2z) - 2. For the partial sums S(N) = sum_{n<=N} a(n), we prove the exact reduction to the Gauss circle problem: S(N) = sum_{k=1}^N G(floor(N/k)) + sum_{m>=1} floor(N/m^2) - 2, where G(M) counts the lattice points inside a closed quadrant disc x^2+y^2<=M. By replacing the traditional direct evaluation of G(M) with a sublinear walk along the convex hull of the lattice points of the disc (a Stern-Brocot descent with an exact quadratic step), the algorithm evaluates G(M) in approximately 1.13 * M^(1/3) iterations, verified numerically up to M=10^18. In Python 3.12, the sum S(10^11) = 2,185,202,802,637 is computed in 1.5 s, G(10^18) in 3.1 s, and a single term a(10^15+7) takes 0.06 s.