diff --git a/src/Numerics.Tests/Random/Xoshiro256StarStarTests.cs b/src/Numerics.Tests/Random/Xoshiro256StarStarTests.cs
new file mode 100644
index 00000000..32b3838c
--- /dev/null
+++ b/src/Numerics.Tests/Random/Xoshiro256StarStarTests.cs
@@ -0,0 +1,54 @@
+// <copyright file="Xoshiro256StarStarTests.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// Copyright (c) 2009-2016 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Numerics.Random;
+using NUnit.Framework;
+
+namespace MathNet.Numerics.UnitTests.Random
+{
+    /// <summary>
+    /// Test for xoshiro256** pseudo random number generator.
+    /// </summary>
+    [TestFixture, Category("Random")]
+    public class Xoshiro256StarStarTests : RandomTests
+    {
+        /// <summary>
+        /// Initializes a new instance of the Wh2006Tests class.
+        /// </summary>
+        public Xoshiro256StarStarTests() : base(typeof (Xoshiro256StarStar))
+        {
+        }
+
+        [Test]
+        public void StaticSamplesConsistent()
+        {
+            Assert.That(Xoshiro256StarStar.Doubles(1000, 1), Is.EqualTo(new Xoshiro256StarStar(1).NextDoubles(1000)).Within(1e-12).AsCollection);
+        }
+    }
+}
diff --git a/src/Numerics/Random/Xoshiro256StarStar.cs b/src/Numerics/Random/Xoshiro256StarStar.cs
new file mode 100644
index 00000000..0500e464
--- /dev/null
+++ b/src/Numerics/Random/Xoshiro256StarStar.cs
@@ -0,0 +1,375 @@
+// <copyright file="Xoshiro256StarStar.cs" company="Math.NET">
+// Math.NET Numerics, part of the Math.NET Project
+// http://numerics.mathdotnet.com
+// http://github.com/mathnet/mathnet-numerics
+//
+// 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
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+/*
+   Original code's copyright and license:
+   Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org)
+
+   To the extent possible under law, the author has dedicated all copyright
+   and related and neighboring rights to this software to the public domain
+   worldwide. This software is distributed without any warranty.
+
+   See <http://creativecommons.org/publicdomain/zero/1.0/>.
+*/
+
+using System.Collections.Generic;
+using System.Runtime.Serialization;
+
+#if !NETSTANDARD1_3
+using System;
+using System.Runtime;
+#endif
+
+namespace MathNet.Numerics.Random
+{
+    /// <summary>
+    /// Xoshiro256** pseudo random number generator.
+    /// A random number generator based on the <see cref="System.Random"/> class in the .NET library.
+    /// </summary>
+    /// <remarks>
+    /// This is xoshiro256** 1.0, our all-purpose, rock-solid generator. It has
+    /// excellent(sub-ns) speed, a state space(256 bits) that is large enough
+    /// for any parallel application, and it passes all tests we are aware of.
+    ///
+    /// For generating just floating-point numbers, xoshiro256+ is even faster.
+    ///
+    /// The state must be seeded so that it is not everywhere zero.If you have
+    /// a 64-bit seed, we suggest to seed a splitmix64 generator and use its
+    /// output to fill s.
+    ///
+    /// For further details see:
+    /// David Blackman &amp; Sebastiano Vigna (2018), "Scrambled Linear Pseudorandom Number Generators".
+    /// https://arxiv.org/abs/1805.01407
+    /// </remarks>
+    [Serializable]
+    [DataContract(Namespace = "urn:MathNet/Numerics/Random")]
+    public class Xoshiro256StarStar : RandomSource
+    {
+        // Constants.
+        const double REAL_UNIT_UINT = 1.0 / (1UL << 53);
+
+        // RNG state.
+        [DataMember(Order = 1)]
+        ulong _s0;
+        [DataMember(Order = 2)]
+        ulong _s1;
+        [DataMember(Order = 3)]
+        ulong _s2;
+        [DataMember(Order = 4)]
+        ulong _s3;
+
+        /// <summary>
+        /// Construct a new random number generator with a random seed.
+        /// </summary>
+        public Xoshiro256StarStar() : this(RandomSeed.Robust())
+        {
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Xoshiro256StarStar(bool threadSafe) : this(RandomSeed.Robust(), threadSafe)
+        {
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        public Xoshiro256StarStar(int seed)
+        {
+            Initialise(seed);
+        }
+
+        /// <summary>
+        /// Construct a new random number generator with random seed.
+        /// </summary>
+        /// <param name="seed">The seed value.</param>
+        /// <param name="threadSafe">if set to <c>true</c> , the class is thread safe.</param>
+        public Xoshiro256StarStar(int seed, bool threadSafe) : base(threadSafe)
+        {
+            Initialise(seed);
+        }
+
+        /// <summary>
+        /// Returns a random double-precision floating point number greater than or equal to 0.0, and less than 1.0.
+        /// </summary>
+        protected sealed override double DoSample()
+        {
+            // Note. Here we generate a random integer between 0 and 2^53-1 (i.e. 53 binary 1s) and multiply
+            // by the fractional unit value 1.0 / 2^53, thus the result has a max value of
+            // 1.0 - (1.0 / 2^53), or 0.99999999999999989 in decimal.
+            return (NextInnerULong() >> 11) * REAL_UNIT_UINT;
+        }
+
+        /// <summary>
+        /// Returns a random 32-bit signed integer greater than or equal to zero and less than <see cref="F:System.Int32.MaxValue"/>
+        /// </summary>
+        protected override int DoSampleInteger()
+        {
+            retry:
+            // Handle the special case where the value int.MaxValue is generated; this is outside 
+            // the range of permitted return values for this method. 
+            ulong rtn = NextInnerULong() & 0x7fff_ffffUL;
+            if (rtn == 0x7fff_ffffUL)
+            {
+                goto retry;
+            }
+            return (int)rtn;
+        }
+
+        /// <summary>
+        /// Fills the elements of a specified array of bytes with random numbers in full range, including zero and 255 (<see cref="F:System.Byte.MaxValue"/>).
+        /// </summary>
+        protected override void DoSampleBytes(byte[] buffer)
+        {
+            // For improved performance the below loop operates on these stack allocated copies of the heap variables.
+            // Notes. doing this means that these heavily used variables are located near to other local/stack variables,
+            // thus they will very likely be cached in the same CPU cache line.
+            ulong s0 = _s0;
+            ulong s1 = _s1;
+            ulong s2 = _s2;
+            ulong s3 = _s3;
+
+            int i = 0;
+
+            // Fill up the bulk of the buffer in chunks of 8 bytes at a time.
+            for (int bound = buffer.Length - 3; i < bound;)
+            {
+                // Generate 64 random bits.
+                ulong x = RotateLeft(s1 * 5, 7) * 9;
+
+                // Update PRNG state.
+                ulong t = s1 << 17;
+                s2 ^= s0;
+                s3 ^= s1;
+                s1 ^= s2;
+                s0 ^= s3;
+                s2 ^= t;
+                s3 = RotateLeft(s3, 45);
+
+                // Assign bits to a segment of 8 bytes.
+                buffer[i++] = (byte)x;
+                buffer[i++] = (byte)(x >> 8);
+                buffer[i++] = (byte)(x >> 16);
+                buffer[i++] = (byte)(x >> 24);
+                buffer[i++] = (byte)(x >> 32);
+                buffer[i++] = (byte)(x >> 40);
+                buffer[i++] = (byte)(x >> 48);
+                buffer[i++] = (byte)(x >> 56);
+            }
+
+            // Fill up any remaining bytes in the buffer.
+            if (i < buffer.Length)
+            {
+                // Generate 64 random bits.
+                ulong x = RotateLeft(s1 * 5, 7) * 9;
+
+                // Update PRNG state.
+                ulong t = s1 << 17;
+                s2 ^= s0;
+                s3 ^= s1;
+                s1 ^= s2;
+                s0 ^= s3;
+                s2 ^= t;
+                s3 = RotateLeft(s3, 45);
+
+                // Allocate one byte at a time until we reach the end of the buffer.
+                while (i < buffer.Length)
+                {
+                    buffer[i++] = (byte)x;
+                    x >>= 8;
+                }
+            }
+
+            // Update the state variables on the heap.
+            _s0 = s0;
+            _s1 = s1;
+            _s2 = s2;
+            _s3 = s3;
+        }
+
+        /// <summary>
+        /// Returns a random N-bit signed integer greater than or equal to zero and less than 2^N.
+        /// N (bit count) is expected to be greater than zero and less than 32 (not verified).
+        /// </summary>
+        protected override int DoSampleInt32WithNBits(int bitCount)
+        {
+            return (int)(NextInnerULong() >> (64 - bitCount));
+        }
+
+        /// <summary>
+        /// Returns a random N-bit signed long integer greater than or equal to zero and less than 2^N.
+        /// N (bit count) is expected to be greater than zero and less than 64 (not verified).
+        /// </summary>
+        protected override long DoSampleInt64WithNBits(int bitCount)
+        {
+            return (long)(NextInnerULong() >> (64 - bitCount));
+        }
+
+        private void Initialise(int seed)
+        {
+            // Notes.
+            // xoroshiro256** requires that at least one of the state variable be non-zero, use of splitmix64
+            // satisfies that requirement because its outputs are equidistributed, i.e. if a zero is output
+            // then the next zero will be after a further 2^64 outputs.
+            // Splitmix will also accept a zero input, thus all possible seeds can be accepted here and will
+            // all generate good initial state for xoshiro256**.
+            ulong longSeed = (ulong)seed;
+
+            // Use the splitmix64 RNG to hash the seed.
+            _s0 = Splitmix64(ref longSeed);
+            _s1 = Splitmix64(ref longSeed);
+            _s2 = Splitmix64(ref longSeed);
+            _s3 = Splitmix64(ref longSeed);
+        }
+        
+        private ulong NextInnerULong()
+        {
+            ulong s0 = _s0;
+            ulong s1 = _s1;
+            ulong s2 = _s2;
+            ulong s3 = _s3;
+
+            ulong result = RotateLeft(s1 * 5, 7) * 9;
+
+            ulong t = s1 << 17;
+            s2 ^= s0;
+            s3 ^= s1;
+            s1 ^= s2;
+            s0 ^= s3;
+            s2 ^= t;
+            s3 = RotateLeft(s3, 45);
+
+            _s0 = s0;
+            _s1 = s1;
+            _s2 = s2;
+            _s3 = s3;
+
+            return result;
+        }
+
+        /// <summary>
+        /// Fills an array with random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        /// <remarks>Supports being called in parallel from multiple threads.</remarks>
+        public static void Doubles(double[] values, int seed)
+        {
+            // Init state.
+            ulong longSeed = (ulong)seed;
+            ulong s0 = Splitmix64(ref longSeed);
+            ulong s1 = Splitmix64(ref longSeed);
+            ulong s2 = Splitmix64(ref longSeed);
+            ulong s3 = Splitmix64(ref longSeed);
+
+            for (int i = 0; i < values.Length; i++)
+            {
+                // Generate sample.
+                values[i] = ((RotateLeft(s1 * 5, 7) * 9) >> 11) * REAL_UNIT_UINT;
+
+                // Update PRNG state.
+                ulong t = s1 << 17;
+                s2 ^= s0;
+                s3 ^= s1;
+                s1 ^= s2;
+                s0 ^= s3;
+                s2 ^= t;
+                s3 = RotateLeft(s3, 45);
+            }
+        }
+
+        /// <summary>
+        /// Returns an array of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        /// <remarks>Supports being called in parallel from multiple threads.</remarks>
+        [TargetedPatchingOptOut("Performance critical to inline this type of method across NGen image boundaries")]
+        public static double[] Doubles(int length, int seed)
+        {
+            var data = new double[length];
+            Doubles(data, seed);
+            return data;
+        }
+
+        /// <summary>
+        /// Returns an infinite sequence of random numbers greater than or equal to 0.0 and less than 1.0.
+        /// </summary>
+        /// <remarks>Supports being called in parallel from multiple threads, but the result must be enumerated from a single thread each.</remarks>
+        public static IEnumerable<double> DoubleSequence(int seed)
+        {
+            // Init state.
+            ulong longSeed = (ulong)seed;
+            ulong s0 = Splitmix64(ref longSeed);
+            ulong s1 = Splitmix64(ref longSeed);
+            ulong s2 = Splitmix64(ref longSeed);
+            ulong s3 = Splitmix64(ref longSeed);
+
+            while (true)
+            {
+                // Generate sample.
+                double x = ((RotateLeft(s1 * 5, 7) * 9) >> 11) * REAL_UNIT_UINT;
+
+                // Update PRNG state.
+                ulong t = s1 << 17;
+                s2 ^= s0;
+                s3 ^= s1;
+                s1 ^= s2;
+                s0 ^= s3;
+                s2 ^= t;
+                s3 = RotateLeft(s3, 45);
+
+                // Yield sample.
+                yield return x;
+            }
+        }
+
+        /// <summary>
+        /// Splitmix64 RNG.
+        /// </summary>
+        /// <param name="x">RNG state. This can take any value, including zero.</param>
+        /// <returns>A new random UInt64.</returns>
+        /// <remarks>
+        /// Splitmix64 produces equidistributed outputs, thus if a zero is generated then the
+        /// next zero will be after a further 2^64 outputs.
+        /// </remarks>
+        private static ulong Splitmix64(ref ulong x)
+        {
+            ulong z = (x += 0x9E3779B97F4A7C15UL);
+            z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9UL;
+            z = (z ^ (z >> 27)) * 0x94D049BB133111EBUL;
+            return z ^ (z >> 31);
+        }
+
+        private static ulong RotateLeft(ulong x, int k)
+        {
+            // Note. RyuJIT will compile this to a single rotate CPU instruction (as of about .NET 4.6.1 and dotnet core 2.0).
+            return (x << k) | (x >> (64 - k));
+        }
+    }
+}