I want to factor $(n-1)^3+1$ for many values of $n$. Meaning, for every value of $n$ from $1$ to $N$ for some large $N$ I want a list of prime factors and corresponding powers for $(n-1)^3+1$. We can see that $n$ divides out and is linear so its contribution is easy to account for with sieving, but then we are left with $n^2-3n+3$.
Obviously this does not factor over the integers, but its outputs can easily factor. For primes larger than 3 we have: $$\begin{align*} p|n^2-3n+3\Rightarrow n^2-3n+3&\equiv 0\mod p\\ 4n^2-12n+12&\equiv 0\mod p\\ (2n-3)^2&\equiv-3\mod p\\ n&\equiv 2^{-1}(3\pm\sqrt{-3})\mod p \end{align*}$$ (to abuse notation a little with the square root). Thus for $p$ to possibly divide $n^2-3n+3$ we need $-3$ to be a quadratic residue, which means $p\equiv 1\mod 6$.
For 2 and 3, 2 can never divide the output because $n^2-3n+3$ is always odd, and for 3 to divide the output 3 must divide $n^2$ so it must divide $n$.
But I'm a little bit stuck here: continuing to set aside 2 and 3, for larger primes we start by finding a square root of $-3$ say $r$. Then we add $p$ to the list of prime factors of the output of the polynomial for $n$ for all values of $n$ of the form $\frac{p+1}{2}(3\pm r)+kp$. I think finding a huge number of square roots will be very expensive using Cipolla's algorithm or Tonelli's algorithm though, so I'm wondering if there is a more efficient way I'm not thinking of. I'm also unsure how to handle finding not just the set of primes that divides $n^2-3n+3$ but the powers of these primes as well. However, working through the above computations again I think Hensel lifting the square root $r$ of $-3$ to $p^k$ and using the same approach would work for higher powers of $p$. I'm not sure if there's a more straightforward way though.