For which values of i is (y_i² − n) divisible by p_j? Define a = floor(sqrt(n)), a constant.
(y_i² ≠ n) = (a+i)² − n. Given p_j, we need to know for which values of i does p_j divide (a+i)² − n. For half of the primes there will be no solutions and for the other half there will be two solutions per period of size p_j. Solving this equation involves taking square root modulo some small prime, not a particularly difficult thing to do. If i solves the equation then so does i+p_j. A simple loop thus augments the counters for all those i’s where p_j divides (y_i² − n).