I don't think there's much hope for a general algorithm if the polynomials you are considering have degree so large that one can't store all the coefficients in memory. Then one can't present the polynomial by giving its coefficients, but has to do it some other way, and the right algorithm might depend a lot on how the polynomial is presented.

For reference, for a general polynomial of degree that may be large but not quite that large, I believe the usual algorithm is to compute the gcd of your polynomial with $x^p-x$, using Euclid's algorithm, noting that for the first step, of finding $x^p-x$ modulo $P$, it suffices to calculate $x^p \mod P$, which can be done rapidly by the exponentiation-by-squaring method.

In your case, the polynomial $R$ has some extra structure - its roots arise from chains of $2020$ consecutive quadratic residues. This algorithm seems less feasible than the direct algorithm of computing whether $x$ is a quadratic residue for each $x$ from $1$ to $p$ and halting if a chain is found, which "merely" takes $2^{2020}$ time and not $2^{2020}$ space.

I am skeptical that an algorithm much better than the direct one is available, but can't give a rigorous argument for this. (If the $2020$ appearing in the sum were smaller, one could express the number of solutions in terms of point counts on hyperelliptic curves of the form $y^2 = \prod_{i\in S} x+i$ for $S$ a subset of $\{1,2020\}$, and apply Kedlaya-type $p$-adic algorithms, but since there are $2^{2020}$ relevant hyperellipti curves this is not really any improvement.)