Probability on Number Theory Problem: Suppose that $a,b,c \in \{1,2,3,\cdots,1000\}$ are randomly selected with replacement. Find the probability that $abc+ab+2a$ is divisible by $5$.
Answer given from the worksheet: $33/125$

My answer: $\frac{641}{3125}$
Attempt: Since $abc+ab+2a = a(bc+b+2)$, either $a \equiv 0 \pmod{5}$ or $bc+b+2 \equiv 0 \pmod{5}$. The first case, which is $a \equiv 0 \pmod{5}$, has a probability of $1/5$. Now on the other case,
$$bc+b+2 \equiv 0 \pmod{5} \Rightarrow b(c+1) \equiv 3 \pmod{5}$$
happens when $b\equiv 1$ and $c \equiv 2$, $b \equiv 2$ and $c \equiv 3$, $b \equiv 3$ and $c \equiv 0$, or $b \equiv 4$ and $c \equiv 1$, total of $4$ solutions. So, this case has a probability of $4 \cdot\frac{1}{5^4} = \frac{4}{625}$. By Inclusion-exclusion principle, the final probability should be $\frac{1}{5} + \frac{4}{625} - \frac{1}{5} \cdot \frac{4}{625} = \frac{641}{3125}.$
This is apparently not the same from the given answer in the worksheet. Where did I go wrong?

Brute forcing this in Python, I get the same answer of $41/125 = 0.328$ previously given more thoroughly by John Omielan, which differs from the worksheet's answer of $33/125$:
import itertools as it
N = 100
combos = list(it.product(range(1,N+1), repeat=3))
exprs = [a*b*c + a*b + 2*a for (a,b,c) in combos]
print(len([x for x in exprs if x%5==0]) / len(exprs))
>>> 0.328

(note: this code only selects $a, b, c \in \{ 1, 2, 3, ..., 100 \}$ for memory reasons, but this should not matter, as long as the upper limit of the set of integers is any multiple of 5.)
