Add pseudorandom number generator WyRand

Author: wigging

Package.swift

@@ -12,10 +12,13 @@ let package = Package(
     targets: [
         .target(
             name: "Numerix",
-            dependencies: ["VectorModule", "MatrixModule", "ShapedArrayModule"],
+            dependencies: ["RandomModule", "VectorModule", "MatrixModule", "ShapedArrayModule"],
             cxxSettings: [.define("ACCELERATE_NEW_LAPACK", to: "1"), .define("ACCELERATE_LAPACK_ILP64", to: "1")],
             linkerSettings: [.linkedFramework("Accelerate")]
         ),
+        .target(
+            name: "RandomModule"
+        ),
         .target(
             name: "VectorModule",
             cxxSettings: [.define("ACCELERATE_NEW_LAPACK", to: "1"), .define("ACCELERATE_LAPACK_ILP64", to: "1")],

Sources/Numerix/Numerix.swift

@@ -2,6 +2,7 @@
 Exported modules for the Numerix package.
 */
 
+@_exported import RandomModule
 @_exported import VectorModule
 @_exported import MatrixModule
 @_exported import ShapedArrayModule

Sources/RandomModule/Random.swift

@@ -0,0 +1,41 @@
+/*
+Protocol for random vector.
+Random vector protocol and init.
+*/
+
+import VectorModule
+
+@_documentation(visibility: private)
+public protocol Random {
+    static func random(size: Int, seed: UInt64?) -> Vector<Self>
+}
+
+@_documentation(visibility: private)
+extension Float: Random {
+    public static func random(size: Int, seed: UInt64?) -> Vector<Float> {
+        var vec = Vector<Float>(size: size)
+        var rng = WyRand(seed: seed)
+        for i in 0..<size {
+            vec[i] = rng.nextUniform()
+        }
+        return vec
+    }
+}
+
+@_documentation(visibility: private)
+extension Double: Random {
+    public static func random(size: Int, seed: UInt64?) -> Vector<Double> {
+        var vec = Vector<Double>(size: size)
+        var rng = WyRand(seed: seed)
+        for i in 0..<size {
+            vec[i] = rng.nextUniform()
+        }
+        return vec
+    }
+}
+
+extension Vector {
+    public static func random(size: Int, seed: UInt64? = nil) -> Vector where Scalar: Random {
+        Scalar.random(size: size, seed: seed)
+    }
+}

Sources/RandomModule/WyRand.swift

@@ -0,0 +1,35 @@
+/*
+WyRand pseudorandom number generator (PRNG) based on wyrand and wyhash functions by Wang Yi.
+
+Original code for wyrand and wyhash functions
+https://github.com/wangyi-fudan/wyhash
+
+Also see Daniel Lemire's article about wyhash
+https://lemire.me/blog/2019/03/19/the-fastest-conventional-random-number-generator-that-can-pass-big-crush/
+*/
+
+import Foundation
+
+public struct WyRand: RandomNumberGenerator {
+    private var state : UInt64
+
+    public init(seed: UInt64? = nil) {
+        state = seed ?? UInt64(abs(UUID().hashValue))
+    }
+
+    public mutating func next() -> UInt64 {
+        state &+= 0x2d358dccaa6c78a5
+        let mul = state.multipliedFullWidth(by: state ^ 0x8bb84b93962eacc9)
+        return mul.low ^ mul.high
+    }
+
+    public mutating func nextUniform() -> Float {
+        let x = self.next()
+        return Float(x >> 40) * 0x1.0p-24
+    }
+
+    public mutating func nextUniform() -> Double {
+        let x = self.next()
+        return Double(x >> 11) * 0x1.0p-53
+    }
+}

Tests/RandomTests.swift

@@ -0,0 +1,36 @@
+/*
+Tests for random numbers and random vectors.
+*/
+
+import Testing
+@testable import Numerix
+
+struct RandomTests {
+
+    @Test func randomWyRand() {
+        var rand = WyRand(seed: 12345)
+        #expect(rand.next() == 13157676964440363053)
+
+        let a: Float = rand.nextUniform()
+        #expect(a == 0.9815675)
+
+        let b: Double = rand.nextUniform()
+        #expect(b ==  0.7952823641045631)
+    }
+
+    @Test func randomVector() {
+        // No seed
+        let a = Vector<Float>.random(size: 3)
+        #expect(a[0] < 1.0)
+
+        let b = Vector<Double>.random(size: 3)
+        #expect(b[0] < 1.0)
+
+        // With seed
+        let c = Vector<Float>.random(size: 4, seed: 12345)
+        #expect(c == Vector<Float>([0.71327907, 0.9815675, 0.79528236, 0.25364798]))
+
+        let d = Vector<Double>.random(size: 4, seed: 12345)
+        #expect(d == Vector<Double>([0.7132790974854359, 0.9815675033121781, 0.7952823641045631, 0.25364800714420743]))
+    }
+}