[proofplan]
With $P(z)$ and $S(A,P,z)$ defined in the statement, the proof rewrites the condition $\gcd(a,P(z)) = 1$ as an exact divisor sum using the Möbius function $\mu$. For each fixed $a \in A$, the sum $\sum_{d \mid \gcd(a,P(z))}\mu(d)$ is an indicator for the event that no prime divisor of the positive integer $P(z)$ divides $a$. Summing this identity over $A$ and exchanging finite sums gives the desired formula, because the condition $d \mid \gcd(a,P(z))$ is equivalent to $d \mid P(z)$ and $a \in A_d$.
[/proofplan]
[step:Prove the Möbius divisor sum detects coprimality]
Let $m$ be a positive integer. We claim that the divisor sum satisfies
\begin{align*}
\sum_{d \mid m} \mu(d) = 1
\end{align*}
when $m=1$, and satisfies
\begin{align*}
\sum_{d \mid m} \mu(d) = 0
\end{align*}
when $m>1$.
If $m=1$, the only positive divisor of $m$ is $1$, and $\mu(1)=1$, so the sum is $1$.
Assume $m>1$. Let $p_1,\dots,p_r$ be the distinct prime divisors of $m$, where $r \ge 1$. Since $\mu(d)=0$ whenever $d$ is not squarefree, only squarefree divisors of $m$ contribute to the sum. Each squarefree divisor contributing to the sum has the form
\begin{align*}
d = \prod_{i \in I} p_i
\end{align*}
for a unique subset $I \subset \{1,\dots,r\}$, and then $\mu(d)=(-1)^{|I|}$. Therefore
\begin{align*}
\sum_{d \mid m} \mu(d) = \sum_{I \subset \{1,\dots,r\}} (-1)^{|I|}.
\end{align*}
By the [binomial theorem](/theorems/750) applied to $(1-1)^r$,
\begin{align*}
\sum_{I \subset \{1,\dots,r\}} (-1)^{|I|} = (1-1)^r = 0.
\end{align*}
This proves the claimed divisor-sum identity.
[guided]
The point of this step is to isolate the elementary Möbius cancellation that drives the whole sieve identity. Let $m$ be a positive integer. We want to compute
\begin{align*}
\sum_{d \mid m} \mu(d).
\end{align*}
If $m=1$, then the divisor set is exactly $\{1\}$. Since the Möbius function satisfies $\mu(1)=1$, the sum is
\begin{align*}
\sum_{d \mid 1} \mu(d) = \mu(1) = 1.
\end{align*}
Now suppose $m>1$. Let $p_1,\dots,p_r$ denote the distinct prime divisors of $m$. Since $m>1$, we have $r \ge 1$. The Möbius function vanishes on every non-squarefree integer, so non-squarefree divisors of $m$ contribute $0$ to the sum. Thus we only need to count squarefree divisors.
Every squarefree divisor $d$ of $m$ is obtained by choosing some subset $I \subset \{1,\dots,r\}$ and multiplying the corresponding primes:
\begin{align*}
d = \prod_{i \in I} p_i.
\end{align*}
For this divisor, the Möbius value is $\mu(d)=(-1)^{|I|}$, because $d$ is a product of exactly $|I|$ distinct primes. Hence
\begin{align*}
\sum_{d \mid m} \mu(d) = \sum_{I \subset \{1,\dots,r\}} (-1)^{|I|}.
\end{align*}
The last sum is the binomial expansion of $(1-1)^r$. Since $r \ge 1$,
\begin{align*}
\sum_{I \subset \{1,\dots,r\}} (-1)^{|I|} = (1-1)^r = 0.
\end{align*}
Therefore the divisor sum equals $1$ when $m=1$ and equals $0$ when $m>1$.
[/guided]
[/step]
[step:Apply the divisor sum to each element of $A$]
Define the sifted subset $E_z \subset A$ by
\begin{align*}
E_z := \{a \in A : \gcd(a,P(z))=1\}.
\end{align*}
Thus $S(A,P,z)=|E_z|$ by the definition of the sifted counting function in the statement. For each $a \in A$, define the positive integer
\begin{align*}
m_a := \gcd(a,P(z)).
\end{align*}
By the identity proved above, the sum
\begin{align*}
\sum_{d \mid m_a} \mu(d)
\end{align*}
equals $1$ when $a \in E_z$ and equals $0$ when $a \notin E_z$. Define the indicator function
\begin{align*}
\mathbb{1}_{E_z}: A &\to \{0,1\}
\end{align*}
by declaring $\mathbb{1}_{E_z}(a)=1$ if $a \in E_z$ and $\mathbb{1}_{E_z}(a)=0$ otherwise. Hence, for every $a \in A$,
\begin{align*}
\mathbb{1}_{E_z}(a) = \sum_{d \mid \gcd(a,P(z))} \mu(d).
\end{align*}
Summing this identity over all $a \in A$ gives
\begin{align*}
S(A,P,z) = |E_z| = \sum_{a \in A} \mathbb{1}_{E_z}(a) = \sum_{a \in A} \sum_{d \mid \gcd(a,P(z))} \mu(d).
\end{align*}
[/step]
[step:Exchange the finite sums and identify the sets $A_d$]
Because $A$ is finite and the sifting modulus $P(z)$ is a positive integer, $P(z)$ has finitely many positive divisors. Hence all sums below are finite. The condition $d \mid \gcd(a,P(z))$ is equivalent to the two simultaneous conditions $d \mid P(z)$ and $d \mid a$. Therefore
\begin{align*}
\sum_{a \in A} \sum_{d \mid \gcd(a,P(z))} \mu(d) = \sum_{d \mid P(z)} \mu(d)|\{a \in A : d \mid a\}|.
\end{align*}
By the definition of $A_d$, the inner cardinality is $|A_d|$. Hence
\begin{align*}
S(A,P,z) = \sum_{d \mid P(z)} \mu(d)|A_d|.
\end{align*}
This is the desired exact sifting identity.
[/step]