I tried this question twice on math.stackexchange but got no answer so I decided to move it here.

Let $M$ be a smoth manifold. Then $$C^\infty(M):=\{f:M\longrightarrow \mathbb R; f\ \textrm{is smooth}\},$$ is an $\mathbb R$-algebra with the pointwise product.

If you don't know anything about smooth manifolds it doesn't matter, all that you need to know is that $C^\infty(M)$ is an $\mathbb R$-algebra for all that follows is purely algebraic.

Let us define the space of vector fields on $M$ by: $$\mathfrak{X}(M):=\mathsf{der}_{\mathbb R}\ C^\infty(M):=\{X\in \mathsf{End}_{\mathbb R}(C^\infty(M)): X(fg)=fX(g)+X(f)g\}.$$ This is a $C^\infty(M)$-module with $$(f\cdot X)(g):=f X(g),$$

where $fX(g)$ is the pointwise product of the functions $f$ and $X(g)$.

We can then define the space of $p$-forms on $M$ as the $C^\infty(M)$-module:

$$\Omega^p(M):=\mathsf{Hom}_{C^\infty(M)}(\Lambda^p \mathfrak{X}(M), C^\infty(M)).$$ The De Rham differential on $M$ is the degree one operator $$d: \Omega^p(M)\longrightarrow \Omega^{p+1}(M)$$ given by $$\begin{eqnarray*} d\varepsilon(X_{1}, \ldots, X_{p+1})&&:=\sum_{\sigma\in\mathsf{Sh}(1, p)} \mathsf{sgn}(\sigma) X_{\sigma(1)}(\varepsilon(X_{\sigma(2)}, \ldots, X_{\sigma(p+1)})\\ &&+\sum_{\sigma\in\mathsf{Sh}(2, p-1)}\mathsf{sgn}(\sigma) \varepsilon([X_{\sigma(1)}, X_{\sigma(2)}], X_{\sigma(3)}, \ldots, X_{\sigma(p+1)}), \end{eqnarray*}$$ where $[\cdot, \cdot]$ is the commutator of vector fields which is defined by $$[X, Y]:=X\circ Y-Y\circ X.$$

Notation: For integers $p, q\geq 1$ let us write $S(p, q)$ as the subset of permutations $\sigma$ of the set $\{1, \ldots, p+q\}$ such that $\sigma(1)<\ldots< \sigma(p)$ and $\sigma(p+1)<\ldots< \sigma(p+q)$. The elements of $\mathsf{Sh}(p, q)$ are known as $(p, q)$-shufles for obvious resons.

Now, we can define a product $$\wedge: \Omega^p(M)\times \Omega^q(M)\longrightarrow \Omega^{p+q}(M)$$ setting $$(\varepsilon\wedge \eta)(X_1, \ldots, X_{p+q}):=\sum_{\sigma\in\mathsf{Sh}(p, q)} \mathsf{sgn}(\sigma) \varepsilon(X_{\sigma(1)}, \ldots, X_{\sigma(p)}) \eta(X_{\sigma(p+1)}, \ldots, X_{\sigma(p+q)}).$$ Can anyone help me to prove $$d(\varepsilon\wedge \eta)=d\varepsilon\wedge \eta+(-1)^p \varepsilon\wedge d\eta,$$ for every $\varepsilon\in \Omega^p(M)$ and $\eta\in \Omega^q(M)$?

I know the property I want to show is a classical one but I can't seem to find the proof using this algebraic formulation. I've already asked this question before and got no answer.

However, by that time things were more obscure so I decided to update with this improved version hoping someone could help me.

Remark:

1. The cardinality of the set $\mathsf{Sh}(p, q)$ is $\binom{p+q}{q}$.

2. There is a bijection $\mathsf{Sh}(p, q)\simeq \mathsf{Sh}(p-1, q)\sqcup \mathsf{Sh}(p, q-1)$. This is clear since:

$$\mathsf{Sh}(p, q)=\{\sigma\in\mathsf{Sh}(p, q): \sigma(1)=1\}\sqcup\{\sigma\in\mathsf{Sh}(p, q): \sigma(p+1)=1\}.$$ In particular:

$$\mathsf{Sh}(p+1, q)\simeq \mathsf{Sh}(p, q)\sqcup \mathsf{Sh}(p+1, q-1)\quad \textrm{and}\quad \mathsf{Sh}(p, q+1)\simeq \mathsf{Sh}(p, q)\sqcup \mathsf{Sh}(p-1, q+1).$$

2. Heuristically, in light of $2)$, when we add the terms $d\varepsilon\wedge \eta$ and $\varepsilon\wedge d\eta$ we get things which can be put into $\sum_{\in \mathsf{Sh}(p, q)}$ which is desirable but the remaining terms $\sum_{\sigma\in \mathsf{Sh}(p+1, q-1)}$ and $\sum_{\sigma\in\mathsf{Sh}(p-1, q+1)}$ don't cancel for the cardinality of the sets $\mathsf{Sh}(p+1, q-1)$ and $\mathsf{Sh}(p-1, q+1)$ are different. Indeed, there is no hope things will cancel when we do $d\varepsilon\wedge \eta+\varepsilon\wedge d\eta$ for $d\varepsilon\wedge \eta$ consists of terms which derivate $\varepsilon(-)$ whereas $\varepsilon\wedge d\eta$ consists of terms which derivate $\eta(-)$.