I am dealing with a polynomial of the form $$p(a,b) = a^n - b^n $$ for integer values $a > b$, and some small integer $n$. I am wanting to factor this polynomial for a large range of values (for example, let $a$ range from 2 to 1000000 and for each $a$ let $b$ range from 1 to $a-1$). At the moment I am factoring each value of $p(a,b)$ independently using a quadratic sieve in sage. It seems like this could be done much faster if I used some sort of sieve (similar to the sieve of Eratosthenes) to factor each value of the polynomial, knowing divisibility properties of different values would be intimately related. Ive tried to implement this but I cant seem to figure out how to do it.
Does anyone have any suggestions?
An obvious factor of $p(a,b)=a^n-b^n$ is of course $\gcd(a,b)^n$. So we might as well assume that $a$ and $b$ are coprime. The polynomial $p(a,b)$ itself factors into a product of cyclotomic polynomials $$p(a,b)=a^n-b^n=b^n\prod_{d\mid n}\Phi_d(\tfrac ab)=\prod_{d\mid n}b^{\varphi(d)}\Phi_d(\tfrac ab).$$ So to factor $p(a,b)=a^n-b^n$ you can first compute $b^{\varphi(d)}\Phi_d(\tfrac ab)$ for every divisor $d$ of $n$; this already yields a partial factorization of $p(a,b)$. The more small divisors $n$ has, the more small factors this produces.
EDIT: For example, if $n=3,5,7$ the polynomial factors as \begin{eqnarray*} a^3-b^3&=&(a-b)(a^2+ab+b^2),\\ a^5-b^5&=&(a-b)(a^4+a^3b+a^2b^2+ab^3+b^4),\\ a^7-b^7&=&(a-b)(a^6+a^5b+a^4b^2+a^3b^3+a^2b^4+ab^5+b^6) \end{eqnarray*} In general, if $n$ is prime this method only yields a factor $a-b$ and its cofactor.
Given that $n$ seems to be prime, it also makes sense to check primes $p\equiv1\pmod{n}$ first as these are heuristically $n$ times more likely to be a factor than other primes.