This answer is a bit late, but the problem has quite an elegant solution based on the Fourier transform and I wanted to add it.

The action of the BSC can be modeled as a$\pmod 2$ sum of the data bit, $X$, and a random noise bit, $N$, who's probability distribution is $P(N=1)=p$ and $P(N=0)=1-p$,

$$Y = X \oplus N.$$

The distribution over $Y$ (the output) is the conditional distribution $P(Y|X,N)$ that is a [convolution of the quantities][1] $P(X)$ and $P(N)$,

$$P(Y=y) = \sum_{x \in \{0,1\}} P(X=x)P(N=x\oplus y)$$

a linear transformation which, by the convolution theorem, is diagonalized by the 2x2 Hadamard DFT matrix,

$$H = \begin{bmatrix}
1 & 1 \\
1 & -1
\end{bmatrix},$$

and becomes a product of their spectra in the transform domain.

Exponentiation in the transform domain followed by an inverse transform,

$$H^{-1}\left(H\begin{bmatrix}1-p\\p\end{bmatrix}\right)^L = \begin{bmatrix}{(1 - 2p)^L}\\{1 - (1 - 2p)^L}\end{bmatrix},$$

is the distribution over the sum of $L$ identically distributed noise bits.

  [1]: https://en.wikipedia.org/wiki/Convolution_of_probability_distributions