According to Theorem 2 of K. Ono, S. Robins, and P. T. Wahl, On the representation of integers as sums of triangular numbers, Aequationes Math. 50 (1995), 73–94., define with $\left( \frac{r}{n} \right)$ being the Jacobi symbol, $$R_3(n) = \begin{cases} 24 \sum_{r=1}^{\lfloor n/4 \rfloor} \left( \frac{r}{n} \right), & \text{if}\ n \equiv 1 \mod 4 \\ 8 \sum_{r=1}^{\lfloor n/2 \rfloor} \left( \frac{r}{n} \right), & \text{if}\ n \equiv 3 \mod 4 \\ \end{cases} $$
The number of ways to represent $n$ as the sum of $3$ squares is $$r_3(n) = \sum_{d^2 \mid n} R_3 \left( \frac{n}{d^2} \right)$$
My question is how to efficiently compute the sum of Jacobi symbols for a fixed denominator and range of numerators, if possible.
It's well known the Jacobi symbol is completely multiplicative: $$\left( \frac {ab} n \right) = \left( \frac a n \right) \left( \frac b n \right)$$ So perhaps it's possible to use a sieve to calculate all values in range based on their prime factors, though the time complexity is close to linear. A sublinear algorithm would be ideal (I'm fairly certain I can take the algorithmic DP approach for the original problem).
