A Closed Formula for the Number of Solutions of 1/x + 1/y = 1/(x+y) mod n (OEIS A376202) and Fast Summation of Its Partial Sums

Let a(n) be the number of pairs 1 <= x <= y <= n-1 with gcd(x,n) = gcd(y,n) = gcd(x+y,n) = 1 and 1/x + 1/y = 1/(x+y) mod n, corresponding to the OEIS sequence A376202. The original entry contains an enumeration program of cost O(n^2) and several unproved observations. This paper shows that the condition is equivalent to the quadratic congruence x^2 + xy + y^2 = 0 mod n. We prove a closed formula: a(n) = 0.5 * phi(n) * rho(n) for n >= 2, n != 3, and a(3) = 2, where phi(n) is Euler's totient function and rho(n) is the number of roots of t^2 + t + 1 mod n. This formula completely determines a(n) from the prime factorization of n, proving the conjectures about zeros and cuban primes (A007645) listed in the OEIS entry. For the partial sums S(N) = sum_{n<=N} a(n), we establish an exact sublinear algorithm using a Min_25-type sieve that runs in O(N^(3/4) / log N) operations. We also prove an asymptotic formula S(N) = c * N^2 + O(N^(3/2+eps)) with an explicit constant c = 0.0792641301... In Python 3.12, computing the partial sum S(10^10) = 7,926,412,625,957,303,606 takes 19.7 s, and a single large term like a(10^15 + 3) takes just 0.57 s.