Transform uniform distribution to normal distribution using Lindeberg–Lévy CLT Currently i am developing a game which is based on many computations of random values and therefore i have implemented many algorithms like the Mersenne-Twister etc. Unfortunately, all generators return uniformly distributed values and i want to modify this to make the min/max (and their surrounding) values rarer than the other ones. Researches brought me up to the point that i have been using now the Lindeberg–Lévy Central Limit Theorem:

Theorem. Suppose {$X_i$} is a sequence of iid random variables with $\mathbb{E}[X_i] = \mu$ and $Var[X_i] = \sigma^2$. Let $S_n=\frac{1}{n}(X_1+X_2+\ldots+X_n)$. Then as $n$ approaches infinity, the random variable $\sqrt{n}(S_n − \mu)$ converges in distribution to a normal $\mathcal{N}(0, \sigma^2)$:
  $$\sqrt{n}\cdot\left(\left(\frac{1}{n}\sum\limits_{i=1}^nX_i\right)-\mu\right)\overset{d}{\rightarrow}\mathcal{N}(0, \sigma^2)$$

My C++ code works fine... however the median is set to 0 and i want to know how i can change the median and the range in combination to get e.g. values from $[1;100]$ as a result will have the property, that the median ones (50 and surrounding) will have a higher probability to be returned than [1;5] or [96;100]. (Hope you understand what i am thinking right now!)
Concrete example
Throwing 10.000 "unfair" dice currently leads to this distribution: (8836,729,341,86,8,1) however i want something that looks like this: (18,341,4418,4568,341,56) 
Current algorithm (old one based on the theorem!)
    n <- 32;
    S <- 0.0;
    for i = 0 to n
        S <- S + mersenneTwister(min, max)
    S <- S / n
    µ <- (min + max) / 2.0
    randval <- sqrt(n) * (S - µ)
    return max(min, floor(randval) mod max)

It seems to me, that something it not implemented correctly...
Solution (i am not yet able to answer my own questions...)
Finally i have found a very simple solution which suits my needs and is no overkill at all.

Marsaglia polar method (more information here)
  
  
*
  
*Generate two iid random values $u_1$ and $u_2$ in [-1;1].
  
*Compute $q = u_1^2 + u_2^2$. In case that $q=0$ or $q > 1$ return to step 1.
  
*Compute $p=\sqrt{-2\cdot \ln(q) / q}$.
  
*$x_{1/2} = u_{1/2}\cdot p$ represent two independant normal distributed random numbers.
  
*If $x\sim\mathcal{N}(0,1)$, then $a\cdot x+b\sim\mathcal{N}(b,a^2)$.
  

Exactly what i need! For throwing my dice i let $a=1$ and $b=\mu=(min+max)/2$.
@cardinal: Thanks for this idea!
 A: The Marsaglia Polar Method is an algorithmic improvement on the Box-Muller Transform. The Box-Muller Transform uses two uniform random variables, $(u,v)$, in $[0,1]\times[0,1]$ to generate two normally distributed variables in $\mathbb{R}\times\mathbb{R}$:
$$
(\cos(2\pi u),\sin(2\pi u))\sqrt{-2\log(v)}
$$
The Marsaglia Polar Method uses the fact that for a uniformly distributed point, $(x,y)$, in the unit disk, $x^2+y^2$ is uniformly distributed in $[0,1]$. Thus, $x^2+y^2$ can take the place of $v$ in the Box-Muller Transform. Since $\frac{(x,y)}{\sqrt{x^2+y^2}}$ is uniformly distributed on the unit circle, it can take the place of $(\cos(2\pi u),\sin(2\pi u))$ in the Box-Muller Transform. With these substitutions, we get that
$$
\frac{(x,y)}{\sqrt{x^2+y^2}}\sqrt{-2\log(x^2+y^2)}
$$
is a normally distributed random variable. The Marsaglia Polar Method uses a Monte-Carlo method to generate a uniformly distributed $(x,y)$ in the unit disk; that is, it picks points in $[-1,1]\times[-1,1]$ until one lands in the unit disk (which happens with probability $\frac{\pi}{4}$). To prevent illegal values, the unlikely value $(0,0)$ is ignored.
I say that the Marsaglia Polar Method is an algorithmic improvement on the Box-Muller Transform because it avoids the computation of a trigonometric function or two, but at the cost of generating $\frac{4}{\pi}$ uniform pairs, on average, per each normal pair.
A: The normal distributed random variable can be just shifted in order to obtain the mean you need. E.g. in you algorithm it holds that
$$
\sqrt{n}(S_n-\mu)+a \to \mathcal N(a,\sigma^2).
$$
But if you just want to obtain the normal distributed random variable from the uniform distribution, why don't just use a random variable $U$ which is uniformly distributed on $[0,1]$ - then $F^{-1}(U)$ will be distributed as a random variable with cumultaive distribution function $F$. 
For example, for normal random variable,
$$
F(x;\mu,\sigma^2) = \frac1{\sigma\sqrt{2\pi}}\int\limits_{-\infty}^x\mathrm e^{-(t-\mu)^2/2\sigma^2}\mathrm dt
$$
and I believe this function you can easily calculate numerically as well as its inverse. Moreover, if you are interested in discrete value, you may want to use Binomial Distribution which is very similar to the normal.
