[proofplan]
The proof is the finite inclusion-exclusion identity for the condition that no prime $p \in \mathcal P$ with $p<z$ divides $a$. We first rewrite the sifted count as a sum over divisors $d \mid P(z)$ weighted by the Möbius function. Substituting the assumed asymptotic formula for $|\mathcal A_d|$ separates the main term from the remainder term. Finally, multiplicativity of $\omega$ turns the divisor sum in the main term into the stated Euler product.
[/proofplan]
[step:Express the sifted count by Möbius inversion over divisors of $P(z)$]
Define the function $I:\mathcal A \to \mathbb Z$ by
\begin{align*}
I: \mathcal A &\to \mathbb Z, \qquad a \mapsto \sum_{d \mid \gcd(a,P(z))}\mu(d).
\end{align*}
Since $P(z)$ is squarefree, every divisor of $\gcd(a,P(z))$ is squarefree. If $\gcd(a,P(z))=1$, then $I(a)=\mu(1)=1$. If $\gcd(a,P(z))>1$, let $Q(a)$ denote the finite set of primes dividing $\gcd(a,P(z))$. Then
\begin{align*}
I(a)=\sum_{E \subseteq Q(a)}(-1)^{|E|}=(1-1)^{|Q(a)|}=0.
\end{align*}
Therefore $I(a)$ is exactly the indicator of the condition $\gcd(a,P(z))=1$. Summing over $a \in \mathcal A$ gives
\begin{align*}
S(\mathcal A,\mathcal P,z)=\sum_{a \in \mathcal A}\sum_{d \mid \gcd(a,P(z))}\mu(d).
\end{align*}
Because $\mathcal A$ is finite and $P(z)$ has finitely many divisors, we may interchange the two finite sums:
\begin{align*}
S(\mathcal A,\mathcal P,z)=\sum_{d \mid P(z)}\mu(d)|\mathcal A_d|.
\end{align*}
[guided]
The sifted condition $\gcd(a,P(z))=1$ says that none of the primes in the finite set $\{p \in \mathcal P : p<z\}$ divides $a$. Inclusion-exclusion encodes this condition using the Möbius function. Define the function $I:\mathcal A \to \mathbb Z$ by
\begin{align*}
I: \mathcal A &\to \mathbb Z, \qquad a \mapsto \sum_{d \mid \gcd(a,P(z))}\mu(d).
\end{align*}
We now verify that this is the indicator of the sifted condition. If $\gcd(a,P(z))=1$, then the only divisor in the sum is $d=1$, so $I(a)=\mu(1)=1$. If $\gcd(a,P(z))>1$, let $Q(a)$ be the finite set of primes dividing $\gcd(a,P(z))$. Since $P(z)$ is a product of distinct primes, every divisor of $\gcd(a,P(z))$ is obtained uniquely by choosing a subset $E \subseteq Q(a)$ and taking the product of the primes in $E$. For such a divisor, the Möbius value is $(-1)^{|E|}$. Hence
\begin{align*}
I(a)=\sum_{E \subseteq Q(a)}(-1)^{|E|}=(1-1)^{|Q(a)|}=0.
\end{align*}
Thus $I(a)=1$ exactly when $a$ survives the sieve, and $I(a)=0$ otherwise. Summing this indicator over all $a \in \mathcal A$ gives
\begin{align*}
S(\mathcal A,\mathcal P,z)=\sum_{a \in \mathcal A}\sum_{d \mid \gcd(a,P(z))}\mu(d).
\end{align*}
The sums are finite, because $\mathcal A$ is finite and $P(z)$ has finitely many divisors. Therefore we may interchange the order of summation. For a fixed divisor $d \mid P(z)$, the condition $d \mid \gcd(a,P(z))$ is equivalent to $d \mid a$, so the inner count is exactly $|\mathcal A_d|$. This gives
\begin{align*}
S(\mathcal A,\mathcal P,z)=\sum_{d \mid P(z)}\mu(d)|\mathcal A_d|.
\end{align*}
[/guided]
[/step]
[step:Substitute the divisor estimates and separate main and remainder terms]
Using the hypothesis $|\mathcal A_d|=\frac{\omega(d)}{d}X+r_d$ for every $d \mid P(z)$, the inclusion-exclusion formula becomes
\begin{align*}
S(\mathcal A,\mathcal P,z)=\sum_{d \mid P(z)}\mu(d)\left(\frac{\omega(d)}{d}X+r_d\right).
\end{align*}
Distributing the finite sum gives
\begin{align*}
S(\mathcal A,\mathcal P,z)=X\sum_{d \mid P(z)}\mu(d)\frac{\omega(d)}{d}+\sum_{d \mid P(z)}\mu(d)r_d.
\end{align*}
[/step]
[step:Factor the main divisor sum into the Euler product]
Define the function $f:\{d \in \mathbb N : d \mid P(z)\} \to \mathbb R$ by
\begin{align*}
f(d) := \mu(d)\frac{\omega(d)}{d}.
\end{align*}
The function $f$ is multiplicative on divisors of $P(z)$, because $\mu$ is multiplicative, $\omega$ is multiplicative on these divisors, and $d \mapsto 1/d$ is multiplicative. Also $f(1)=\mu(1)\omega(1)=1$ by the normalization hypothesis $\omega(1)=1$. Since the divisors of $P(z)$ are exactly the products of subsets of the finite set $\{p \in \mathcal P : p<z\}$, the finite [Euler product identity](/theorems/1694) gives
\begin{align*}
\sum_{d \mid P(z)}\mu(d)\frac{\omega(d)}{d}=\prod_{\substack{p<z, p\in \mathcal P}}\left(1+\mu(p)\frac{\omega(p)}{p}\right).
\end{align*}
For every prime $p$, $\mu(p)=-1$, so
\begin{align*}
\sum_{d \mid P(z)}\mu(d)\frac{\omega(d)}{d}=\prod_{\substack{p<z, p\in \mathcal P}}\left(1-\frac{\omega(p)}{p}\right).
\end{align*}
Substituting this identity into the expression from the previous step yields
\begin{align*}
S(\mathcal A,\mathcal P,z)=X\prod_{\substack{p<z, p\in \mathcal P}}\left(1-\frac{\omega(p)}{p}\right)+\sum_{d \mid P(z)}\mu(d)r_d.
\end{align*}
This is the claimed exact sieve expansion with remainders.
[/step]
