Extreme Polynomial Division problem

Note that $x^4 + x^3 + 2x^2 + x + 1 = (x^2+1)(x^2+x+1).$

Then modulo $x^2+1$ , $x^{1001} - 1 \equiv x(x^2)^{500} - 1 \equiv x(-1)^{500} - 1 \equiv x-1 \pmod{x^2+1},$ and modulo $x^3-1 = (x-1)(x^2+x+1)$ , $x^{1001} - 1 \equiv x^2(x^3)^{333} - 1 \equiv x^2(1)^{333} - 1 \equiv x^2 - 1 \pmod{x^3-1},$ so modulo $x^2+x+1$ , $x^{1001} - 1 \equiv x^2 - 1 \equiv -x - 2 \pmod{x^2+x+1}.$

---

We can conclude that if $R(x)$ is the remainder $x^{1001} - 1 \equiv R(x) \pmod{x^4+x^3+2x^2+x+1},$ then there is a polynomial $P(x)$ such that $R(x) = -x - 2 + P(x)(x^2 + x + 1) \equiv x-1 \pmod{x^2+1} \\ xP(x) \equiv 2x+1 \equiv 2x+1 - (x^2+1) \equiv 2x-x^2 \pmod{x^2+1} \\ P(x) \equiv 2-x \pmod{x^2 + 1}$

Therefore $R(x) = -x-2 + (2-x)(x^2+x+1) = x^2 - x^3$

We can factor the divisor as

$\begin{array}{rl} x^4 + x^3 + 2x^2 + x + 1 &= \ \ (x^4 + 2x^2 + 1) + (x^3 + x) \\ &= \ \ (x^2 + 1)^2 + x(x^2 + 1) \\ &= \ \ (x^2 + 1)(x^2 + x + 1) \end{array}$

We know the roots (zeroes) of $x^2 + 1$ are $i$ and $\text{-}i$ ; also, the roots of $x^2 + x + 1$ are $\omega$ and $\omega^2$ , the complex cube roots of unity. [ See Note 1 below. ] We also know that since our divisor is of degree $4$ , the remainder $R(x)$ must be of degree $3$ or lower. Then we can write

$\begin{array}{rl} P(x) = x^{1001} - 1 &= \ \ Q(x) \cdot (x^4 + x^3 + 2x^2 + x + 1) + R(x) \\ &= \ \ Q(x) \cdot (x - i)(x + i)(x - \omega)(x + \omega) + Ax^3 + Bx^2 + Cx + D \end{array}$

Now we evaluate this at two of the four roots mentioned above.

$\begin{array}{rllll} P(i) = i^{1001} - 1 &= \ \ Q(i) \cdot 0 + Ai^3 + Bi^2 + Ci + D \\ \implies (i^4)^{250} \cdot i - 1 &= \ \ Ai^3 + Bi^2 + Ci + D \\ \implies i - 1 &= \ \ \text{-}Ai - B + Ci + D & & & \small \text{[ Since } i^4 = 1 \ \& \ i^2 = \text{-}1 \ ] \\ \\ \implies C - A &= \ \ 1 \quad \& \quad B - D \ \ = \ \ 1 \\ \\ P(\omega) = \omega^{1001} - 1 &= \ \ Q(\omega) \cdot 0 + A\omega^3 + B\omega^2 + C\omega + D \\ \implies (\omega^3)^{333} \cdot \omega^2 - 1 &= \ \ A\omega^3 + B\omega^2 + C\omega + D \\ \implies \omega^2 - 1 &= \ \ A + B \omega^2 + C \omega + D & & & \small \text{[ Since } \omega^3 = 1; \text{ see Note 1 below ]} \\ \\ \implies B \ \ = \ \ 1 \quad \& \quad C &= \ \ 0 \quad \& \quad A + D \ \ = \ \ \text{-}1& & & \small \text{[ See Note 3 below ]} \end{array}$

Combining these results, we get $A = \text{-}1, B = 1, C = D = 0$ . Thus our remainder is $R(x) = Ax^3 + Bx^2 + Cx + D = \boxed{\text{-}x^3 + x^2}$

Note 1: The equation $x^3 = 1$ can be re-written as $(x - 1)(x^2 + x + 1) = 0$ ; the second factor yields the complex roots $\frac{-1 \pm \sqrt{3}i}{2}$ . These complex cube roots of unity are usually called $\omega$ and $\omega^2$ as it's easily verified that each is the square of the other. This also, of course, means that $\omega^3 = 1$ .

Note 2: We could have also evaluated at - $i$ and $\omega^2$ , but we would have gotten the exact same equations for the coefficients as with $i$ and $\omega$ .

Note 3: We can equate coefficients as we did because it is easy to verify that $B \omega^2 + C \omega = \omega^2$ can only be true if $B = 1$ and $C = 0$ .

Note that $x^4 + x^3 + 2x^2 + x + 1$ divides evenly into $x^{12} - 1$ , which means in a long division setup of $x^{1001} - 1$ divided by $x^4 + x^3 + 2x^2 + x + 1$ , the same steps would be cycled through every $12$ terms, and so its remainder would be the same as the remainder of $x^{1001 - 12n} - 1$ divided by $x^4 + x^3 + 2x^2 + x + 1$ for any integer $n$ .

Letting $n = 83$ simplifies the problem to finding the remainder when $x^5 - 1$ is divided by $x^4 + x^3 + 2x^2 + x + 1$ , which is $\boxed{-x^3 + x^2}$ .

Let's begin by factoring the divisor $x^4 + x^3 + 2x^2 + x + 1 = \left( x - i \right) \left( x + i \right) \left( x - \omega \right) ( x - \omega^2 )$ where $\omega = \dfrac{ -1 + i \sqrt{3} }{ 2 }$

Now by simple substitution and the Remainder Factor Theorem $\begin{aligned} x^{1001} - 1 &\equiv i - 1 &\pmod{ x - i } \\ x^{1001} - 1 &\equiv -i - 1 &\pmod{ x + i } \\ x^{1001} - 1 &\equiv \omega^2 - 1 &\pmod{ x - \omega } \\ x^{1001} - 1 &\equiv \omega - 1 &\pmod{ x - \omega^2 } \\ \end{aligned}$

Now this can just be solved like one would solve system of linear congruences in number theory (with help of Remainder Factor Theorem ) $\begin{aligned} x^{1001} - 1 &\equiv \left( i - 1 \right) \dfrac{ \left( x + i \right) \left( x - \omega \right) \left( x - \omega^2 \right) }{ \left( i + i \right) \left( i - \omega \right) \left( i - \omega^2 \right) }&\pmod{ x - i } \\\\ x^{1001} - 1 &\equiv \left( -i - 1 \right) \dfrac{ \left( x - i \right) \left( x - \omega \right) \left( x - \omega^2 \right) }{ \left( -i - i \right) \left( -i - \omega \right) \left( -i - \omega^2 \right) } &\pmod{ x + i } \\\\ x^{1001} - 1 &\equiv \left( \omega^2 - 1 \right) \dfrac{ \left( x - i \right) \left( x + i \right) \left( x - \omega^2 \right) }{ \left( \omega - i \right) \left( \omega + i \right) \left( \omega - \omega^2 \right) } &\pmod{ x - \omega } \\\\ x^{1001} - 1 &\equiv \left( \omega - 1 \right) \dfrac{ \left( x - i \right) \left( x + i \right) \left( x - \omega \right) }{ \left( \omega^2 - i \right) \left( \omega^2 + i \right) \left( \omega^2 - \omega \right) } &\pmod{ x - \omega^2 } \\\\ \end{aligned}$

Finally $\begin{aligned} x^{1001} - 1 &\equiv \left( i - 1 \right) \dfrac{ \left( x + i \right) \left( x - \omega \right) \left( x + \omega \right) }{ \left( i + i \right) \left( i - \omega \right) \left( i + \omega \right) } \\\\ &+ \left( -i - 1 \right) \dfrac{ \left( x - i \right) \left( x - \omega \right) \left( x + \omega \right) }{ \left( -i - i \right) \left( -i - \omega \right) \left( -i + \omega \right) } \\\\ &+ \left( \omega^2 - 1 \right) \dfrac{ \left( x - i \right) \left( x + i \right) \left( x + \omega \right) }{ \left( \omega - i \right) \left( \omega + i \right) \left( \omega + \omega \right) } \\\\ &+ \left( \omega - 1 \right) \dfrac{ \left( x - i \right) \left( x + i \right) \left( x - \omega \right) }{ \left( \omega^2 - i \right) \left( \omega^2 + i \right) \left( \omega^2 - \omega \right) } \\\\ &\equiv x^2 - x^3 \pmod{ x^4 + x^3 + 2x^2 + x + 1 } \end{aligned}$

In general, we can find $P(z) \mod D(z)$ for any $D \neq 0$ if $D$ can be factored, by using this method,
We just need to construct the expression in the last step and evaluate it.
(Note the similarity with the Lagrange Interpolation )

By the way, WolframAlpha can do this for you.