Wrapping up my series on the, xoshiro/xoroshiro algorithms, in this article I present a package to return 32-bit values using the following variants.
- xoroshiro64*
- xoroshiro64**
- xoshiro128+
- xoshiro128++
- xoshiro128**
The numbers in the names refer to the size of the state array that each method operates on. The 64s operate on two 32-bit values, while the 128s operate on four 32-bit values.
Using each of them follows the same pattern:
- Choose the algorithm mode (the package defaults to xoshiro64*)
- Seed the state array
- Retrieve the next number from the chosen algorithm as many times as needed
For example…
SQL> BEGIN
2 xoshiro32.set_mode('xoroshiro64**');
3 xoshiro32.set_state(xoshiro32.seed_tab(12345, 67890));
4
5 FOR i IN 0 .. 9
6 LOOP
7 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.get_next, 10));
8 END LOOP;
9 END;
10 /
0 3157960260
1 4142509522
2 1831851427
3 506054173
4 2910589752
5 1819521659
6 3282141937
7 2257682835
8 2133372007
9 3757018772
While all of the 64-bit variants supported jump and long_jump functions, only those with 128-bit state arrays support them in the 32-bit variants. in the state by simulating many (from 2^64 to 2^768 depending on the mode and jump size) calls of the get_next function. As noted by Blackman and Vigna, these are typically used to generate non-overlapping sequences for parallel processing. In single threaded use there would be no need for them.
The 64-bit variants included splitmix64, which I ported from Vigna’s original c implementation. No corresponding method was provided for the 32-bit variants, so I went back to the original SplitMix paper by Steele, Lea, and Flood. Mirroring the logic of splitmix64, I also used a 32-bit version of the MurmurHash3 by Appleby to calculate the Golden Gamma.
SQL> BEGIN 2 xoshiro32.splitmix32_set_state(12345); 3 4 FOR i IN 0 .. 9 5 LOOP 6 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.splitmix32_next, 10)); 7 END LOOP; 8 END; 9 / 0 1200724404 1 818072533 2 996137225 3 2397394836 4 4079075752 5 2274189806 6 2795887828 7 4161515127 8 3291005408 9 722528451
Just as with splitmix64, using splitmix32 directly to generate pseudorandom numbers is not the typical usage . Rather, splitmix32 is best used as a means of seeding the state arrays of the other algorithms. Below, I set the splitmix32 state, and then create a seed of 4 values for the xoshiro128** generator.
SQL> DECLARE
2 v_seed xoshiro32.seed_tab := xoshiro32.seed_tab();
3 BEGIN
4 xoshiro32.splitmix32_set_state(12345);
5 v_seed.EXTEND(4);
6
7 FOR i IN 1 .. 4
8 LOOP
9 v_seed(i) := xoshiro32.splitmix32_next;
10 END LOOP;
11
12 xoshiro32.set_mode('xoshiro128**');
13 xoshiro32.set_state(v_seed);
14
15 FOR i IN 0 .. 9
16 LOOP
17 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.get_next, 10));
18 END LOOP;
19 END;
20 /
0 518667457
1 440444462
2 4232892992
3 3757857622
4 3939018813
5 1334683535
6 3795058715
7 2092637810
8 2829112157
9 779180383
Following is the complete package for all of the 32-bit algorithms.
CREATE OR REPLACE PACKAGE xoshiro32
IS
-- Based on xoshiro/xoroshiro family of pseudorandom number generators
-- by by David Blackman and Sebastiano Vigna
-- original https://prng.di.unimi.it
-- .///.
-- (0 o)
---------------0000--(_)--0000---------------
--
-- Sean D. Stuber
-- sean.stuber@gmail.com
--
-- oooO Oooo
--------------( )-----( )---------------
-- \ ( ) /
-- \_) (_/
-- It would be nice to define an explicit RANGE 0..4294967295 for the subtype
-- but as of 23ai, only subtypes of PLS_INTEGER may have a range.
SUBTYPE uint32_t IS NUMBER(10, 0);
TYPE seed_tab IS TABLE OF uint32_t;
SUBTYPE bits_t IS SIMPLE_INTEGER RANGE 0 .. 31;
SUBTYPE mode_t IS SIMPLE_INTEGER RANGE 13 .. 17;
-- These do not support jump and long_jump operations
c_64star CONSTANT mode_t := 13; -- xoroshiro64*
c_64starstar CONSTANT mode_t := 14; -- xoroshiro64**
-- These do support jump and long_jump operations
c_128plus CONSTANT mode_t := 15; -- xoshiro128+
c_128plusplus CONSTANT mode_t := 16; -- xoshiro128++
c_128starstar CONSTANT mode_t := 17; -- xoshiro128**
--Exceptions
invalid_mode EXCEPTION;
-- Select a mode by numeric value or by name
PROCEDURE set_mode(p_mode IN mode_t);
PROCEDURE set_mode(p_mode IN VARCHAR2);
-- Return the current mode number or name
FUNCTION get_mode
RETURN mode_t;
FUNCTION get_mode_name
RETURN VARCHAR2;
-- Given a set of seed values, initialize the state array for the current mode size
-- Different modes for a particular size can share the same state array.
PROCEDURE set_state(p_seed seed_tab);
-- Return the next pseudorandom generator value for the current mode
FUNCTION get_next
RETURN uint32_t;
-- Simulate many calls to "get_next".
-- Typically used to create non-overlapping sequences from the same initial state.
-- xoroshiro64* and xoroshiro64** do not support jump or long_jump operations
PROCEDURE jump;
-- Simulate even more calls to "get_next".
--Also used to create non-overlapping sequences from the same initial state.
-- xoroshiro64* and xoroshiro64** do not support jump or long_jump operations
PROCEDURE long_jump;
-- SplitMix32 isn't part of the xoshiro/xoroshiro family
-- However using a function like it to fill the state arrays instead of a simple
-- integer is recommended by the original authors.
PROCEDURE splitmix32_set_state(p_seed uint32_t);
FUNCTION splitmix32_next
RETURN uint32_t;
END;
/
As a comparison reference for myself and readers, I included snippets of the original c code within the comments of the package body for the various procedures and functions of each mode.
CREATE OR REPLACE PACKAGE BODY xoshiro32
IS
c_32bitmask CONSTANT INTEGER := POWER(2, 32) - 1;
SUBTYPE array_64index_t IS PLS_INTEGER RANGE 0 .. 1;
SUBTYPE array_128index_t IS PLS_INTEGER RANGE 0 .. 3;
TYPE state64_array_t IS TABLE OF uint32_t
INDEX BY array_64index_t;
TYPE state128_array_t IS TABLE OF uint32_t
INDEX BY array_128index_t;
g_mode mode_t := c_64star;
g_splitmix32_state uint32_t;
g_64_state state64_array_t;
g_128_state state128_array_t;
FUNCTION bitor32(p_a IN uint32_t, p_b IN uint32_t)
RETURN uint32_t
IS
BEGIN
RETURN BITAND(p_a + p_b - BITAND(p_a, p_b), c_32bitmask);
END bitor32;
FUNCTION bitxor32(p_a IN uint32_t, p_b IN uint32_t)
RETURN uint32_t
IS
BEGIN
RETURN BITAND((p_a + p_b) - (BITAND(p_a, p_b) * 2), c_32bitmask);
END bitxor32;
FUNCTION shl32(p_n IN uint32_t, p_bits IN bits_t)
RETURN uint32_t
IS
BEGIN
RETURN CASE p_bits WHEN 0 THEN p_n ELSE BITAND(p_n * (2 ** p_bits), c_32bitmask) END;
END shl32;
FUNCTION shr32(p_n IN uint32_t, p_bits IN bits_t)
RETURN uint32_t
IS
BEGIN
RETURN CASE p_bits WHEN 0 THEN p_n ELSE BITAND(FLOOR(p_n / (2 ** p_bits)), c_32bitmask) END;
END shr32;
FUNCTION rotl32(p_x uint32_t, p_k IN bits_t)
RETURN uint32_t
IS
-- static inline uint32_t rotl(const uint32_t x, int k) {
-- return (x << k) | (x >> (32 - k));
-- }
BEGIN
RETURN CASE WHEN p_k = 0 THEN p_x ELSE bitor32(shl32(p_x, p_k), shr32(p_x, 32 - p_k)) END;
END rotl32;
-----------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------
FUNCTION next64
RETURN uint32_t
IS
-- uint32_t next(void) {
-- const uint32_t s0 = s[0];
-- uint32_t s1 = s[1];
-- const uint32_t result = s0 * 0x9E3779BB; -- xoroshiro64*
-- const uint32_t result = rotl(s0 * 0x9E3779BB, 5) * 5; -- xoroshiro64**
--
-- s1 ^= s0;
-- s[0] = rotl(s0, 26) ^ s1 ^ (s1 << 9); // a, b
-- s[1] = rotl(s1, 13); // c
--
-- return result;
-- }
v_s0 uint32_t := g_64_state(0);
v_s1 uint32_t := g_64_state(1);
v_result uint32_t
:= BITAND(
CASE g_mode
WHEN c_64star
THEN
v_s0 * TO_NUMBER('9E3779BB', 'xxxxxxxx')
WHEN c_64starstar
THEN
rotl32(BITAND(v_s0 * TO_NUMBER('9E3779BB', 'xxxxxxxx'), c_32bitmask), 5) * 5
END,
c_32bitmask);
BEGIN
v_s1 := bitxor32(v_s1, v_s0);
g_64_state(0) := bitxor32(bitxor32(rotl32(v_s0, 26), v_s1), shl32(v_s1, 9));
g_64_state(1) := rotl32(v_s1, 13);
RETURN v_result;
END next64;
FUNCTION next128
RETURN uint32_t
IS
-- uint32_t next(void) {
-- const uint32_t result = s[0] + s[3]; -- xoshiro128+
-- const uint32_t result = s[0] + s[3]; -- xoshiro128++
-- const uint32_t result = rotl(s[1] * 5, 7) * 9; -- xoshiro128**
--
-- const uint32_t t = s[1] << 9;
--
-- s[2] ^= s[0];
-- s[3] ^= s[1];
-- s[1] ^= s[2];
-- s[0] ^= s[3];
--
-- s[2] ^= t;
--
-- s[3] = rotl(s[3], 11);
--
-- return result;
-- }
v_result uint32_t
:= BITAND(
CASE g_mode
WHEN c_128plus
THEN
g_128_state(0) + g_128_state(3)
WHEN c_128plusplus
THEN
rotl32(BITAND(g_128_state(0) + g_128_state(3), c_32bitmask), 7) + g_128_state(0)
WHEN c_128starstar
THEN
rotl32(BITAND(g_128_state(1) * 5, c_32bitmask), 7) * 9
END,
c_32bitmask);
v_t uint32_t := shl32(g_128_state(1), 9);
BEGIN
g_128_state(2) := bitxor32(g_128_state(2), g_128_state(0));
g_128_state(3) := bitxor32(g_128_state(3), g_128_state(1));
g_128_state(1) := bitxor32(g_128_state(1), g_128_state(2));
g_128_state(0) := bitxor32(g_128_state(0), g_128_state(3));
g_128_state(2) := bitxor32(g_128_state(2), v_t);
g_128_state(3) := rotl32(g_128_state(3), 11);
RETURN v_result;
END next128;
PROCEDURE jump
IS
-- /* This is the jump function for the generator. It is equivalent
-- to 2^64 calls to next(); it can be used to generate 2^64
-- non-overlapping subsequences for parallel computations. */
--
-- void jump(void) {
-- static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
--
-- uint32_t s0 = 0;
-- uint32_t s1 = 0;
-- uint32_t s2 = 0;
-- uint32_t s3 = 0;
-- for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
-- for(int b = 0; b < 32; b++) {
-- if (JUMP[i] & UINT32_C(1) << b) {
-- s0 ^= s[0];
-- s1 ^= s[1];
-- s2 ^= s[2];
-- s3 ^= s[3];
-- }
-- next();
-- }
--
-- s[0] = s0;
-- s[1] = s1;
-- s[2] = s2;
-- s[3] = s3;
-- }
v_jump state128_array_t
:= state128_array_t(0 => TO_NUMBER('8764000b', 'xxxxxxxx'),
1 => TO_NUMBER('f542d2d3', 'xxxxxxxx'),
2 => TO_NUMBER('6fa035c3', 'xxxxxxxx'),
3 => TO_NUMBER('77f2db5b', 'xxxxxxxx'));
v_s0 uint32_t := 0;
v_s1 uint32_t := 0;
v_s2 uint32_t := 0;
v_s3 uint32_t := 0;
v_dummy uint32_t;
BEGIN
CASE
WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
THEN
FOR i IN 0 .. 3
LOOP
FOR b IN 0 .. 31
LOOP
IF BITAND(v_jump(i), shl32(1, b)) != 0
THEN
v_s0 := bitxor32(v_s0, g_128_state(0));
v_s1 := bitxor32(v_s1, g_128_state(1));
v_s2 := bitxor32(v_s2, g_128_state(2));
v_s3 := bitxor32(v_s3, g_128_state(3));
END IF;
v_dummy := next128;
END LOOP;
END LOOP;
g_128_state(0) := v_s0;
g_128_state(1) := v_s1;
g_128_state(2) := v_s2;
g_128_state(3) := v_s3;
ELSE
-- 64-bit state modes don't support jump and long_jump operations
RAISE invalid_mode;
END CASE;
END jump;
PROCEDURE long_jump
IS
-- /* This is the long-jump function for the generator. It is equivalent to
-- 2^96 calls to next(); it can be used to generate 2^32 starting points,
-- from each of which jump() will generate 2^32 non-overlapping
-- subsequences for parallel distributed computations. */
--
-- void long_jump(void) {
-- static const uint32_t LONG_JUMP[] = { 0xb523952e, 0x0b6f099f, 0xccf5a0ef, 0x1c580662 };
--
-- uint32_t s0 = 0;
-- uint32_t s1 = 0;
-- uint32_t s2 = 0;
-- uint32_t s3 = 0;
-- for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++)
-- for(int b = 0; b < 32; b++) {
-- if (LONG_JUMP[i] & UINT32_C(1) << b) {
-- s0 ^= s[0];
-- s1 ^= s[1];
-- s2 ^= s[2];
-- s3 ^= s[3];
-- }
-- next();
-- }
--
-- s[0] = s0;
-- s[1] = s1;
-- s[2] = s2;
-- s[3] = s3;
-- }
v_jump state128_array_t
:= state128_array_t(0 => TO_NUMBER('b523952e', 'xxxxxxxx'),
1 => TO_NUMBER('0b6f099f', 'xxxxxxxx'),
2 => TO_NUMBER('ccf5a0ef', 'xxxxxxxx'),
3 => TO_NUMBER('1c580662', 'xxxxxxxx'));
v_s0 uint32_t := 0;
v_s1 uint32_t := 0;
v_s2 uint32_t := 0;
v_s3 uint32_t := 0;
v_dummy uint32_t;
BEGIN
CASE
WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
THEN
FOR i IN 0 .. 3
LOOP
FOR b IN 0 .. 31
LOOP
IF BITAND(v_jump(i), shl32(1, b)) != 0
THEN
v_s0 := bitxor32(v_s0, g_128_state(0));
v_s1 := bitxor32(v_s1, g_128_state(1));
v_s2 := bitxor32(v_s2, g_128_state(2));
v_s3 := bitxor32(v_s3, g_128_state(3));
END IF;
v_dummy := next128;
END LOOP;
END LOOP;
g_128_state(0) := v_s0;
g_128_state(1) := v_s1;
g_128_state(2) := v_s2;
g_128_state(3) := v_s3;
ELSE
-- 64-bit state modes don't support jump and long_jump operations
RAISE invalid_mode;
END CASE;
END long_jump;
-----------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------
PROCEDURE set_mode(p_mode IN mode_t)
IS
BEGIN
g_mode := p_mode;
END set_mode;
PROCEDURE set_mode(p_mode IN VARCHAR2)
IS
BEGIN
set_mode(
CASE LOWER(p_mode)
WHEN 'xoroshiro64*' THEN c_64star
WHEN 'xoroshiro64**' THEN c_64starstar
WHEN 'xoshiro128+' THEN c_128plus
WHEN 'xoshiro128++' THEN c_128plusplus
WHEN 'xoshiro128**' THEN c_128starstar
END);
END set_mode;
FUNCTION get_mode
RETURN mode_t
IS
BEGIN
RETURN g_mode;
END get_mode;
FUNCTION get_mode_name
RETURN VARCHAR2
IS
BEGIN
RETURN CASE g_mode
WHEN c_64star THEN 'xoroshiro64*'
WHEN c_64starstar THEN 'xoroshiro64**'
WHEN c_128plus THEN 'xoshiro128+'
WHEN c_128plusplus THEN 'xoshiro128++'
WHEN c_128starstar THEN 'xoshiro128**'
END;
END get_mode_name;
FUNCTION get_next
RETURN uint32_t
IS
BEGIN
RETURN CASE
WHEN g_mode IN (c_64star, c_64starstar) THEN next64
WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN next128
END;
END get_next;
PROCEDURE set_state(p_seed IN seed_tab)
IS
BEGIN
CASE
WHEN g_mode IN (c_64star, c_64starstar)
THEN
g_64_state := state64_array_t();
FOR i IN 0 .. 1
LOOP
g_64_state(i) := p_seed(i + 1);
END LOOP;
WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
THEN
g_128_state := state128_array_t();
FOR i IN 0 .. 3
LOOP
g_128_state(i) := p_seed(i + 1);
END LOOP;
ELSE
RAISE invalid_mode;
END CASE;
END set_state;
-----------------------------------------------------------------------------------------------------------------
-----------------------------------------------------------------------------------------------------------------
PROCEDURE splitmix32_set_state(p_seed uint32_t)
IS
BEGIN
g_splitmix32_state := p_seed;
END splitmix32_set_state;
FUNCTION splitmix32_next
RETURN uint32_t
IS
-- splitmix32 is implemented using MurmurHash3
-- with the extra step of adding the "Golden Gamma"
-- This modification is base on the SplitMix algorithm described by
-- Guy L. Steele Jr., Doug Lea, and Christine H. Flood
-- https://gee.cs.oswego.edu/dl/papers/oopsla14.pdf
--
-- The MurmurHash3 algorithm by Austin Appleby can be found here.
-- https://github.com/aappleby/smhasher/blob/master/src/MurmurHash3.cpp
--
-- The Golden Gamma is the closest odd number to 2^32/golden_ratio
-- Where the golden ratio is (1+sqrt(5))/2
c_golden_gamma CONSTANT INTEGER := 2654435769; -- 0x9e3779b9
v_z NUMBER;
BEGIN
-- In Appleby's code, the state is passed in as a parameter.
-- In this package, the state is maintained as a package variable.
g_splitmix32_state := BITAND(g_splitmix32_state + c_golden_gamma, c_32bitmask);
v_z := g_splitmix32_state;
v_z := bitxor32(v_z, shr32(v_z, 16));
v_z := v_z * TO_NUMBER('85ebca6b', 'xxxxxxxx');
v_z := BITAND(v_z, c_32bitmask);
v_z := bitxor32(v_z, shr32(v_z, 13));
v_z := v_z * TO_NUMBER('c2b2ae35', 'xxxxxxxx');
v_z := BITAND(v_z, c_32bitmask);
RETURN BITAND(bitxor32(v_z, shr32(v_z, 16)), c_32bitmask);
END splitmix32_next;
END;
/
Because the seed numbers and intermediate results are smaller, these algorithm variants can use the native BITAND function and don’t requite the BIG_BITAND function for extremely large values.
This concludes the series and I hope you find the code interesting and useful. As mentioned in my 64-bit article, I intentionally split the 32-bit and 64-bit algorithms into distinct packages on the assumption they would be unlikely to be used together. But, if you find you do want them to be combined, I defined the mode constants such that they are mutually exclusive thus could be combined without fear of ambiguity from overlapping values.
Questions and comments, as always, are welcome.