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?
 A: You've made a good attempt. The only mistake I can determine is that, modulo $5$, there are only $5$ possibilities for each of $b$ and $c$, so there are $5^2 = 25$ total potential combinations, not $5^4$. Thus, your $4 \cdot\frac{1}{5^4} = \frac{4}{625}$ should be $4 \cdot\frac{1}{5^2} = \frac{4}{25}$ instead. The final probability using your methodology (i.e., the inclusion-exclusion principle) would then be
$$\frac{1}{5} + \frac{4}{25} - \frac{1}{5} \cdot \frac{4}{25} = \frac{41}{125} \tag{1}\label{eq1A}$$
However, this is also different from what the worksheet states. It seems the problem statement has an issue (e.g., part of it is incorrect, it's missing something, etc.), the worksheet answer is incorrect, or we're both missing something. I assume no explanation is provided regarding how the worksheet's answer was determined.
Also, note that regarding how to calculate this, I would instead consider the $2$ events of $a \equiv 0 \pmod{5}$ and $a \not\equiv 0 \pmod{5}$. The first one has a $\frac{1}{5}$ chance, with $b$ and $c$ being anything so their probability is $1$, giving an overall probability of $\frac{1}{5}$, as you determined. Otherwise, there's a $\frac{4}{5}$ chance that $a \not\equiv 0 \pmod{5}$, with a $\frac{4}{25}$ probability that $b$ and $c$ will be appropriate values. This then gives the probability to be
$$\frac{1}{5} + \frac{4}{5} \cdot \frac{4}{25} = \frac{41}{125} \tag{2}\label{eq2A}$$
i.e., the same as in \eqref{eq1A}. That's because $\frac{4}{25} - \frac{1}{5} \cdot \frac{4}{25} = \frac{5}{5} \cdot \frac{4}{25} - \frac{1}{5} \cdot \frac{4}{25} = \frac{5-1}{5}\cdot \frac{4}{25} = \frac{4}{5} \cdot \frac{4}{25}$.
A: 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.)
