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.