How can I compute the probability that each set of n flips resulted in the same number of heads?

I've spent the last hour coming up with the general formula below. My primary question is: how would I compute the result for say, n = 10000? (As a last resort I could write C++ logic to compute the formula, but I imagine there's an easier free math tool out there).

Also, I do wonder if the general formula is correct.

So, for n flips, with a probability p of landing heads each flip.

p^{2n} + \sum_{k = 1}^n \left[ \left( \prod_{i=1}^k \left[\frac{n-(k-i)}{i} \right] \times p^{n-k} (1-p)^{k} \right)^2 \right]


If the general formula is well known or you could provide a link to it, I'd appreciate that.

Short explanation of how I came up with it: The \prod part calculates the binomial coefficient. So, if you flipped a coin 3 times, there are 3 ways it to get one head. That is multiplied by the probability of getting n-k heads and k tails.

Finally, I want to know the chance that n-k heads occurred in each set of flips, so I square it, then sum the whole thing from k=1 to n. The leading p^{2n} is there because the multiplicative binomial formula (the \prod part) does not work when k = 0.