I know I have made a mistake! But I can't figure out where, and I'm just hoping someone can explain where/why.
The Legendre Sieve, if we use the floor function, gives the formula
$$\pi_L(x) = \pi(x) - \pi(\sqrt{x}) + 1 = \sum_{d \mid P} \mu(d)\left[\frac{x}{d}\right] \tag{1}$$
Where $d \mid P$ indicates the set of divisors (including $1$) of $P = p_n\text{#} = p_1p_2 \cdots p_n$, with $p_1 = 2$ and $p_n$ being the largest prime $p_n \leq \sqrt{x}$.
Of course the floor function is untidy (and not multiplicative), so we want to replace it with:
$$\left[\frac{x}{d} \right] = \frac{x}{d} - \varepsilon(d)$$
Where $\varepsilon(d)$ is the error between the division and the floor function, bounded by $0 \leq \varepsilon(d) \leq (1 - \frac{1}{d})$. The bound on this error term does not depend on $x$; if we choose a large $x$ and divide it by $3$, the error $\varepsilon(3)$ will always be in $\{0, \frac13, \frac23\}$. Now we can rewrite $(1)$ as
$$\pi_L(x) = \sum_{d \mid P} \mu(d)\left(\frac{x}{d} - \varepsilon(d) \right) = \color{red}{x \sum_{d \mid P} \mu(d)\frac{1}{d}} - \color{blue}{\sum_{d \mid P} \mu(d) \varepsilon(d)} \tag{2}$$
The Euler product version of this is then:
$$\pi_L(x) = \color{red}{x \prod_{2 \leq p}^{p_n} \left(1 - \frac{1}{p}\right)} - \color{blue}{(-1)^n \prod_{2 \leq p}^{p_n} \varepsilon(p)} \tag{3}$$
Since $\varepsilon(1) = 0$, the Euler product of the second term results in simply $-\varepsilon(p)$ for each prime used, giving the alternating sign. Because we know $0 \leq \varepsilon(p) \leq (1 - \frac{1}{p}) < 1$, the total error term (in blue) should be bounded by $|\varepsilon| < 1$. We know this is not the case, and that this error term becomes large quickly. Therefore, I've done something wrong in the midst of my math.
Can someone tell me what I've done wrong? Is there some substitution that fails due to some non-obvious (to me) assumption? Is the problem that $\varepsilon(x)$ is not multiplicative? If not, it should still be bounded by $|\varepsilon| < 2^{n-1}$ (where $\varepsilon = \sum \mu(d)\varepsilon(d)$ in equation $(2)$), and on average much closer to $0$ than to the bounds, with the chance of being closer to $0$ higher as $x$ increases. Though maybe not provably closer? (Empirically, $|\varepsilon| < 1$ for $x < 60$.)
(The bound on non-multiplicative $\varepsilon$ comes from having alternating values of $\mu(d)$ as more primes are used. There are $\binom{n}{1} + \binom{n}{3} + \cdots = 2^{n-1}$ error terms with $\mu(d) = -1$, and $\binom{n}{2} + \binom{n}{4} + \cdots = 2^{n-1} - 1$ error terms with $\mu(d) = 1$. The largest possible error is then $-2^{n-1}$, if all the positive error terms are $0$. (We can safely assume there's no error term on $\frac{x}{1}$).)