diff --git a/src/FSharp/Distributions.fs b/src/FSharp/Distributions.fs
index 266c7839..83ba715c 100644
--- a/src/FSharp/Distributions.fs
+++ b/src/FSharp/Distributions.fs
@@ -84,8 +84,8 @@ module Sample =
     let discreteUniformSeq lower upper (rng:System.Random) = DiscreteUniform.Samples(rng, lower, upper)
 
     /// Erlang with shape (k) and rate or inverse scale (λ).
-    let erlang shape rate (rng:System.Random) = Erlang.Sample(rng, shape, rate)
-    let erlangSeq shape rate (rng:System.Random) = Erlang.Samples(rng, shape, rate)
+    let erlang (shape:int) rate (rng:System.Random) = Erlang.Sample(rng, shape, rate)
+    let erlangSeq (shape:int) rate (rng:System.Random) = Erlang.Samples(rng, shape, rate)
 
     /// Exponential with rate (λ).
     let exponential rate (rng:System.Random) = Exponential.Sample(rng, rate)
diff --git a/src/Numerics/Distributions/Bernoulli.cs b/src/Numerics/Distributions/Bernoulli.cs
index c49a6243..3f454bcd 100644
--- a/src/Numerics/Distributions/Bernoulli.cs
+++ b/src/Numerics/Distributions/Bernoulli.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -288,6 +290,23 @@ namespace MathNet.Numerics.Distributions
             return 0;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double p)
+        {
+            var uniform = rnd.NextDoubles(values.Length);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = uniform[i] < p ? 1 : 0;
+                }
+            });
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double p)
+        {
+            return rnd.NextDoubleSequence().Select(r => r < p ? 1 : 0);
+        }
+
         /// <summary>
         /// Samples a Bernoulli distributed random variable.
         /// </summary>
@@ -297,6 +316,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _p);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _p);
+        }
+
         /// <summary>
         /// Samples an array of Bernoulli distributed random variables.
         /// </summary>
@@ -318,6 +345,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, p);
         }
 
@@ -330,10 +358,22 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p);
-            }
+
+            return SamplesUnchecked(rnd, p);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The probability (p) of generating one. Range: 0 ≤ p ≤ 1.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, int[] values, double p)
+        {
+            if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, p);
         }
 
         /// <summary>
@@ -344,6 +384,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, p);
         }
 
@@ -355,12 +396,21 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, p);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The probability (p) of generating one. Range: 0 ≤ p ≤ 1.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(int[] values, double p)
+        {
+            if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, p);
+        }
     }
 }
diff --git a/src/Numerics/Distributions/Beta.cs b/src/Numerics/Distributions/Beta.cs
index 26bd6258..b4523205 100644
--- a/src/Numerics/Distributions/Beta.cs
+++ b/src/Numerics/Distributions/Beta.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -33,6 +33,7 @@ using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.RootFinding;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -142,7 +143,7 @@ namespace MathNet.Numerics.Distributions
         {
             get
             {
-                if (_shapeA == 0.0 && _shapeB == 0.0)  return 0.5;
+                if (_shapeA == 0.0 && _shapeB == 0.0) return 0.5;
                 if (_shapeA == 0.0) return 0.0;
                 if (_shapeB == 0.0) return 1.0;
 
@@ -352,16 +353,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _shapeA, _shapeB);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _shapeA, _shapeB);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Beta distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _shapeA, _shapeB);
-            }
+            return SamplesUnchecked(_random, _shapeA, _shapeB);
         }
 
         /// <summary>
@@ -375,7 +381,29 @@ namespace MathNet.Numerics.Distributions
         {
             var x = Gamma.SampleUnchecked(rnd, a, 1.0);
             var y = Gamma.SampleUnchecked(rnd, b, 1.0);
-            return x / (x + y);
+            return x/(x + y);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double a, double b)
+        {
+            var y = new double[values.Length];
+            Gamma.SamplesUnchecked(rnd, values, a, 1.0);
+            Gamma.SamplesUnchecked(rnd, y, b, 1.0);
+            CommonParallel.For(0, values.Length, 4096, (aa, bb) =>
+            {
+                for (int i = aa; i < bb; i++)
+                {
+                    values[i] = values[i]/(values[i] + y[i]);
+                }
+            });
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double a, double b)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, a, b);
+            }
         }
 
         /// <summary>
@@ -421,8 +449,8 @@ namespace MathNet.Numerics.Distributions
             if (b == 0.0) return x == 1.0 ? Double.PositiveInfinity : 0.0;
             if (a == 1.0 && b == 1.0) return 1.0;
 
-            var bb = SpecialFunctions.Gamma(a + b) / (SpecialFunctions.Gamma(a) * SpecialFunctions.Gamma(b));
-            return bb * Math.Pow(x, a - 1.0) * Math.Pow(1.0 - x, b - 1.0);
+            var bb = SpecialFunctions.Gamma(a + b)/(SpecialFunctions.Gamma(a)*SpecialFunctions.Gamma(b));
+            return bb*Math.Pow(x, a - 1.0)*Math.Pow(1.0 - x, b - 1.0);
         }
 
         /// <summary>
@@ -559,10 +587,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (a < 0.0 || b < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, a, b);
-            }
+            return SamplesUnchecked(rnd, a, b);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="a">The α shape parameter of the Beta distribution. Range: α ≥ 0.</param>
+        /// <param name="b">The β shape parameter of the Beta distribution. Range: β ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double a, double b)
+        {
+            if (a < 0.0 || b < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, a, b);
         }
 
         /// <summary>
@@ -573,7 +613,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double a, double b)
         {
-            return Sample(SystemRandomSource.Default, a, b);
+            if (a < 0.0 || b < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, a, b);
         }
 
         /// <summary>
@@ -584,7 +626,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double a, double b)
         {
-            return Samples(SystemRandomSource.Default, a, b);
+            if (a < 0.0 || b < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, a, b);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="a">The α shape parameter of the Beta distribution. Range: α ≥ 0.</param>
+        /// <param name="b">The β shape parameter of the Beta distribution. Range: β ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double a, double b)
+        {
+            if (a < 0.0 || b < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, a, b);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Binomial.cs b/src/Numerics/Distributions/Binomial.cs
index 028350f5..9f019613 100644
--- a/src/Numerics/Distributions/Binomial.cs
+++ b/src/Numerics/Distributions/Binomial.cs
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -222,8 +223,8 @@ namespace MathNet.Numerics.Distributions
         {
             get
             {
-                if (_p == 1.0) return new [] {_trials};
-                if (_p == 0.0) return new [] {0};
+                if (_p == 1.0) return new[] { _trials };
+                if (_p == 0.0) return new[] { 0 };
 
                 double td = (_trials + 1)*_p;
                 int t = (int)Math.Floor(td);
@@ -336,6 +337,32 @@ namespace MathNet.Numerics.Distributions
             return k;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double p, int n)
+        {
+            var uniform = rnd.NextDoubles(values.Length*n);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    int k = i*n;
+                    int sum = 0;
+                    for (int j = 0; j < n; j++)
+                    {
+                        sum += uniform[k + j] < p ? 1 : 0;
+                    }
+                    values[i] = sum;
+                }
+            });
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double p, int n)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, p, n);
+            }
+        }
+
         /// <summary>
         /// Samples a Binomially distributed random variable.
         /// </summary>
@@ -345,16 +372,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _p, _trials);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _p, _trials);
+        }
+
         /// <summary>
         /// Samples an array of Binomially distributed random variables.
         /// </summary>
         /// <returns>a sequence of successes in N trials.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _p, _trials);
-            }
+            return SamplesUnchecked(_random, _p, _trials);
         }
 
         /// <summary>
@@ -367,6 +399,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double p, int n)
         {
             if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, p, n);
         }
 
@@ -380,10 +413,23 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double p, int n)
         {
             if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p, n);
-            }
+
+            return SamplesUnchecked(rnd, p, n);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The success probability (p) in each trial. Range: 0 ≤ p ≤ 1.</param>
+        /// <param name="n">The number of trials (n). Range: n ≥ 0.</param>
+        /// <returns>a sequence of successes in <paramref name="n"/> trials.</returns>
+        public static void Samples(System.Random rnd, int[] values, double p, int n)
+        {
+            if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, p, n);
         }
 
         /// <summary>
@@ -395,6 +441,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double p, int n)
         {
             if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, p, n);
         }
 
@@ -407,11 +454,22 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double p, int n)
         {
             if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p, n);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, p, n);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The success probability (p) in each trial. Range: 0 ≤ p ≤ 1.</param>
+        /// <param name="n">The number of trials (n). Range: n ≥ 0.</param>
+        /// <returns>a sequence of successes in <paramref name="n"/> trials.</returns>
+        public static void Samples(int[] values, double p, int n)
+        {
+            if (!(p >= 0.0 && p <= 1.0 && n >= 0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, p, n);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Categorical.cs b/src/Numerics/Distributions/Categorical.cs
index 7acee836..ad9bd6c5 100644
--- a/src/Numerics/Distributions/Categorical.cs
+++ b/src/Numerics/Distributions/Categorical.cs
@@ -34,6 +34,7 @@ using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.Statistics;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -232,7 +233,7 @@ namespace MathNet.Numerics.Distributions
                 var sum = 0.0;
                 for (int i = 0; i < _pmfNormalized.Length; i++)
                 {
-                    sum += i * _pmfNormalized[i];
+                    sum += i*_pmfNormalized[i];
                 }
                 return sum;
             }
@@ -372,7 +373,7 @@ namespace MathNet.Numerics.Distributions
                 return 1.0;
             }
 
-            return _cdfUnnormalized[(int) Math.Floor(x)]/_cdfUnnormalized[_cdfUnnormalized.Length - 1];
+            return _cdfUnnormalized[(int)Math.Floor(x)]/_cdfUnnormalized[_cdfUnnormalized.Length - 1];
         }
 
         /// <summary>
@@ -388,7 +389,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentOutOfRangeException("probability");
             }
 
-            var denormalizedProbability = probability * _cdfUnnormalized[_cdfUnnormalized.Length - 1];
+            var denormalizedProbability = probability*_cdfUnnormalized[_cdfUnnormalized.Length - 1];
             int idx = Array.BinarySearch(_cdfUnnormalized, denormalizedProbability);
             if (idx < 0)
             {
@@ -473,7 +474,7 @@ namespace MathNet.Numerics.Distributions
             }
 
             var cdfUnnormalized = ProbabilityMassToCumulativeDistribution(probabilityMass);
-            var denormalizedProbability = probability * cdfUnnormalized[cdfUnnormalized.Length - 1];
+            var denormalizedProbability = probability*cdfUnnormalized[cdfUnnormalized.Length - 1];
             int idx = Array.BinarySearch(cdfUnnormalized, denormalizedProbability);
             if (idx < 0)
             {
@@ -502,7 +503,7 @@ namespace MathNet.Numerics.Distributions
                 throw new ArgumentOutOfRangeException("probability");
             }
 
-            var denormalizedProbability = probability * cdfUnnormalized[cdfUnnormalized.Length - 1];
+            var denormalizedProbability = probability*cdfUnnormalized[cdfUnnormalized.Length - 1];
             int idx = Array.BinarySearch(cdfUnnormalized, denormalizedProbability);
             if (idx < 0)
             {
@@ -551,6 +552,34 @@ namespace MathNet.Numerics.Distributions
             return idx;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double[] cdfUnnormalized)
+        {
+            // TODO : use binary search to speed up this procedure.
+            var uniform = rnd.NextDoubles(values.Length);
+            double w = cdfUnnormalized[cdfUnnormalized.Length - 1];
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    var u = uniform[i]*w;
+                    var idx = 0;
+                    while (u > cdfUnnormalized[idx])
+                    {
+                        idx++;
+                    }
+                    values[i] = idx;
+                }
+            });
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double[] cdfUnnormalized)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, cdfUnnormalized);
+            }
+        }
+
         /// <summary>
         /// Samples a Binomially distributed random variable.
         /// </summary>
@@ -560,16 +589,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _cdfUnnormalized);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _cdfUnnormalized);
+        }
+
         /// <summary>
         /// Samples an array of Bernoulli distributed random variables.
         /// </summary>
         /// <returns>a sequence of successful trial counts.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _cdfUnnormalized);
-            }
+            return SamplesUnchecked(_random, _cdfUnnormalized);
         }
 
         /// <summary>
@@ -581,9 +615,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double[] probabilityMass)
         {
             if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
             var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
             return SampleUnchecked(rnd, cdf);
@@ -598,15 +630,69 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double[] probabilityMass)
         {
             if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
             var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, cdf);
-            }
+            return SamplesUnchecked(rnd, cdf);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="probabilityMass">An array of nonnegative ratios. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static void Samples(System.Random rnd, int[] values, double[] probabilityMass)
+        {
+            if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
+            SamplesUnchecked(rnd, values, cdf);
+        }
+
+        /// <summary>
+        /// Samples one categorical distributed random variable; also known as the Discrete distribution.
+        /// </summary>
+        /// <param name="probabilityMass">An array of nonnegative ratios. Not assumed to be normalized.</param>
+        /// <returns>One random integer between 0 and the size of the categorical (exclusive).</returns>
+        public static int Sample(double[] probabilityMass)
+        {
+            if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
+            return SampleUnchecked(SystemRandomSource.Default, cdf);
+        }
+
+        /// <summary>
+        /// Samples a categorically distributed random variable.
+        /// </summary>
+        /// <param name="probabilityMass">An array of nonnegative ratios. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static IEnumerable<int> Samples(double[] probabilityMass)
+        {
+            if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
+            return SamplesUnchecked(SystemRandomSource.Default, cdf);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="probabilityMass">An array of nonnegative ratios. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static void Samples(int[] values, double[] probabilityMass)
+        {
+            if (Control.CheckDistributionParameters && !IsValidProbabilityMass(probabilityMass))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var cdf = ProbabilityMassToCumulativeDistribution(probabilityMass);
+            SamplesUnchecked(SystemRandomSource.Default, values, cdf);
         }
 
         /// <summary>
@@ -618,9 +704,7 @@ namespace MathNet.Numerics.Distributions
         public static int SampleWithCumulativeDistribution(System.Random rnd, double[] cdfUnnormalized)
         {
             if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
             return SampleUnchecked(rnd, cdfUnnormalized);
         }
@@ -634,14 +718,64 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> SamplesWithCumulativeDistribution(System.Random rnd, double[] cdfUnnormalized)
         {
             if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, cdfUnnormalized);
-            }
+            return SamplesUnchecked(rnd, cdfUnnormalized);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="cdfUnnormalized">An array of the cumulative distribution. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static void SamplesWithCumulativeDistribution(System.Random rnd, int[] values, double[] cdfUnnormalized)
+        {
+            if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, cdfUnnormalized);
+        }
+
+        /// <summary>
+        /// Samples one categorical distributed random variable; also known as the Discrete distribution.
+        /// </summary>
+        /// <param name="cdfUnnormalized">An array of the cumulative distribution. Not assumed to be normalized.</param>
+        /// <returns>One random integer between 0 and the size of the categorical (exclusive).</returns>
+        public static int SampleWithCumulativeDistribution(double[] cdfUnnormalized)
+        {
+            if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, cdfUnnormalized);
+        }
+
+        /// <summary>
+        /// Samples a categorically distributed random variable.
+        /// </summary>
+        /// <param name="cdfUnnormalized">An array of the cumulative distribution. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static IEnumerable<int> SamplesWithCumulativeDistribution(double[] cdfUnnormalized)
+        {
+            if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, cdfUnnormalized);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="cdfUnnormalized">An array of the cumulative distribution. Not assumed to be normalized.</param>
+        /// <returns>random integers between 0 and the size of the categorical (exclusive).</returns>
+        public static void SamplesWithCumulativeDistribution(int[] values, double[] cdfUnnormalized)
+        {
+            if (Control.CheckDistributionParameters && !IsValidCumulativeDistribution(cdfUnnormalized))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, cdfUnnormalized);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Cauchy.cs b/src/Numerics/Distributions/Cauchy.cs
index 38b91037..9b70b9b9 100644
--- a/src/Numerics/Distributions/Cauchy.cs
+++ b/src/Numerics/Distributions/Cauchy.cs
@@ -261,6 +261,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _location, _scale);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _location, _scale);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Cauchy distribution.
         /// </summary>
@@ -275,14 +283,6 @@ namespace MathNet.Numerics.Distributions
             return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
         }
 
-        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
-        {
-            while (true)
-            {
-                yield return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
-            }
-        }
-
         static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale)
         {
             rnd.NextDoubles(values);
@@ -295,6 +295,14 @@ namespace MathNet.Numerics.Distributions
             });
         }
 
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
+        {
+            while (true)
+            {
+                yield return location + scale*Math.Tan(Constants.Pi*(rnd.NextDouble() - 0.5));
+            }
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
diff --git a/src/Numerics/Distributions/Chi.cs b/src/Numerics/Distributions/Chi.cs
index 1b8d04e2..54c3e511 100644
--- a/src/Numerics/Distributions/Chi.cs
+++ b/src/Numerics/Distributions/Chi.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -237,7 +238,15 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return SampleUnchecked(_random, (int) _freedom);
+            return SampleUnchecked(_random, (int)_freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, (int)_freedom);
         }
 
         /// <summary>
@@ -246,11 +255,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            var freedom = (int)_freedom;
-            while (true)
-            {
-                yield return SampleUnchecked(_random, freedom);
-            }
+            return SamplesUnchecked(_random, (int)_freedom);
         }
 
         /// <summary>
@@ -266,10 +271,36 @@ namespace MathNet.Numerics.Distributions
             {
                 sum += Math.Pow(Normal.Sample(rnd, 0.0, 1.0), 2);
             }
-
             return Math.Sqrt(sum);
         }
 
+        static void SamplesUnchecked(System.Random rnd, double[] values, int freedom)
+        {
+            var standard = new double[values.Length*freedom];
+            Normal.SamplesUnchecked(rnd, standard, 0.0, 1.0);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    int k = i*freedom;
+                    double sum = 0;
+                    for (int j = 0; j < freedom; j++)
+                    {
+                        sum += standard[k + j]*standard[k + j];
+                    }
+                    values[i] = Math.Sqrt(sum);
+                }
+            });
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, int freedom)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, freedom);
+            }
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
@@ -335,10 +366,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (freedom <= 0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, freedom);
-            }
+            return SamplesUnchecked(rnd, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, int freedom)
+        {
+            if (freedom <= 0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, freedom);
         }
 
         /// <summary>
@@ -348,7 +390,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(int freedom)
         {
-            return Sample(SystemRandomSource.Default, freedom);
+            if (freedom <= 0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, freedom);
         }
 
         /// <summary>
@@ -358,7 +402,22 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(int freedom)
         {
-            return Samples(SystemRandomSource.Default, freedom);
+            if (freedom <= 0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, int freedom)
+        {
+            if (freedom <= 0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, freedom);
         }
     }
 }
diff --git a/src/Numerics/Distributions/ChiSquared.cs b/src/Numerics/Distributions/ChiSquared.cs
index 16c4214a..ffe7c1e9 100644
--- a/src/Numerics/Distributions/ChiSquared.cs
+++ b/src/Numerics/Distributions/ChiSquared.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -133,7 +134,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double StdDev
         {
-            get { return Math.Sqrt(2.0 * _freedom); }
+            get { return Math.Sqrt(2.0*_freedom); }
         }
 
         /// <summary>
@@ -149,7 +150,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Skewness
         {
-            get { return Math.Sqrt(8.0 / _freedom); }
+            get { return Math.Sqrt(8.0/_freedom); }
         }
 
         /// <summary>
@@ -165,7 +166,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public double Median
         {
-            get { return _freedom - (2.0 / 3.0); }
+            get { return _freedom - (2.0/3.0); }
         }
 
         /// <summary>
@@ -226,16 +227,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _freedom);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _freedom);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the <c>ChiSquare</c> distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _freedom);
-            }
+            return SamplesUnchecked(_random, _freedom);
         }
 
         /// <summary>
@@ -246,7 +252,7 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a random number from the distribution.</returns>
         static double SampleUnchecked(System.Random rnd, double freedom)
         {
-            // Use the simple method if the degrees if freedom is an integer anyway
+            // Use the simple method if the degrees of freedom is an integer anyway
             if (Math.Floor(freedom) == freedom && freedom < Int32.MaxValue)
             {
                 double sum = 0;
@@ -258,9 +264,46 @@ namespace MathNet.Numerics.Distributions
                 return sum;
             }
 
-            //Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations
-            //for a justification)
-            return Gamma.SampleUnchecked(rnd, freedom / 2.0, .5);
+            // Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations
+            // for a justification)
+            return Gamma.SampleUnchecked(rnd, freedom/2.0, .5);
+        }
+
+        internal static void SamplesUnchecked(System.Random rnd, double[] values, double freedom)
+        {
+            // Use the simple method if the degrees of freedom is an integer anyway
+            if (Math.Floor(freedom) == freedom && freedom < Int32.MaxValue)
+            {
+                var n = (int)freedom;
+                var standard = new double[values.Length*n];
+                Normal.SamplesUnchecked(rnd, standard, 0.0, 1.0);
+                CommonParallel.For(0, values.Length, 4096, (a, b) =>
+                {
+                    for (int i = a; i < b; i++)
+                    {
+                        int k = i*n;
+                        double sum = 0;
+                        for (int j = 0; j < n; j++)
+                        {
+                            sum += standard[k + j]*standard[k + j];
+                        }
+                        values[i] = Math.Sqrt(sum);
+                    }
+                });
+                return;
+            }
+
+            // Call the gamma function (see http://en.wikipedia.org/wiki/Gamma_distribution#Specializations
+            // for a justification)
+            Gamma.SamplesUnchecked(rnd, values, freedom/2.0, .5);
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double freedom)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, freedom);
+            }
         }
 
         /// <summary>
@@ -328,10 +371,21 @@ namespace MathNet.Numerics.Distributions
         {
             if (freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, freedom);
-            }
+            return SamplesUnchecked(rnd, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
+        /// <returns>a sample from the distribution. </returns>
+        public static void Samples(System.Random rnd, double[] values, double freedom)
+        {
+            if (freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, freedom);
         }
 
         /// <summary>
@@ -341,7 +395,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution. </returns>
         public static double Sample(double freedom)
         {
-            return Sample(SystemRandomSource.Default, freedom);
+            if (freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, freedom);
         }
 
         /// <summary>
@@ -351,7 +407,22 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution. </returns>
         public static IEnumerable<double> Samples(double freedom)
         {
-            return Samples(SystemRandomSource.Default, freedom);
+            if (freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="freedom">The degrees of freedom (k) of the distribution. Range: k > 0.</param>
+        /// <returns>a sample from the distribution. </returns>
+        public static void Samples(double[] values, double freedom)
+        {
+            if (freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, freedom);
         }
     }
 }
diff --git a/src/Numerics/Distributions/ContinuousUniform.cs b/src/Numerics/Distributions/ContinuousUniform.cs
index b4e0363a..877e51b0 100644
--- a/src/Numerics/Distributions/ContinuousUniform.cs
+++ b/src/Numerics/Distributions/ContinuousUniform.cs
@@ -265,6 +265,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _lower, _upper);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _lower, _upper);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the <c>ContinuousUniform</c> distribution.
         /// </summary>
@@ -288,7 +296,7 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
-        static void SamplesUnchecked(System.Random rnd, double[] values, double lower, double upper)
+        internal static void SamplesUnchecked(System.Random rnd, double[] values, double lower, double upper)
         {
             rnd.NextDoubles(values);
             var difference = upper - lower;
diff --git a/src/Numerics/Distributions/ConwayMaxwellPoisson.cs b/src/Numerics/Distributions/ConwayMaxwellPoisson.cs
index 5b92cf53..01284511 100644
--- a/src/Numerics/Distributions/ConwayMaxwellPoisson.cs
+++ b/src/Numerics/Distributions/ConwayMaxwellPoisson.cs
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -505,6 +506,36 @@ namespace MathNet.Numerics.Distributions
             return i;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double lambda, double nu, double z)
+        {
+            var uniform = rnd.NextDoubles(values.Length);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    var u = uniform[i];
+                    var p = 1.0/z;
+                    var cdf = p;
+                    var k = 0;
+                    while (u > cdf)
+                    {
+                        k++;
+                        p = p*lambda/Math.Pow(k, nu);
+                        cdf += p;
+                    }
+                    values[i] = k;
+                }
+            });
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double lambda, double nu, double z)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, lambda, nu, z);
+            }
+        }
+
         /// <summary>
         /// Samples a Conway-Maxwell-Poisson distributed random variable.
         /// </summary>
@@ -514,6 +545,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _lambda, _nu, Z);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _lambda, _nu, Z);
+        }
+
         /// <summary>
         /// Samples a sequence of a Conway-Maxwell-Poisson distributed random variables.
         /// </summary>
@@ -522,10 +561,7 @@ namespace MathNet.Numerics.Distributions
         /// </returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _lambda, _nu, Z);
-            }
+            return SamplesUnchecked(_random, _lambda, _nu, Z);
         }
 
         /// <summary>
@@ -537,6 +573,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double lambda, double nu)
         {
             if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             var z = Normalization(lambda, nu);
             return SampleUnchecked(rnd, lambda, nu, z);
         }
@@ -550,11 +587,24 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double lambda, double nu)
         {
             if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             var z = Normalization(lambda, nu);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lambda, nu, z);
-            }
+            return SamplesUnchecked(rnd, lambda, nu, z);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lambda">The lambda (λ) parameter. Range: λ > 0.</param>
+        /// <param name="nu">The rate of decay (ν) parameter. Range: ν ≥ 0.</param>
+        public static void Samples(System.Random rnd, int[] values, double lambda, double nu)
+        {
+            if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var z = Normalization(lambda, nu);
+            SamplesUnchecked(rnd, values, lambda, nu, z);
         }
 
         /// <summary>
@@ -565,6 +615,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double lambda, double nu)
         {
             if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             var z = Normalization(lambda, nu);
             return SampleUnchecked(SystemRandomSource.Default, lambda, nu, z);
         }
@@ -577,12 +628,23 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double lambda, double nu)
         {
             if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            var z = Normalization(lambda, nu);
+            return SamplesUnchecked(SystemRandomSource.Default, lambda, nu, z);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lambda">The lambda (λ) parameter. Range: λ > 0.</param>
+        /// <param name="nu">The rate of decay (ν) parameter. Range: ν ≥ 0.</param>
+        public static void Samples(int[] values, double lambda, double nu)
+        {
+            if (!(lambda > 0.0 && nu >= 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             var z = Normalization(lambda, nu);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lambda, nu, z);
-            }
+            SamplesUnchecked(SystemRandomSource.Default, values, lambda, nu, z);
         }
     }
 }
diff --git a/src/Numerics/Distributions/DiscreteUniform.cs b/src/Numerics/Distributions/DiscreteUniform.cs
index cf6479e2..8efbbd62 100644
--- a/src/Numerics/Distributions/DiscreteUniform.cs
+++ b/src/Numerics/Distributions/DiscreteUniform.cs
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -51,8 +52,8 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// Initializes a new instance of the DiscreteUniform class.
         /// </summary>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         public DiscreteUniform(int lower, int upper)
         {
             if (!IsValidParameterSet(lower, upper))
@@ -68,8 +69,8 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// Initializes a new instance of the DiscreteUniform class.
         /// </summary>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <param name="randomSource">The random number generator which is used to draw random samples.</param>
         public DiscreteUniform(int lower, int upper, System.Random randomSource)
         {
@@ -97,15 +98,15 @@ namespace MathNet.Numerics.Distributions
         /// <summary>
         /// Tests whether the provided values are valid parameters for this distribution.
         /// </summary>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         public static bool IsValidParameterSet(int lower, int upper)
         {
             return lower <= upper;
         }
 
         /// <summary>
-        /// Gets or sets the lower bound of the probability distribution.
+        /// Gets the inclusive lower bound of the probability distribution.
         /// </summary>
         public int LowerBound
         {
@@ -113,7 +114,7 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Gets or sets the upper bound of the probability distribution.
+        /// Gets the inclusive upper bound of the probability distribution.
         /// </summary>
         public int UpperBound
         {
@@ -235,8 +236,8 @@ namespace MathNet.Numerics.Distributions
         /// Computes the probability mass (PMF) at k, i.e. P(X = k).
         /// </summary>
         /// <param name="k">The location in the domain where we want to evaluate the probability mass function.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>the probability mass at location <paramref name="k"/>.</returns>
         public static double PMF(int lower, int upper, int k)
         {
@@ -249,8 +250,8 @@ namespace MathNet.Numerics.Distributions
         /// Computes the log probability mass (lnPMF) at k, i.e. ln(P(X = k)).
         /// </summary>
         /// <param name="k">The location in the domain where we want to evaluate the log probability mass function.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>the log probability mass at location <paramref name="k"/>.</returns>
         public static double PMFLn(int lower, int upper, int k)
         {
@@ -263,8 +264,8 @@ namespace MathNet.Numerics.Distributions
         /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
         /// </summary>
         /// <param name="x">The location at which to compute the cumulative distribution function.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
         /// <seealso cref="CumulativeDistribution"/>
         public static double CDF(int lower, int upper, double x)
@@ -281,12 +282,22 @@ namespace MathNet.Numerics.Distributions
         /// Generates one sample from the discrete uniform distribution. This method does not do any parameter checking.
         /// </summary>
         /// <param name="rnd">The random source to use.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>A random sample from the discrete uniform distribution.</returns>
         static int SampleUnchecked(System.Random rnd, int lower, int upper)
         {
-            return (rnd.Next()%(upper - lower + 1)) + lower;
+            return rnd.Next(lower, upper + 1);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, int[] values, int lower, int upper)
+        {
+            rnd.NextInt32s(values, lower, upper + 1);
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, int lower, int upper)
+        {
+            return rnd.NextInt32Sequence(lower, upper + 1);
         }
 
         /// <summary>
@@ -298,28 +309,34 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _lower, _upper);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _lower, _upper);
+        }
+
         /// <summary>
         /// Samples an array of uniformly distributed random variables.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _lower, _upper);
-            }
+            return SamplesUnchecked(_random, _lower, _upper);
         }
 
         /// <summary>
         /// Samples a uniformly distributed random variable.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>A sample from the discrete uniform distribution.</returns>
         public static int Sample(System.Random rnd, int lower, int upper)
         {
             if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, lower, upper);
         }
 
@@ -327,44 +344,69 @@ namespace MathNet.Numerics.Distributions
         /// Samples a sequence of uniformly distributed random variables.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>a sequence of samples from the discrete uniform distribution.</returns>
         public static IEnumerable<int> Samples(System.Random rnd, int lower, int upper)
         {
             if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lower, upper);
-            }
+
+            return SamplesUnchecked(rnd, lower, upper);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
+        /// <returns>a sequence of samples from the discrete uniform distribution.</returns>
+        public static void Samples(System.Random rnd, int[] values, int lower, int upper)
+        {
+            if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, lower, upper);
         }
 
         /// <summary>
         /// Samples a uniformly distributed random variable.
         /// </summary>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>A sample from the discrete uniform distribution.</returns>
         public static int Sample(int lower, int upper)
         {
             if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, lower, upper);
         }
 
         /// <summary>
         /// Samples a sequence of uniformly distributed random variables.
         /// </summary>
-        /// <param name="lower">Lower bound. Range: lower ≤ upper.</param>
-        /// <param name="upper">Upper bound. Range: lower ≤ upper.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
         /// <returns>a sequence of samples from the discrete uniform distribution.</returns>
         public static IEnumerable<int> Samples(int lower, int upper)
         {
             if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lower, upper);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, lower, upper);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lower">Lower bound, inclusive. Range: lower ≤ upper.</param>
+        /// <param name="upper">Upper bound, inclusive. Range: lower ≤ upper.</param>
+        /// <returns>a sequence of samples from the discrete uniform distribution.</returns>
+        public static void Samples(int[] values, int lower, int upper)
+        {
+            if (!(lower <= upper)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, lower, upper);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Erlang.cs b/src/Numerics/Distributions/Erlang.cs
index e6e6eae3..b6abee39 100644
--- a/src/Numerics/Distributions/Erlang.cs
+++ b/src/Numerics/Distributions/Erlang.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -45,7 +45,7 @@ namespace MathNet.Numerics.Distributions
     {
         System.Random _random;
 
-        readonly double _shape;
+        readonly int _shape;
         readonly double _rate;
 
         /// <summary>
@@ -121,9 +121,15 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
+        public static bool IsValidParameterSet(int shape, double rate)
+        {
+            return shape >= 0 && rate >= 0.0;
+        }
+
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
         public static bool IsValidParameterSet(double shape, double rate)
         {
-            return shape >= 0.0 && rate >= 0.0;
+            return IsValidParameterSet((int)shape, rate);
         }
 
         /// <summary>
@@ -131,7 +137,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         public int Shape
         {
-            get { return (int)_shape; }
+            get { return _shape; }
         }
 
         /// <summary>
@@ -319,7 +325,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="x">The location at which to compute the density.</param>
         /// <returns>the density at <paramref name="x"/>.</returns>
-        /// <seealso cref="PDF"/>
+        /// <seealso cref="PDF(int, double, double)"/>
         public double Density(double x)
         {
             return PDF(_shape, _rate, x);
@@ -330,7 +336,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="x">The location at which to compute the log density.</param>
         /// <returns>the log density at <paramref name="x"/>.</returns>
-        /// <seealso cref="PDFLn"/>
+        /// <seealso cref="PDFLn(int, double, double)"/>
         public double DensityLn(double x)
         {
             return PDFLn(_shape, _rate, x);
@@ -341,7 +347,7 @@ namespace MathNet.Numerics.Distributions
         /// </summary>
         /// <param name="x">The location at which to compute the cumulative distribution function.</param>
         /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
-        /// <seealso cref="CDF"/>
+        /// <seealso cref="CDF(int, double, double)"/>
         public double CumulativeDistribution(double x)
         {
             return CDF(_shape, _rate, x);
@@ -353,72 +359,26 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public double Sample()
         {
-            return SampleUnchecked(_random, _shape, _rate);
+            return Gamma.SampleUnchecked(_random, _shape, _rate);
         }
 
         /// <summary>
-        /// Generates a sequence of samples from the Erlang distribution.
+        /// Fills an array with samples generated from the distribution.
         /// </summary>
-        /// <returns>a sequence of samples from the distribution.</returns>
-        public IEnumerable<double> Samples()
+        public void Samples(double[] values)
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _shape, _rate);
-            }
+            Gamma.SamplesUnchecked(_random, values, _shape, _rate);
         }
 
         /// <summary>
-        /// <para>Sampling implementation based on:
-        /// "A Simple Method for Generating Erlang Variables" - Marsaglia &amp; Tsang
-        /// ACM Transactions on Mathematical Software, Vol. 26, No. 3, September 2000, Pages 363–372.</para>
-        /// <para>This method performs no parameter checks.</para>
+        /// Generates a sequence of samples from the Erlang distribution.
         /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
-        /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
-        /// <returns>A sample from a Erlang distributed random variable.</returns>
-        static double SampleUnchecked(System.Random rnd, double shape, double rate)
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public IEnumerable<double> Samples()
         {
-            if (Double.IsPositiveInfinity(rate))
-            {
-                return shape;
-            }
-
-            var a = shape;
-            var alphafix = 1.0;
-
-            // Fix when alpha is less than one.
-            if (shape < 1.0)
-            {
-                a = shape + 1.0;
-                alphafix = Math.Pow(rnd.NextDouble(), 1.0 / shape);
-            }
-
-            var d = a - (1.0 / 3.0);
-            var c = 1.0 / Math.Sqrt(9.0 * d);
             while (true)
             {
-                var x = Normal.Sample(rnd, 0.0, 1.0);
-                var v = 1.0 + (c * x);
-                while (v <= 0.0)
-                {
-                    x = Normal.Sample(rnd, 0.0, 1.0);
-                    v = 1.0 + (c * x);
-                }
-
-                v = v * v * v;
-                var u = rnd.NextDouble();
-                x = x * x;
-                if (u < 1.0 - (0.0331 * x * x))
-                {
-                    return alphafix * d * v / rate;
-                }
-
-                if (Math.Log(u) < (0.5 * x) + (d * (1.0 - v + Math.Log(v))))
-                {
-                    return alphafix * d * v / rate;
-                }
+                yield return Gamma.SampleUnchecked(_random, _shape, _rate);
             }
         }
 
@@ -430,7 +390,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="x">The location at which to compute the density.</param>
         /// <returns>the density at <paramref name="x"/>.</returns>
         /// <seealso cref="Density"/>
-        public static double PDF(double shape, double rate, double x)
+        public static double PDF(int shape, double rate, double x)
         {
             if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
@@ -441,6 +401,12 @@ namespace MathNet.Numerics.Distributions
             return Math.Pow(rate, shape)*Math.Pow(x, shape - 1.0)*Math.Exp(-rate*x)/SpecialFunctions.Gamma(shape);
         }
 
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
+        public static double PDF(double shape, double rate, double x)
+        {
+            return PDF((int)shape, rate, x);
+        }
+
         /// <summary>
         /// Computes the log probability density of the distribution (lnPDF) at x, i.e. ln(∂P(X ≤ x)/∂x).
         /// </summary>
@@ -449,7 +415,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="x">The location at which to compute the density.</param>
         /// <returns>the log density at <paramref name="x"/>.</returns>
         /// <seealso cref="DensityLn"/>
-        public static double PDFLn(double shape, double rate, double x)
+        public static double PDFLn(int shape, double rate, double x)
         {
             if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
@@ -460,6 +426,12 @@ namespace MathNet.Numerics.Distributions
             return (shape*Math.Log(rate)) + ((shape - 1.0)*Math.Log(x)) - (rate*x) - SpecialFunctions.GammaLn(shape);
         }
 
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
+        public static double PDFLn(double shape, double rate, double x)
+        {
+            return PDFLn((int)shape, rate, x);
+        }
+
         /// <summary>
         /// Computes the cumulative distribution (CDF) of the distribution at x, i.e. P(X ≤ x).
         /// </summary>
@@ -468,7 +440,7 @@ namespace MathNet.Numerics.Distributions
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
         /// <returns>the cumulative distribution at location <paramref name="x"/>.</returns>
         /// <seealso cref="CumulativeDistribution"/>
-        public static double CDF(double shape, double rate, double x)
+        public static double CDF(int shape, double rate, double x)
         {
             if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
@@ -478,6 +450,12 @@ namespace MathNet.Numerics.Distributions
             return SpecialFunctions.GammaLowerRegularized(shape, x*rate);
         }
 
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
+        public static double CDF(double shape, double rate, double x)
+        {
+            return CDF((int)shape, rate, x);
+        }
+
         /// <summary>
         /// Generates a sample from the distribution.
         /// </summary>
@@ -485,11 +463,15 @@ namespace MathNet.Numerics.Distributions
         /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
         /// <returns>a sample from the distribution.</returns>
-        public static double Sample(System.Random rnd, double shape, double rate)
+        public static double Sample(System.Random rnd, int shape, double rate)
         {
-            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+            return Gamma.Sample(rnd, shape, rate);
+        }
 
-            return SampleUnchecked(rnd, shape, rate);
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
+        public static double Sample(System.Random rnd, double shape, double rate)
+        {
+            return Sample(rnd, (int)shape, rate);
         }
 
         /// <summary>
@@ -499,14 +481,28 @@ namespace MathNet.Numerics.Distributions
         /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
+        public static IEnumerable<double> Samples(System.Random rnd, int shape, double rate)
+        {
+            return Gamma.Samples(rnd, shape, rate);
+        }
+
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
         public static IEnumerable<double> Samples(System.Random rnd, double shape, double rate)
         {
-            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+            return Samples(rnd, (int)shape, rate);
+        }
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, shape, rate);
-            }
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
+        /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, int shape, double rate)
+        {
+            Gamma.Samples(rnd, values, shape, rate);
         }
 
         /// <summary>
@@ -515,9 +511,15 @@ namespace MathNet.Numerics.Distributions
         /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
         /// <returns>a sample from the distribution.</returns>
+        public static double Sample(int shape, double rate)
+        {
+            return Gamma.Sample(shape, rate);
+        }
+
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
         public static double Sample(double shape, double rate)
         {
-            return Sample(SystemRandomSource.Default, shape, rate);
+            return Sample((int)shape, rate);
         }
 
         /// <summary>
@@ -526,9 +528,27 @@ namespace MathNet.Numerics.Distributions
         /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
         /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
+        public static IEnumerable<double> Samples(int shape, double rate)
+        {
+            return Gamma.Samples(shape, rate);
+        }
+
+        [Obsolete("Use the variant that expects an int shape, or use Gamma instead. Will be dropped in v4.")]
         public static IEnumerable<double> Samples(double shape, double rate)
         {
-            return Samples(SystemRandomSource.Default, shape, rate);
+            return Samples((int)shape, rate);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (k) of the Erlang distribution. Range: k ≥ 0.</param>
+        /// <param name="rate">The rate or inverse scale (λ) of the Erlang distribution. Range: λ ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, int shape, double rate)
+        {
+            Gamma.Samples(values, shape, rate);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Exponential.cs b/src/Numerics/Distributions/Exponential.cs
index 4f210c5d..62a58d20 100644
--- a/src/Numerics/Distributions/Exponential.cs
+++ b/src/Numerics/Distributions/Exponential.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -240,6 +240,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _rate);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _rate);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Exponential distribution.
         /// </summary>
@@ -249,12 +257,6 @@ namespace MathNet.Numerics.Distributions
             return SamplesUnchecked(_random, _rate);
         }
 
-        /// <summary>
-        /// Samples the distribution.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
-        /// <returns>a random number from the distribution.</returns>
         static double SampleUnchecked(System.Random rnd, double rate)
         {
             var r = rnd.NextDouble();
@@ -266,12 +268,7 @@ namespace MathNet.Numerics.Distributions
             return -Math.Log(r)/rate;
         }
 
-        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double rate)
-        {
-            return rnd.NextDoubleSequence().Where(r => r != 0.0).Select(r => -Math.Log(r)/rate);
-        }
-
-        static void SamplesUnchecked(System.Random rnd, double[] values, double rate)
+        internal static void SamplesUnchecked(System.Random rnd, double[] values, double rate)
         {
             rnd.NextDoubles(values);
             CommonParallel.For(0, values.Length, 4096, (a, b) =>
@@ -289,6 +286,11 @@ namespace MathNet.Numerics.Distributions
             });
         }
 
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double rate)
+        {
+            return rnd.NextDoubleSequence().Where(r => r != 0.0).Select(r => -Math.Log(r)/rate);
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
@@ -360,30 +362,30 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Generates a sequence of samples from the Exponential distribution.
+        /// Fills an array with samples generated from the distribution.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
         /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
-        public static IEnumerable<double> Samples(System.Random rnd, double rate)
+        public static void Samples(System.Random rnd, double[] values, double rate)
         {
             if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return SamplesUnchecked(rnd, rate);
+            SamplesUnchecked(rnd, values, rate);
         }
 
         /// <summary>
-        /// Fills an array with samples generated from the distribution.
+        /// Generates a sequence of samples from the Exponential distribution.
         /// </summary>
         /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="values">The array to fill with the samples.</param>
         /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
-        public static void Samples(System.Random rnd, double[] values, double rate)
+        public static IEnumerable<double> Samples(System.Random rnd, double rate)
         {
             if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            SamplesUnchecked(rnd, values, rate);
+            return SamplesUnchecked(rnd, rate);
         }
 
         /// <summary>
@@ -399,28 +401,28 @@ namespace MathNet.Numerics.Distributions
         }
 
         /// <summary>
-        /// Generates a sequence of samples from the Exponential distribution.
+        /// Fills an array with samples generated from the distribution.
         /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
         /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
-        public static IEnumerable<double> Samples(double rate)
+        public static void Samples(double[] values, double rate)
         {
             if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return SamplesUnchecked(SystemRandomSource.Default, rate);
+            SamplesUnchecked(SystemRandomSource.Default, values, rate);
         }
 
         /// <summary>
-        /// Fills an array with samples generated from the distribution.
+        /// Generates a sequence of samples from the Exponential distribution.
         /// </summary>
-        /// <param name="values">The array to fill with the samples.</param>
         /// <param name="rate">The rate (λ) parameter of the distribution. Range: λ ≥ 0.</param>
         /// <returns>a sequence of samples from the distribution.</returns>
-        public static void Samples(double[] values, double rate)
+        public static IEnumerable<double> Samples(double rate)
         {
             if (rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            SamplesUnchecked(SystemRandomSource.Default, values, rate);
+            return SamplesUnchecked(SystemRandomSource.Default, rate);
         }
     }
 }
diff --git a/src/Numerics/Distributions/FisherSnedecor.cs b/src/Numerics/Distributions/FisherSnedecor.cs
index 0ad7af65..db074e6f 100644
--- a/src/Numerics/Distributions/FisherSnedecor.cs
+++ b/src/Numerics/Distributions/FisherSnedecor.cs
@@ -33,6 +33,7 @@ using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
 using MathNet.Numerics.RootFinding;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -263,7 +264,7 @@ namespace MathNet.Numerics.Distributions
         {
             return SpecialFunctions.BetaRegularized(_freedom1/2.0, _freedom2/2.0, _freedom1*x/((_freedom1*x) + _freedom2));
         }
-        
+
         /// <summary>
         /// Computes the inverse of the cumulative distribution function (InvCDF) for the distribution
         /// at the given probability. This is also known as the quantile or percent point function.
@@ -286,16 +287,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _freedom1, _freedom2);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _freedom1, _freedom2);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the <c>FisherSnedecor</c> distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _freedom1, _freedom2);
-            }
+            return SamplesUnchecked(_random, _freedom1, _freedom2);
         }
 
         /// <summary>
@@ -307,7 +313,29 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a <c>FisherSnedecor</c> distributed random number.</returns>
         static double SampleUnchecked(System.Random rnd, double d1, double d2)
         {
-            return (ChiSquared.Sample(rnd, d1) / d1) / (ChiSquared.Sample(rnd, d2) / d2);
+            return (ChiSquared.Sample(rnd, d1)*d2)/(ChiSquared.Sample(rnd, d2)*d1);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double d1, double d2)
+        {
+            var values2 = new double[values.Length];
+            ChiSquared.SamplesUnchecked(rnd, values, d1);
+            ChiSquared.SamplesUnchecked(rnd, values2, d2);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = (values[i]*d2)/(values2[i]*d1);
+                }
+            });
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double d1, double d2)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, d1, d2);
+            }
         }
 
         /// <summary>
@@ -397,10 +425,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, d1, d2);
-            }
+            return SamplesUnchecked(rnd, d1, d2);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="d1">The first degree of freedom (d1) of the distribution. Range: d1 > 0.</param>
+        /// <param name="d2">The second degree of freedom (d2) of the distribution. Range: d2 > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double d1, double d2)
+        {
+            if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, d1, d2);
         }
 
         /// <summary>
@@ -411,7 +451,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double d1, double d2)
         {
-            return Sample(SystemRandomSource.Default, d1, d2);
+            if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, d1, d2);
         }
 
         /// <summary>
@@ -422,7 +464,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double d1, double d2)
         {
-            return Samples(SystemRandomSource.Default, d1, d2);
+            if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, d1, d2);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="d1">The first degree of freedom (d1) of the distribution. Range: d1 > 0.</param>
+        /// <param name="d2">The second degree of freedom (d2) of the distribution. Range: d2 > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double d1, double d2)
+        {
+            if (d1 <= 0.0 || d2 <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, d1, d2);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Gamma.cs b/src/Numerics/Distributions/Gamma.cs
index bd0f0da1..4a02265c 100644
--- a/src/Numerics/Distributions/Gamma.cs
+++ b/src/Numerics/Distributions/Gamma.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -372,16 +372,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _shape, _rate);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _shape, _rate);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Gamma distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _shape, _rate);
-            }
+            return SamplesUnchecked(_random, _shape, _rate);
         }
 
         /// <summary>
@@ -438,6 +443,22 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        internal static void SamplesUnchecked(System.Random rnd, double[] values, double shape, double rate)
+        {
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = SampleUnchecked(rnd, shape, rate);
+            }
+        }
+
+        internal static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, location, scale);
+            }
+        }
+
         /// <summary>
         /// Computes the probability density of the distribution (PDF) at x, i.e. ∂P(X ≤ x)/∂x.
         /// </summary>
@@ -535,10 +556,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, shape, rate);
-            }
+            return SamplesUnchecked(rnd, shape, rate);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (k, α) of the Gamma distribution. Range: α ≥ 0.</param>
+        /// <param name="rate">The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double shape, double rate)
+        {
+            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, shape, rate);
         }
 
         /// <summary>
@@ -549,7 +582,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double shape, double rate)
         {
-            return Sample(SystemRandomSource.Default, shape, rate);
+            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, shape, rate);
         }
 
         /// <summary>
@@ -560,7 +595,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double shape, double rate)
         {
-            return Samples(SystemRandomSource.Default, shape, rate);
+            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, shape, rate);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (k, α) of the Gamma distribution. Range: α ≥ 0.</param>
+        /// <param name="rate">The rate or inverse scale (β) of the Gamma distribution. Range: β ≥ 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double shape, double rate)
+        {
+            if (shape < 0.0 || rate < 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, shape, rate);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Geometric.cs b/src/Numerics/Distributions/Geometric.cs
index c68afe39..6a43e236 100644
--- a/src/Numerics/Distributions/Geometric.cs
+++ b/src/Numerics/Distributions/Geometric.cs
@@ -30,8 +30,10 @@
 
 using System;
 using System.Collections.Generic;
+using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -268,6 +270,42 @@ namespace MathNet.Numerics.Distributions
             return p == 1.0 ? 1 : (int)Math.Ceiling(Math.Log(1.0 - rnd.NextDouble(), 1.0 - p));
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double p)
+        {
+            if (p == 1.0)
+            {
+                CommonParallel.For(0, values.Length, 4096, (a, b) =>
+                {
+                    for (int i = a; i < b; i++)
+                    {
+                        values[i] = 1;
+                    }
+                });
+                return;
+            }
+
+            var uniform = rnd.NextDoubles(values.Length);
+            double rp = 1.0 - p;
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = (int)Math.Ceiling(Math.Log(1.0 - uniform[i], rp));
+                }
+            });
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double p)
+        {
+            if (p == 1.0)
+            {
+                return Generate.RepeatSequence(1);
+            }
+
+            double rp = 1.0 - p;
+            return rnd.NextDoubleSequence().Select(r => (int)Math.Ceiling(Math.Log(1.0 - r, rp)));
+        }
+
         /// <summary>
         /// Samples a Geometric distributed random variable.
         /// </summary>
@@ -277,16 +315,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _p);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _p);
+        }
+
         /// <summary>
         /// Samples an array of Geometric distributed random variables.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _p);
-            }
+            return SamplesUnchecked(_random, _p);
         }
 
         /// <summary>
@@ -297,6 +340,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, p);
         }
 
@@ -308,10 +352,21 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p);
-            }
+
+            return SamplesUnchecked(rnd, p);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The probability (p) of generating one. Range: 0 ≤ p ≤ 1.</param>
+        public static void Samples(System.Random rnd, int[] values, double p)
+        {
+            if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, p);
         }
 
         /// <summary>
@@ -321,6 +376,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, p);
         }
 
@@ -331,11 +387,20 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double p)
         {
             if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, p);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, p);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="p">The probability (p) of generating one. Range: 0 ≤ p ≤ 1.</param>
+        public static void Samples(int[] values, double p)
+        {
+            if (!(p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, p);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Hypergeometric.cs b/src/Numerics/Distributions/Hypergeometric.cs
index 9897063d..bb37e41a 100644
--- a/src/Numerics/Distributions/Hypergeometric.cs
+++ b/src/Numerics/Distributions/Hypergeometric.cs
@@ -340,6 +340,22 @@ namespace MathNet.Numerics.Distributions
             return x;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, int population, int success, int draws)
+        {
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = SampleUnchecked(rnd, population, success, draws);
+            }
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, int population, int success, int draws)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, population, success, draws);
+            }
+        }
+
         /// <summary>
         /// Samples a Hypergeometric distributed random variable.
         /// </summary>
@@ -349,16 +365,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _population, _success, _draws);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _population, _success, _draws);
+        }
+
         /// <summary>
         /// Samples an array of Hypergeometric distributed random variables.
         /// </summary>
         /// <returns>a sequence of successes in n trials.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _population, _success, _draws);
-            }
+            return SamplesUnchecked(_random, _population, _success, _draws);
         }
 
         /// <summary>
@@ -371,9 +392,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, int population, int success, int draws)
         {
             if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
             return SampleUnchecked(rnd, population, success, draws);
         }
@@ -388,14 +407,25 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, int population, int success, int draws)
         {
             if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, population, success, draws);
-            }
+            return SamplesUnchecked(rnd, population, success, draws);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="population">The size of the population (N).</param>
+        /// <param name="success">The number successes within the population (K, M).</param>
+        /// <param name="draws">The number of draws without replacement (n).</param>
+        public static void Samples(System.Random rnd, int[] values, int population, int success, int draws)
+        {
+            if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, population, success, draws);
         }
 
         /// <summary>
@@ -407,9 +437,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(int population, int success, int draws)
         {
             if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
             return SampleUnchecked(SystemRandomSource.Default, population, success, draws);
         }
@@ -423,15 +451,24 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(int population, int success, int draws)
         {
             if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
-            {
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
-            }
 
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, population, success, draws);
-            }
+            return SamplesUnchecked(SystemRandomSource.Default, population, success, draws);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="population">The size of the population (N).</param>
+        /// <param name="success">The number successes within the population (K, M).</param>
+        /// <param name="draws">The number of draws without replacement (n).</param>
+        public static void Samples(int[] values, int population, int success, int draws)
+        {
+            if (!(population >= 0 && success >= 0 && draws >= 0 && success <= population && draws <= population))
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, population, success, draws);
         }
     }
 }
diff --git a/src/Numerics/Distributions/IContinuousDistribution.cs b/src/Numerics/Distributions/IContinuousDistribution.cs
index 28802418..0bf3e0b5 100644
--- a/src/Numerics/Distributions/IContinuousDistribution.cs
+++ b/src/Numerics/Distributions/IContinuousDistribution.cs
@@ -73,6 +73,11 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         double Sample();
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        void Samples(double[] values);
+
         /// <summary>
         /// Draws a sequence of random samples from the distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/IDiscreteDistribution.cs b/src/Numerics/Distributions/IDiscreteDistribution.cs
index 1949e23c..4f646fb7 100644
--- a/src/Numerics/Distributions/IDiscreteDistribution.cs
+++ b/src/Numerics/Distributions/IDiscreteDistribution.cs
@@ -73,6 +73,11 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         int Sample();
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        void Samples(int[] values);
+
         /// <summary>
         /// Draws a sequence of random samples from the distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/InverseGamma.cs b/src/Numerics/Distributions/InverseGamma.cs
index 774d7eb5..3f151964 100644
--- a/src/Numerics/Distributions/InverseGamma.cs
+++ b/src/Numerics/Distributions/InverseGamma.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -33,6 +33,7 @@ using System.Collections.Generic;
 using System.Linq;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -264,7 +265,15 @@ namespace MathNet.Numerics.Distributions
         /// <returns>A random number from this distribution.</returns>
         public double Sample()
         {
-            return 1.0/Gamma.Sample(_random, _shape, _scale);
+            return SampleUnchecked(_random, _shape, _scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _shape, _scale);
         }
 
         /// <summary>
@@ -273,7 +282,29 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            return Gamma.Samples(_random, _shape, _scale).Select(z => 1.0/z);
+            return SamplesUnchecked(_random, _shape, _scale);
+        }
+
+        static double SampleUnchecked(System.Random rnd, double shape, double scale)
+        {
+            return 1.0/Gamma.SampleUnchecked(rnd, shape, scale);
+        }
+
+        static void SamplesUnchecked(System.Random rnd, double[] values, double shape, double scale)
+        {
+            Gamma.SamplesUnchecked(rnd, values, shape, scale);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    values[i] = 1.0/values[i];
+                }
+            });
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double shape, double scale)
+        {
+            return Gamma.SamplesUnchecked(rnd, shape, scale).Select(z => 1.0/z);
         }
 
         /// <summary>
@@ -330,7 +361,7 @@ namespace MathNet.Numerics.Distributions
         {
             if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return 1.0/Gamma.Sample(rnd, shape, scale);
+            return SampleUnchecked(rnd, shape, scale);
         }
 
         /// <summary>
@@ -344,7 +375,22 @@ namespace MathNet.Numerics.Distributions
         {
             if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            return Gamma.Samples(rnd, shape, scale).Select(z => 1.0/z);
+            return SamplesUnchecked(rnd, shape, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <param name="scale">The scale (β) of the distribution. Range: β > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double shape, double scale)
+        {
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, shape, scale);
         }
 
         /// <summary>
@@ -355,7 +401,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double shape, double scale)
         {
-            return Sample(SystemRandomSource.Default, shape, scale);
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, shape, scale);
         }
 
         /// <summary>
@@ -366,7 +414,23 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double shape, double scale)
         {
-            return Samples(SystemRandomSource.Default, shape, scale);
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, shape, scale);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="shape">The shape (α) of the distribution. Range: α > 0.</param>
+        /// <param name="scale">The scale (β) of the distribution. Range: β > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double shape, double scale)
+        {
+            if (shape <= 0.0 || scale <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, shape, scale);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Laplace.cs b/src/Numerics/Distributions/Laplace.cs
index 43869fe1..bafa61c1 100644
--- a/src/Numerics/Distributions/Laplace.cs
+++ b/src/Numerics/Distributions/Laplace.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -32,6 +32,7 @@ using System;
 using System.Collections.Generic;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics.Distributions
 {
@@ -252,6 +253,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _location, _scale);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _location, _scale);
+        }
+
         /// <summary>
         /// Generates a sample from the Laplace distribution.
         /// </summary>
@@ -261,35 +270,30 @@ namespace MathNet.Numerics.Distributions
             return SamplesUnchecked(_random, _location, _scale);
         }
 
-        /// <summary>
-        /// Samples the distribution.
-        /// </summary>
-        /// <param name="rnd">The random number generator to use.</param>
-        /// <param name="location">The location (μ) of the distribution.</param>
-        /// <param name="scale">The scale (b) of the distribution. Range: b > 0.</param>
-        /// <returns>a random number from the distribution.</returns>
         static double SampleUnchecked(System.Random rnd, double location, double scale)
         {
             var u = rnd.NextDouble() - 0.5;
             return location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
         }
 
-        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
+        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale)
         {
-            while (true)
+            rnd.NextDoubles(values);
+            CommonParallel.For(0, values.Length, 4096, (a, b) =>
             {
-                var u = rnd.NextDouble() - 0.5;
-                yield return location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
-            }
+                for (int i = a; i < b; i++)
+                {
+                    var u = values[i] - 0.5;
+                    values[i] = location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
+                }
+            });
         }
 
-        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale)
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale)
         {
-            rnd.NextDoubles(values);
-            for (int i = 0; i < values.Length; i++)
+            while (true)
             {
-                var u = values[i] - 0.5;
-                values[i] = location - (scale*Math.Sign(u)*Math.Log(1.0 - (2.0*Math.Abs(u))));
+                yield return SampleUnchecked(rnd, location, scale);
             }
         }
 
diff --git a/src/Numerics/Distributions/LogNormal.cs b/src/Numerics/Distributions/LogNormal.cs
index dd8c821a..f10f02b3 100644
--- a/src/Numerics/Distributions/LogNormal.cs
+++ b/src/Numerics/Distributions/LogNormal.cs
@@ -323,6 +323,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _mu, _sigma);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _mu, _sigma);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the log-normal distribution using the <i>Box-Muller</i> algorithm.
         /// </summary>
diff --git a/src/Numerics/Distributions/NegativeBinomial.cs b/src/Numerics/Distributions/NegativeBinomial.cs
index dae87356..c7182366 100644
--- a/src/Numerics/Distributions/NegativeBinomial.cs
+++ b/src/Numerics/Distributions/NegativeBinomial.cs
@@ -298,6 +298,22 @@ namespace MathNet.Numerics.Distributions
             return k - 1;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double r, double p)
+        {
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = SampleUnchecked(rnd, r, p);
+            }
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double r, double p)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, r, p);
+            }
+        }
+
         /// <summary>
         /// Samples a <c>NegativeBinomial</c> distributed random variable.
         /// </summary>
@@ -307,16 +323,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _trials, _p);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _trials, _p);
+        }
+
         /// <summary>
         /// Samples an array of <c>NegativeBinomial</c> distributed random variables.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _trials, _p);
-            }
+            return SamplesUnchecked(_random, _trials, _p);
         }
 
         /// <summary>
@@ -328,6 +349,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double r, double p)
         {
             if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, r, p);
         }
 
@@ -340,10 +362,22 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double r, double p)
         {
             if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, r, p);
-            }
+
+            return SamplesUnchecked(rnd, r, p);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="r">The number of failures (r) until the experiment stopped. Range: r ≥ 0.</param>
+        /// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
+        public static void Samples(System.Random rnd, int[] values, double r, double p)
+        {
+            if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, r, p);
         }
 
         /// <summary>
@@ -354,6 +388,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double r, double p)
         {
             if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, r, p);
         }
 
@@ -365,11 +400,21 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double r, double p)
         {
             if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, r, p);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, r, p);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="r">The number of failures (r) until the experiment stopped. Range: r ≥ 0.</param>
+        /// <param name="p">The probability (p) of a trial resulting in success. Range: 0 ≤ p ≤ 1.</param>
+        public static void Samples(int[] values, double r, double p)
+        {
+            if (!(r >= 0.0 && p >= 0.0 && p <= 1.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, r, p);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Normal.cs b/src/Numerics/Distributions/Normal.cs
index dd1d9c1c..f145a72e 100644
--- a/src/Numerics/Distributions/Normal.cs
+++ b/src/Numerics/Distributions/Normal.cs
@@ -319,6 +319,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _mean, _stdDev);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _mean, _stdDev);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the normal distribution using the <i>Box-Muller</i> algorithm.
         /// </summary>
diff --git a/src/Numerics/Distributions/Pareto.cs b/src/Numerics/Distributions/Pareto.cs
index c368bdb4..44e5e0eb 100644
--- a/src/Numerics/Distributions/Pareto.cs
+++ b/src/Numerics/Distributions/Pareto.cs
@@ -274,6 +274,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _scale, _shape);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _scale, _shape);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Pareto distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/Poisson.cs b/src/Numerics/Distributions/Poisson.cs
index 2775a3e6..85d4097a 100644
--- a/src/Numerics/Distributions/Poisson.cs
+++ b/src/Numerics/Distributions/Poisson.cs
@@ -271,6 +271,72 @@ namespace MathNet.Numerics.Distributions
             return (lambda < 30.0) ? DoSampleShort(rnd, lambda) : DoSampleLarge(rnd, lambda);
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double lambda)
+        {
+            if (lambda < 30.0)
+            {
+                var limit = Math.Exp(-lambda);
+                for (int i = 0; i < values.Length; i++)
+                {
+                    var count = 0;
+                    for (var product = rnd.NextDouble(); product >= limit; product *= rnd.NextDouble())
+                    {
+                        count++;
+                    }
+                    values[i] = count;
+                }
+            }
+            else
+            {
+                var c = 0.767 - (3.36/lambda);
+                var beta = Math.PI/Math.Sqrt(3.0*lambda);
+                var alpha = beta*lambda;
+                var k = Math.Log(c) - lambda - Math.Log(beta);
+                for (int i = 0; i < values.Length; i++)
+                {
+                    for (;;)
+                    {
+                        var u = rnd.NextDouble();
+                        var x = (alpha - Math.Log((1.0 - u)/u))/beta;
+                        var n = (int)Math.Floor(x + 0.5);
+                        if (n < 0)
+                        {
+                            continue;
+                        }
+
+                        var v = rnd.NextDouble();
+                        var y = alpha - (beta*x);
+                        var temp = 1.0 + Math.Exp(y);
+                        var lhs = y + Math.Log(v/(temp*temp));
+                        var rhs = k + (n*Math.Log(lambda)) - SpecialFunctions.FactorialLn(n);
+                        if (lhs <= rhs)
+                        {
+                            values[i] = n;
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double lambda)
+        {
+            if (lambda < 30.0)
+            {
+                while (true)
+                {
+                    yield return DoSampleShort(rnd, lambda);
+                }
+            }
+            else
+            {
+                while (true)
+                {
+                    yield return DoSampleLarge(rnd, lambda);
+                }
+            }
+        }
+
         /// <summary>
         /// Generates one sample from the Poisson distribution by Knuth's method.
         /// </summary>
@@ -336,16 +402,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _lambda);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _lambda);
+        }
+
         /// <summary>
         /// Samples an array of Poisson distributed random variables.
         /// </summary>
         /// <returns>a sequence of successes in N trials.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _lambda);
-            }
+            return SamplesUnchecked(_random, _lambda);
         }
 
         /// <summary>
@@ -357,6 +428,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double lambda)
         {
             if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, lambda);
         }
 
@@ -369,10 +441,22 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double lambda)
         {
             if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lambda);
-            }
+
+            return SamplesUnchecked(rnd, lambda);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lambda">The lambda (λ) parameter of the Poisson distribution. Range: λ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, int[] values, double lambda)
+        {
+            if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, lambda);
         }
 
         /// <summary>
@@ -383,6 +467,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double lambda)
         {
             if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, lambda);
         }
 
@@ -394,11 +479,21 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double lambda)
         {
             if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, lambda);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, lambda);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="lambda">The lambda (λ) parameter of the Poisson distribution. Range: λ > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(int[] values, double lambda)
+        {
+            if (!(lambda > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, lambda);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Rayleigh.cs b/src/Numerics/Distributions/Rayleigh.cs
index 374d256d..0323c6c3 100644
--- a/src/Numerics/Distributions/Rayleigh.cs
+++ b/src/Numerics/Distributions/Rayleigh.cs
@@ -245,6 +245,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _scale);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _scale);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Rayleigh distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/Stable.cs b/src/Numerics/Distributions/Stable.cs
index 27b23a89..7c5b71f6 100644
--- a/src/Numerics/Distributions/Stable.cs
+++ b/src/Numerics/Distributions/Stable.cs
@@ -356,6 +356,53 @@ namespace MathNet.Numerics.Distributions
             }
         }
 
+        static void SamplesUnchecked(System.Random rnd, double[] values, double alpha, double beta, double scale, double location)
+        {
+            var randThetas = new double[values.Length];
+            var randWs = new double[values.Length];
+            ContinuousUniform.SamplesUnchecked(rnd, randThetas, -Constants.PiOver2, Constants.PiOver2);
+            Exponential.SamplesUnchecked(rnd, randWs, 1.0);
+
+            if (!1.0.AlmostEqual(alpha))
+            {
+                for (int i = 0; i < values.Length; i++)
+                {
+                    var randTheta = randThetas[i];
+
+                    var theta = (1.0/alpha)*Math.Atan(beta*Math.Tan(Constants.PiOver2*alpha));
+                    var angle = alpha*(randTheta + theta);
+                    var part1 = beta*Math.Tan(Constants.PiOver2*alpha);
+
+                    var factor = Math.Pow(1.0 + (part1*part1), 1.0/(2.0*alpha));
+                    var factor1 = Math.Sin(angle)/Math.Pow(Math.Cos(randTheta), (1.0/alpha));
+                    var factor2 = Math.Pow(Math.Cos(randTheta - angle)/randWs[i], (1 - alpha)/alpha);
+
+                    values[i] = location + scale*(factor*factor1*factor2);
+                }
+            }
+            else
+            {
+                for (int i = 0; i < values.Length; i++)
+                {
+                    var randTheta = randThetas[i];
+
+                    var part1 = Constants.PiOver2 + (beta*randTheta);
+                    var summand = part1*Math.Tan(randTheta);
+                    var subtrahend = beta*Math.Log(Constants.PiOver2*randWs[i]*Math.Cos(randTheta)/part1);
+
+                    values[i] = location + scale*Constants.TwoInvPi*(summand - subtrahend);
+                }
+            }
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double alpha, double beta, double scale, double location)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, alpha, beta, scale, location);
+            }
+        }
+
         /// <summary>
         /// Draws a random sample from the distribution.
         /// </summary>
@@ -365,16 +412,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _alpha, _beta, _scale, _location);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _alpha, _beta, _scale, _location);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Stable distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _alpha, _beta, _scale, _location);
-            }
+            return SamplesUnchecked(_random, _alpha, _beta, _scale, _location);
         }
 
         /// <summary>
@@ -486,10 +538,25 @@ namespace MathNet.Numerics.Distributions
             if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0)
                 throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, alpha, beta, scale, location);
-            }
+            return SamplesUnchecked(rnd, alpha, beta, scale, location);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="alpha">The stability (α) of the distribution. Range: 2 ≥ α > 0.</param>
+        /// <param name="beta">The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1.</param>
+        /// <param name="scale">The scale (c) of the distribution. Range: c > 0.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double alpha, double beta, double scale, double location)
+        {
+            if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0)
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, alpha, beta, scale, location);
         }
 
         /// <summary>
@@ -502,7 +569,10 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double alpha, double beta, double scale, double location)
         {
-            return Sample(SystemRandomSource.Default, alpha, beta, scale, location);
+            if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0)
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, alpha, beta, scale, location);
         }
 
         /// <summary>
@@ -515,7 +585,27 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double alpha, double beta, double scale, double location)
         {
-            return Samples(SystemRandomSource.Default, alpha, beta, scale, location);
+            if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0)
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, alpha, beta, scale, location);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="alpha">The stability (α) of the distribution. Range: 2 ≥ α > 0.</param>
+        /// <param name="beta">The skewness (β) of the distribution. Range: 1 ≥ β ≥ -1.</param>
+        /// <param name="scale">The scale (c) of the distribution. Range: c > 0.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double alpha, double beta, double scale, double location)
+        {
+            if (alpha <= 0.0 || alpha > 2.0 || beta < -1.0 || beta > 1.0 || scale <= 0.0)
+                throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, alpha, beta, scale, location);
         }
     }
 }
diff --git a/src/Numerics/Distributions/StudentT.cs b/src/Numerics/Distributions/StudentT.cs
index 5991b08c..e16b5d01 100644
--- a/src/Numerics/Distributions/StudentT.cs
+++ b/src/Numerics/Distributions/StudentT.cs
@@ -4,7 +4,7 @@
 // http://github.com/mathnet/mathnet-numerics
 // http://mathnetnumerics.codeplex.com
 //
-// Copyright (c) 2009-2013 Math.NET
+// Copyright (c) 2009-2014 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -339,6 +339,23 @@ namespace MathNet.Numerics.Distributions
             return Normal.Sample(rnd, location, scale*Math.Sqrt(freedom/gamma));
         }
 
+        static void SamplesUnchecked(System.Random rnd, double[] values, double location, double scale, double freedom)
+        {
+            Gamma.SamplesUnchecked(rnd, values, 0.5*freedom, 0.5);
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = Normal.Sample(rnd, location, scale*Math.Sqrt(freedom/values[i]));
+            }
+        }
+
+        static IEnumerable<double> SamplesUnchecked(System.Random rnd, double location, double scale, double freedom)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, location, scale, freedom);
+            }
+        }
+
         /// <summary>
         /// Generates a sample from the Student t-distribution.
         /// </summary>
@@ -348,16 +365,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _location, _scale, _freedom);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _location, _scale, _freedom);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Student t-distribution.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<double> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _location, _scale, _freedom);
-            }
+            return SamplesUnchecked(_random, _location, _scale, _freedom);
         }
 
         /// <summary>
@@ -484,10 +506,23 @@ namespace MathNet.Numerics.Distributions
         {
             if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
 
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, location, scale, freedom);
-            }
+            return SamplesUnchecked(rnd, location, scale, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
+        /// <param name="freedom">The degrees of freedom (ν) for the distribution. Range: ν > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(System.Random rnd, double[] values, double location, double scale, double freedom)
+        {
+            if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, location, scale, freedom);
         }
 
         /// <summary>
@@ -499,7 +534,9 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sample from the distribution.</returns>
         public static double Sample(double location, double scale, double freedom)
         {
-            return Sample(SystemRandomSource.Default, location, scale, freedom);
+            if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SampleUnchecked(SystemRandomSource.Default, location, scale, freedom);
         }
 
         /// <summary>
@@ -511,7 +548,24 @@ namespace MathNet.Numerics.Distributions
         /// <returns>a sequence of samples from the distribution.</returns>
         public static IEnumerable<double> Samples(double location, double scale, double freedom)
         {
-            return Samples(SystemRandomSource.Default, location, scale, freedom);
+            if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            return SamplesUnchecked(SystemRandomSource.Default, location, scale, freedom);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="location">The location (μ) of the distribution.</param>
+        /// <param name="scale">The scale (σ) of the distribution. Range: σ > 0.</param>
+        /// <param name="freedom">The degrees of freedom (ν) for the distribution. Range: ν > 0.</param>
+        /// <returns>a sequence of samples from the distribution.</returns>
+        public static void Samples(double[] values, double location, double scale, double freedom)
+        {
+            if (scale <= 0.0 || freedom <= 0.0) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, location, scale, freedom);
         }
     }
 }
diff --git a/src/Numerics/Distributions/Triangular.cs b/src/Numerics/Distributions/Triangular.cs
index 169fb710..d835f6b4 100644
--- a/src/Numerics/Distributions/Triangular.cs
+++ b/src/Numerics/Distributions/Triangular.cs
@@ -289,6 +289,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _lower, _upper, _mode);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _lower, _upper, _mode);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the <c>Triangular</c> distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/Weibull.cs b/src/Numerics/Distributions/Weibull.cs
index 517df6c9..4ca49099 100644
--- a/src/Numerics/Distributions/Weibull.cs
+++ b/src/Numerics/Distributions/Weibull.cs
@@ -290,6 +290,14 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _shape, _scale);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(double[] values)
+        {
+            SamplesUnchecked(_random, values, _shape, _scale);
+        }
+
         /// <summary>
         /// Generates a sequence of samples from the Weibull distribution.
         /// </summary>
diff --git a/src/Numerics/Distributions/Zipf.cs b/src/Numerics/Distributions/Zipf.cs
index 5a998d46..2949e32e 100644
--- a/src/Numerics/Distributions/Zipf.cs
+++ b/src/Numerics/Distributions/Zipf.cs
@@ -337,6 +337,22 @@ namespace MathNet.Numerics.Distributions
             return i;
         }
 
+        static void SamplesUnchecked(System.Random rnd, int[] values, double s, int n)
+        {
+            for (int i = 0; i < values.Length; i++)
+            {
+                values[i] = SampleUnchecked(rnd, s, n);
+            }
+        }
+
+        static IEnumerable<int> SamplesUnchecked(System.Random rnd, double s, int n)
+        {
+            while (true)
+            {
+                yield return SampleUnchecked(rnd, s, n);
+            }
+        }
+
         /// <summary>
         /// Draws a random sample from the distribution.
         /// </summary>
@@ -346,16 +362,21 @@ namespace MathNet.Numerics.Distributions
             return SampleUnchecked(_random, _s, _n);
         }
 
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        public void Samples(int[] values)
+        {
+            SamplesUnchecked(_random, values, _s, _n);
+        }
+
         /// <summary>
         /// Samples an array of zipf distributed random variables.
         /// </summary>
         /// <returns>a sequence of samples from the distribution.</returns>
         public IEnumerable<int> Samples()
         {
-            while (true)
-            {
-                yield return SampleUnchecked(_random, _s, _n);
-            }
+            return SamplesUnchecked(_random, _s, _n);
         }
 
         /// <summary>
@@ -367,6 +388,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(System.Random rnd, double s, int n)
         {
             if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(rnd, s, n);
         }
 
@@ -379,10 +401,22 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(System.Random rnd, double s, int n)
         {
             if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, s, n);
-            }
+
+            return SamplesUnchecked(rnd, s, n);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="rnd">The random number generator to use.</param>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="s">The s parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        public static void Samples(System.Random rnd, int[] values, double s, int n)
+        {
+            if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(rnd, values, s, n);
         }
 
         /// <summary>
@@ -393,6 +427,7 @@ namespace MathNet.Numerics.Distributions
         public static int Sample(double s, int n)
         {
             if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
             return SampleUnchecked(SystemRandomSource.Default, s, n);
         }
 
@@ -404,11 +439,21 @@ namespace MathNet.Numerics.Distributions
         public static IEnumerable<int> Samples(double s, int n)
         {
             if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
-            SystemRandomSource rnd = SystemRandomSource.Default;
-            while (true)
-            {
-                yield return SampleUnchecked(rnd, s, n);
-            }
+
+            return SamplesUnchecked(SystemRandomSource.Default, s, n);
+        }
+
+        /// <summary>
+        /// Fills an array with samples generated from the distribution.
+        /// </summary>
+        /// <param name="values">The array to fill with the samples.</param>
+        /// <param name="s">The s parameter of the distribution.</param>
+        /// <param name="n">The n parameter of the distribution.</param>
+        public static void Samples(int[] values, double s, int n)
+        {
+            if (!(n > 0 && s > 0.0)) throw new ArgumentException(Resources.InvalidDistributionParameters);
+
+            SamplesUnchecked(SystemRandomSource.Default, values, s, n);
         }
     }
 }
diff --git a/src/Numerics/Generate.cs b/src/Numerics/Generate.cs
index 475528d0..4a70e8ab 100644
--- a/src/Numerics/Generate.cs
+++ b/src/Numerics/Generate.cs
@@ -34,6 +34,7 @@ using System.Linq;
 using MathNet.Numerics.Distributions;
 using MathNet.Numerics.Properties;
 using MathNet.Numerics.Random;
+using MathNet.Numerics.Threading;
 
 namespace MathNet.Numerics
 {
@@ -422,6 +423,36 @@ namespace MathNet.Numerics
             }
         }
 
+        /// <summary>
+        /// Create an array with each field set to the same value.
+        /// </summary>
+        /// <param name="length">The number of samples to generate.</param>
+        /// <param name="value">The value that each field should be set to.</param>
+        public static T[] Repeat<T>(int length, T value)
+        {
+            var data = new T[length];
+            CommonParallel.For(0, data.Length, 4096, (a, b) =>
+            {
+                for (int i = a; i < b; i++)
+                {
+                    data[i] = value;
+                }
+            });
+            return data;
+        }
+
+        /// <summary>
+        /// Create an infinite sequence where each element has the same value.
+        /// </summary>
+        /// <param name="value">The value that each element should be set to.</param>
+        public static IEnumerable<T> RepeatSequence<T>(T value)
+        {
+            while (true)
+            {
+                yield return value;
+            }
+        }
+
         /// <summary>
         /// Create a Heaviside Step sample vector.
         /// </summary>