Random Gaussian Variables

Jarrett's suggestion of using a Box-Muller transform is good for a quick-and-dirty solution. A simple implementation:

Random rand = new Random(); //reuse this if you are generating many
double u1 = 1.0-rand.NextDouble(); //uniform(0,1] random doubles
double u2 = 1.0-rand.NextDouble();
double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) *
             Math.Sin(2.0 * Math.PI * u2); //random normal(0,1)
double randNormal =
             mean + stdDev * randStdNormal; //random normal(mean,stdDev^2)

I created a request for such a feature on Microsoft Connect. If this is something you're looking for, please vote for it and increase its visibility.

https://connect.microsoft.com/VisualStudio/feedback/details/634346/guassian-normal-distribution-random-numbers

This feature is included in the Java SDK. Its implementation is available as part of the documentation and is easily ported to C# or other .NET languages.

If you're looking for pure speed, then the Zigorat Algorithm is generally recognised as the fastest approach.

I'm not an expert on this topic though -- I came across the need for this while implementing a particle filter for my RoboCup 3D simulated robotic soccer library and was surprised when this wasn't included in the framework.

In the meanwhile, here's a wrapper for Random that provides an efficient implementation of the Box Muller polar method:

public sealed class GaussianRandom
{
    private bool _hasDeviate;
    private double _storedDeviate;
    private readonly Random _random;

    public GaussianRandom(Random random = null)
    {
        _random = random ?? new Random();
    }

    /// <summary>
    /// Obtains normally (Gaussian) distributed random numbers, using the Box-Muller
    /// transformation.  This transformation takes two uniformly distributed deviates
    /// within the unit circle, and transforms them into two independently
    /// distributed normal deviates.
    /// </summary>
    /// <param name="mu">The mean of the distribution.  Default is zero.</param>
    /// <param name="sigma">The standard deviation of the distribution.  Default is one.</param>
    /// <returns></returns>
    public double NextGaussian(double mu = 0, double sigma = 1)
    {
        if (sigma <= 0)
            throw new ArgumentOutOfRangeException("sigma", "Must be greater than zero.");

        if (_hasDeviate)
        {
            _hasDeviate = false;
            return _storedDeviate*sigma + mu;
        }

        double v1, v2, rSquared;
        do
        {
            // two random values between -1.0 and 1.0
            v1 = 2*_random.NextDouble() - 1;
            v2 = 2*_random.NextDouble() - 1;
            rSquared = v1*v1 + v2*v2;
            // ensure within the unit circle
        } while (rSquared >= 1 || rSquared == 0);

        // calculate polar tranformation for each deviate
        var polar = Math.Sqrt(-2*Math.Log(rSquared)/rSquared);
        // store first deviate
        _storedDeviate = v2*polar;
        _hasDeviate = true;
        // return second deviate
        return v1*polar*sigma + mu;
    }
}

Math.NET Iridium also claims to implement "non-uniform random generators (normal, poisson, binomial, ...)".

Expanding off of @Noakes and @Hameer's answers, I have also implemented a 'Gaussian' class, but to simplify memory space, I made it a child of the Random class so that you can also call the basic Next(), NextDouble(), etc from the Gaussian class as well without having to create an additional Random object to handle it. I also eliminated the _available, and _nextgauss global class properties, as I didn't see them as necessary since this class is instance based, it should be thread-safe, if you give each thread its own Gaussian object. I also moved all of the run-time allocated variables out of the function and made them class properties, this will reduce the number of calls to the memory manager since the 4 doubles should theoretically never be de-allocated until the object is destroyed.

public class Gaussian : Random
{

    private double u1;
    private double u2;
    private double temp1;
    private double temp2;

    public Gaussian(int seed):base(seed)
    {
    }

    public Gaussian() : base()
    {
    }

    /// <summary>
    /// Obtains normally (Gaussian) distrubuted random numbers, using the Box-Muller
    /// transformation.  This transformation takes two uniformly distributed deviates
    /// within the unit circle, and transforms them into two independently distributed normal deviates.
    /// </summary>
    /// <param name="mu">The mean of the distribution.  Default is zero</param>
    /// <param name="sigma">The standard deviation of the distribution.  Default is one.</param>
    /// <returns></returns>

    public double RandomGauss(double mu = 0, double sigma = 1)
    {
        if (sigma <= 0)
            throw new ArgumentOutOfRangeException("sigma", "Must be greater than zero.");

        u1 = base.NextDouble();
        u2 = base.NextDouble();
        temp1 = Math.Sqrt(-2 * Math.Log(u1));
        temp2 = 2 * Math.PI * u2;

        return mu + sigma*(temp1 * Math.Cos(temp2));
    }
}

The other implementations are a little too heady for me. My implementation relies on the fact that throwing three six sided dice (D6) and recording the sum of the dice will creates a normal distribution. Just replace the three D6 with three calls to your your language's "rand()" function and you're just about there. In the implementation below, I subtract by two divided by the number of "dice," so 1.5 in this case.

using System;

public class CustomMath
{
    private static readonly Random _random = new Random();
    public static double GaussianRandom() =>
        _random.NextDouble() + _random.NextDouble() + _random.NextDouble() - 1.5;
}