Maths Spl

Theory
It represent the results of two coin tosses (e.g. half the difference between the total number of heads and tails) with the expression A + A*cos (Pi*n/N). N is the number of samples in your buffer - a binomial expansion of even order O will have O+1 coefficients and require a buffer of N = O/2 + 1 samples - n is the sample number being generated, and A is a scale factor that will usually be either 2 (for generating binomial coefficients) or 0.5 (for generating a binomial probability distribution).
In frequency, this expression resembles the binomial distribution of those two coin tosses. There are three symmetrical spikes at positions corresponding to the number (heads-tails)/2. Since, modelling the overall probability distribution of independent events requires convolving their distributions. We want to convolve our expression in the frequency domain, which is equivalent to multiplication in the time domain.
In other words, by raising our cosine expression for the result of two tosses to a power (e.g. to simulate 500 tosses, raise it to the power of 250 since it already represents a pair), we can arrange for the binomial distribution for a large number to appear in the frequency domain. Since, this is all real and even, we can substitute the DCT-I for the DFT to improve efficiency.
Accuracy:
There's a limit to how high O can be before accumulated floating-point rounding errors rob you of accurate integer values for the coefficients, but I'd guess the number is pretty high. Double-precision floating-point can represent 53-bit integers with complete accuracy, and I'm going to ignore the rounding loss involved in the use of pow() because the generating expression will take place in FP registers.