The formulas are straightforward but they are not as simple as intimated in the question.
Let $Y$ be the previous EWMA and let $X = x_n$, which is presumed independent of $Y$. By definition, the new weighted average is $Z = \alpha X + (1 - \alpha)Y$ for a constant value $\alpha$. For notational convenience, set $\beta = 1-\alpha$. Let $F$ denote the CDF of a random variable and $\phi$ denote its moment generating function, so that
$$\phi_X(t) = \mathbb{E}_F[\exp(t X)] = \int_\mathbb{R}{\exp(t x) dF_X(x)}.$$
With Kendall and Stuart, let $\mu_k^{'}(Z)$ denote the non-central moment of order $k$ for the random variable $Z$; that is, $\mu_k^{'}(Z) = \mathbb{E}[Z^k]$. The skewness and kurtosis are expressible in terms of the $\mu_k^{'}$ for $k = 1,2,3,4$; for example, the skewness is defined as $\mu_3 / \mu_2^{3/2}$ where
$$\mu_3 = \mu_3^{'} - 3 \mu_2^{'}\mu_1^{'} + 2{\mu_1^{'}}^3 \text{ and }\mu_2 = \mu_2^{'} - {\mu_1^{'}}^2$$
are the third and second central moments, respectively.
By standard elementary results,
$$\eqalign{
&1 + \mu_1^{'}(Z) t + \frac{1}{2!} \mu_2^{'}(Z) t^2 + \frac{1}{3!} \mu_3^{'}(Z) t^3 + \frac{1}{4!} \mu_4^{'}(Z) t^4 +O(t^5) \cr
&= \phi_Z(t) \cr
&= \phi_{\alpha X}(t) \phi_{\beta Y}(t) \cr
&= \phi_X(\alpha t) \phi_Y(\beta t) \cr
&= (1 + \mu_1^{'}(X) \alpha t + \frac{1}{2!} \mu_2^{'}(X) \alpha^2 t^2 + \cdots)
(1 + \mu_1^{'}(Y) \beta t + \frac{1}{2!} \mu_2^{'}(Y) \beta^2 t^2 + \cdots).
}
$$
To obtain the desired non-central moments, multiply the latter power series through fourth order in $t$ and equate the result term-by-term with the terms in $\phi_Z(t)$.
I agree with @NickCox : I think the mistake is in the first line of your post, where you define "bimodality coefficient". I Googled and found Pfister et al (which references SAS/STAT from 1990). That paper indicates problems with BC that are quite similar to the ones you found and recommends Hartigan's dip test, instead of BC (or in addition to it). The dip test is available in R through the diptest package. In addition, the kurtosis in the formula is supposed to be excess kurtosis and you appear to not have adjusted for that (although I am not certain of this)
The SAS documentation also mentions problems with BC, in particular
Very heavy-tailed distributions have small values of regardless of the
number of modes.
The long tail of your second distribution is probably lowering the value of BC.
In short, the problem is in the formula, not in your code. There is, as far as I know, no perfect measure of the number of modes.
Best Answer
Yes, examples with skewness and excess kurtosis both zero are relatively easy to construct. (Indeed examples (a) to (d) below also have Pearson mean-median skewness 0)
(a) For example, in this answer an example is given by taking a 50-50 mixture of a gamma variate, (which I call $X$), and the negative of a second one, which has a density that looks like this:
Clearly the result is symmetric and not normal. The scale parameter is unimportant here, so we can make it 1. Careful choice of the shape parameter of the gamma yields the required kurtosis:
The variance of this double-gamma ($Y$) is easy to work out in terms of the gamma variate it's based on: $\text{Var}(Y)=E(X^2)=\text{Var}(X)+E(X)^2=\alpha+\alpha^2$.
The fourth central moment of the variable $Y$ is the same as $E(X^4)$, which for a gamma($\alpha$) is $\alpha(\alpha+1)(\alpha+2)(\alpha+3)$
As a result the kurtosis is $\frac{\alpha(\alpha+1)(\alpha+2)(\alpha+3)}{\alpha^2(\alpha+1)^2}=\frac{(\alpha+2)(\alpha+3)}{\alpha(\alpha+1)}$. This is $3$ when $(\alpha+2)(\alpha+3)=3\alpha(\alpha+1)$, which happens when $\alpha=(\sqrt{13}+1)/2\approx 2.303$.
(b) We could also create an example as a scale mixture of two uniforms. Let $U_1\sim U(-1,1)$ and let $U_2\sim U(-a,a)$, and let $M=\frac12 U_1+\frac12 U_2$. Clearly by considering that $M$ is symmetric and has finite range, we must have $E(M)=0$; the skewness will also be 0 and central moments and raw moments will be the same.
$\text{Var}(M)=E(M^2)=\frac12\text{Var}(U1)+\frac12\text{Var}(U_2)=\frac16[1+a^2]$.
Similarly, $E(M^4)=\frac{1}{10} (1+a^4)$ and so the kurtosis is $\frac{\frac{1}{10} (1+a^4)}{[\frac16 (1+a^2)]^2}=3.6\frac{1+a^4}{(1+a^2)^2}$
If we choose $a=\sqrt{5+\sqrt{24}}\approx 3.1463$, then kurtosis is 3, and the density looks like this:
(c) here's a fun example. Let $X_i\stackrel{_\text{iid}}{\sim}\text{Pois}(\lambda)$, for $i=1,2$.
Let $Y$ be a 50-50 mixture of $\sqrt{X_1}$ and $-\sqrt{X_2}$:
by symmetry $E(Y)=0$ (we also need $E(|Y|)$ to be finite but given $E(X_1)$ is finite, we have that)
$Var(Y)=E(Y^2)=E(X_1)=\lambda$
by symmetry (and the fact that the absolute 3rd moment exists) skew=0
4th moment: $E(Y^4) = E(X_1^2) = \lambda+\lambda^2$
kurtosis = $\frac{\lambda+\lambda^2}{\lambda^2}= 1+1/\lambda$
so when $\lambda=\frac12$, kurtosis is 3. This is the case illustrated above.
(d) all my examples so far have been symmetric, since symmetric answers are easier to create -- but asymmetric solutions are also possible. Here's a discrete example.
As you see, none of these examples look particularly "normal". It would be a simple matter to make any number of discrete, continuous or mixed variables with the same properties. While most of my examples were constructed as mixtures, there's nothing special about mixtures, other than they're often a convenient way to make distributions with properties the way you want, a bit like building things with Lego.
This answer gives some additional details on kurtosis that should make some of the considerations involved in constructing other examples a little clearer.
You could match more moments in similar fashion, though it requires more effort to do so. However, because the MGF of the normal exists, you can't match all integer moments of a normal with some non-normal distribution, since that would mean their MGFs match, implying the second distribution was normal as well.