Continuing in the path of previous two posts, the package below supports all of the xoshiro/xoroshiro algorithms returning 64-bit values.
- xoroshiro128+
- xoroshiro128++
- xoroshiro128**
- xoshiro256+
- xoshiro256++
- xoshiro256**
- xoshiro512+
- xoshiro512++
- xoshiro512**
- xoroshiro1024*
- xoroshiro1024++
- xoroshiro1024**
The numbers in the names refer to the size of the state array that each method operates on.
Using each of them follows the same pattern:
- Choose the algorithm mode (the package defaults to xoshiro256++)
- Seed the state array
- Retrieve the next number from the chosen algorithm as many times as needed
For example…
SQL> BEGIN 2 xoshiro.set_mode('xoroshiro128++'); 3 xoshiro.set_state(xoshiro.seed_tab(12345, 67890)); 4 5 FOR i IN 0 .. 9 6 LOOP 7 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro.get_next, 10)); 8 END LOOP; 9 END; 10 / 0 1051657426 1 9791286188 2 9989489894 3 6590593871 4 6705930516 5 1618779258 6 8800234758 7 2334047907 8 8330851765 9 3477959413
Optionally, you can call the jump or long_jump functions to jump ahead 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.
In addition to these I also included the splitmix64 function which is, itself, a pseudo-random number generator and can be used in a similar manner. Set the seed, and then call the splitmix64_next to retrieve the next number.
SQL> BEGIN 2 xoshiro.splitmix64_set_state(12345); 3 4 FOR i IN 0 .. 9 5 LOOP 6 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro.splitmix64_next, 10)); 7 END LOOP; 8 END; 9 / 0 2454886589 1 3778200017 2 2205171434 3 3248800117 4 9350289611 5 6217189988 6 2262534019 7 7959005890 8 8850488307 9 1600295491
This is not the typical usage though. Rather, splitmix64 is suggested as a means of seeding the state arrays of the other algorithms. Below, I set the splitmix64 state, and then create a seed of 4 values for the xoshiro256** generator.
SQL> DECLARE 2 v_seed xoshiro.seed_tab := xoshiro.seed_tab(); 3 BEGIN 4 xoshiro.splitmix64_set_state(12345); 5 v_seed.EXTEND(4); 6 7 FOR i IN 1 .. 4 8 LOOP 9 v_seed(i) := xoshiro.splitmix64_next; 10 END LOOP; 11 12 xoshiro.set_mode('xoshiro256**'); 13 xoshiro.set_state(v_seed); 14 15 FOR i IN 0 .. 9 16 LOOP 17 DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro.get_next, 10)); 18 END LOOP; 19 END; 20 / 0 1372083882 1 2398916695 2 1777038484 3 8917177268 4 1024131604 5 1969754298 6 2947371003 7 5456629693 8 7119811276 9 1679784635
Following is the complete package for all of the 64-bit algorithms.
CREATE OR REPLACE PACKAGE xoshiro 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..18446744073709551615 for the subtype -- but as of 23ai, only subtypes of PLS_INTEGER may have a range. SUBTYPE uint64_t IS NUMBER(20, 0); TYPE seed_tab IS TABLE OF uint64_t; SUBTYPE bits_t IS SIMPLE_INTEGER RANGE 0 .. 63; SUBTYPE mode_t IS SIMPLE_INTEGER RANGE 1 .. 12; c_128plus CONSTANT mode_t := 1; -- xoroshiro128+ c_128plusplus CONSTANT mode_t := 2; -- xoroshiro128++ c_128starstar CONSTANT mode_t := 3; -- xoroshiro128** c_256plus CONSTANT mode_t := 4; -- xoshiro256+ c_256plusplus CONSTANT mode_t := 5; -- xoshiro256++ c_256starstar CONSTANT mode_t := 6; -- xoshiro256** c_512plus CONSTANT mode_t := 7; -- xoshiro512+ c_512plusplus CONSTANT mode_t := 8; -- xoshiro512++ c_512starstar CONSTANT mode_t := 9; -- xoshiro512** c_1024star CONSTANT mode_t := 10; -- xoroshiro1024* c_1024plusplus CONSTANT mode_t := 11; -- xoroshiro1024++ c_1024starstar CONSTANT mode_t := 12; -- xoroshiro1024** --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 IN seed_tab); -- Return the next pseudorandom generator value for the current mode FUNCTION get_next RETURN uint64_t; -- Simulate many calls to "get_next". Typically used to create non-overlapping sequences from the same initial state. PROCEDURE jump; -- Simulate even more calls to "get_next". Also used to create non-overlapping sequences from the same initial state. PROCEDURE long_jump; -- SplitMix64 isn't part of the xoshiro/xoroshiro family -- However using it to fill the state arrays instead of a simple integer -- is recommended by the original authors. PROCEDURE splitmix64_set_state(p_seed uint64_t); FUNCTION splitmix64_next RETURN uint64_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 xoshiro IS c_64bitmask CONSTANT uint64_t := POWER(2, 64) - 1; SUBTYPE array_128index_t IS PLS_INTEGER RANGE 0 .. 1; SUBTYPE array_256index_t IS PLS_INTEGER RANGE 0 .. 3; SUBTYPE array_512index_t IS PLS_INTEGER RANGE 0 .. 7; SUBTYPE array_1024index_t IS PLS_INTEGER RANGE 0 .. 15; TYPE state128_array_t IS TABLE OF uint64_t INDEX BY array_128index_t; TYPE state256_array_t IS TABLE OF uint64_t INDEX BY array_256index_t; TYPE state512_array_t IS TABLE OF uint64_t INDEX BY array_512index_t; TYPE state1024_array_t IS TABLE OF uint64_t INDEX BY array_1024index_t; g_mode mode_t := c_256plusplus; g_splitmix64_state uint64_t; g_128_state state128_array_t; -- xoroshiro128+, xoroshiro128++, xoroshiro128** g_256_state state256_array_t; -- xoshiro256+, xoshiro256++, xoshiro256** g_512_state state512_array_t; -- xoshiro512+, xoshiro512++, xoshiro512** g_1024_state state1024_array_t; -- xoroshiro1024*, xoroshiro1024++, xoroshiro1024** -- Used to index into the state array for the 1024bit modes g_1024_p array_1024index_t := 0; FUNCTION bitor(p_a IN uint64_t, p_b IN uint64_t) RETURN uint64_t IS BEGIN RETURN BITAND(p_a + p_b - BITAND(p_a, p_b), c_64bitmask); END bitor; FUNCTION bitxor(p_a IN uint64_t, p_b IN uint64_t) RETURN uint64_t IS BEGIN RETURN BITAND((p_a + p_b) - (BITAND(p_a, p_b) * 2), c_64bitmask); END bitxor; FUNCTION big_bitand(p_a IN NUMBER, p_b IN NUMBER) RETURN NUMBER IS v_result NUMBER; v_temp NUMBER; v_a_shift NUMBER; v_b_shift NUMBER; v_a_first_bit NUMBER; v_b_first_bit NUMBER; BEGIN -- Oracle's native BITAND function has an upper limit of power(2,127)-1 CASE WHEN p_a < POWER(2, 127) AND p_b < POWER(2, 127) THEN v_result := BITAND(p_a, p_b); ELSE -- Divide by 2 (shift right by 1 bit) to put the values in range -- BITAND the reduced values, and then shift back left by multiplying by 2. -- Then add the 1-bit that was dropped off by shifting right. -- Since we can't use bitand, we have to use modulo division to find the 1-bit value. v_a_first_bit := MOD(p_a, 2); v_b_first_bit := MOD(p_b, 2); -- for large numbers, MOD can fail and return -1 instead of 1, so check for that. -- when MOD fails that also means division fails CASE WHEN v_a_first_bit = -1 THEN v_a_first_bit := 1; v_a_shift := FLOOR(p_a / 2) - 1; ELSE v_a_shift := FLOOR(p_a / 2); END CASE; CASE WHEN v_b_first_bit = -1 THEN v_b_first_bit := 1; v_b_shift := FLOOR(p_b / 2) - 1; ELSE v_b_shift := FLOOR(p_b / 2); END CASE; v_temp := BITAND(v_a_shift, v_b_shift) * 2; v_result := v_temp + CASE WHEN v_a_first_bit = 1 AND v_b_first_bit = 1 THEN 1 ELSE 0 END; END CASE; RETURN v_result; END big_bitand; FUNCTION shl64(p_n IN uint64_t, p_bits IN bits_t) RETURN uint64_t IS BEGIN -- Apply 64-bit mask to ensure results are within proper range RETURN CASE p_bits WHEN 0 THEN p_n ELSE big_bitand(p_n * (2 ** p_bits), c_64bitmask) END; END shl64; FUNCTION shr64(p_n IN uint64_t, p_bits IN bits_t) RETURN uint64_t IS BEGIN -- Apply 64-bit mask to ensure results are within proper range RETURN CASE p_bits WHEN 0 THEN p_n ELSE BITAND(FLOOR(p_n / (2 ** p_bits)), c_64bitmask) END; END shr64; FUNCTION rotl(p_x uint64_t, p_k IN bits_t) RETURN uint64_t IS -- static inline uint64_t rotl(const uint64_t x, int k) { -- return (x << k) | (x >> (64 - k)); -- } BEGIN RETURN CASE WHEN p_k = 0 THEN p_x ELSE bitor(shl64(p_x, p_k), shr64(p_x, 64 - p_k)) END; END rotl; ----------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------- FUNCTION next128plus RETURN uint64_t IS -- uint64_t next(void) { -- const uint64_t s0 = s[0]; -- uint64_t s1 = s[1]; -- const uint64_t result = s0 + s1; -- -- s1 ^= s0; -- s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b -- s[1] = rotl(s1, 37); // c -- -- return result; -- } v_s0 uint64_t := g_128_state(0); v_s1 uint64_t := g_128_state(1); v_result uint64_t := BITAND(v_s0 + v_s1, c_64bitmask); BEGIN v_s1 := bitxor(v_s1, v_s0); g_128_state(0) := bitxor(bitxor(rotl(v_s0, 24), v_s1), shl64(v_s1, 16)); g_128_state(1) := rotl(v_s1, 37); RETURN v_result; END next128plus; FUNCTION next128plusplus RETURN uint64_t IS -- uint64_t next(void) { -- const uint64_t s0 = s[0]; -- uint64_t s1 = s[1]; -- const uint64_t result = rotl(s0 + s1, 17) + s0; -- -- s1 ^= s0; -- s[0] = rotl(s0, 49) ^ s1 ^ (s1 << 21); // a, b -- s[1] = rotl(s1, 28); // c -- -- return result; -- } v_s0 uint64_t := g_128_state(0); v_s1 uint64_t := g_128_state(1); v_result uint64_t := BITAND(rotl(BITAND(v_s0 + v_s1, c_64bitmask), 17) + v_s0, c_64bitmask); BEGIN v_s1 := bitxor(v_s1, v_s0); g_128_state(0) := bitxor(bitxor(rotl(v_s0, 49), v_s1), shl64(v_s1, 21)); g_128_state(1) := rotl(v_s1, 28); RETURN v_result; END next128plusplus; FUNCTION next128starstar RETURN uint64_t IS -- uint64_t next(void) { -- const uint64_t s0 = s[0]; -- uint64_t s1 = s[1]; -- const uint64_t result = rotl(s0 * 5, 7) * 9; -- -- s1 ^= s0; -- s[0] = rotl(s0, 24) ^ s1 ^ (s1 << 16); // a, b -- s[1] = rotl(s1, 37); // c -- -- return result; -- } -- v_s0 uint64_t := g_128_state(0); v_s1 uint64_t := g_128_state(1); v_result uint64_t := BITAND(rotl(BITAND(v_s0 * 5, c_64bitmask), 7) * 9, c_64bitmask); BEGIN v_s1 := bitxor(v_s1, v_s0); g_128_state(0) := bitxor(bitxor(rotl(v_s0, 24), v_s1), shl64(v_s1, 16)); g_128_state(1) := rotl(v_s1, 37); RETURN v_result; END next128starstar; FUNCTION next128 RETURN uint64_t IS BEGIN RETURN CASE g_mode WHEN c_128plus THEN next128plus WHEN c_128plusplus THEN next128plusplus WHEN c_128starstar THEN next128starstar END; END next128; PROCEDURE jump128 IS -- xoroshiro128+ and xoroshiro128** use the same jump array -- -- /* 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 uint64_t JUMP[] = { 0xdf900294d8f554a5, 0x170865df4b3201fc }; -- -- uint64_t s0 = 0; -- uint64_t s1 = 0; -- for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (JUMP[i] & UINT64_C(1) << b) { -- s0 ^= s[0]; -- s1 ^= s[1]; -- } -- next(); -- } -- -- s[0] = s0; -- s[1] = s1; -- } -- xoroshiro128++ uses a different jump array than the other two 128 bit modes -- -- static const uint64_t JUMP[] = { 0x2bd7a6a6e99c2ddc, 0x0992ccaf6a6fca05 }; v_jump state128_array_t := CASE WHEN g_mode IN (c_128plus, c_128starstar) THEN state128_array_t(0 => TO_NUMBER('df900294d8f554a5', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('170865df4b3201fc', 'xxxxxxxxxxxxxxxx')) WHEN g_mode = c_128plusplus THEN state128_array_t(0 => TO_NUMBER('2bd7a6a6e99c2ddc', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('0992ccaf6a6fca05', 'xxxxxxxxxxxxxxxx')) END; v_s0 uint64_t := 0; v_s1 uint64_t := 0; v_dummy uint64_t; BEGIN FOR i IN 0 .. 1 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_jump(i), shl64(1, b)) != 0 THEN v_s0 := bitxor(v_s0, g_128_state(0)); v_s1 := bitxor(v_s1, g_128_state(1)); END IF; v_dummy := next128; END LOOP; END LOOP; g_128_state(0) := v_s0; g_128_state(1) := v_s1; END jump128; PROCEDURE long_jump128 IS -- xoroshiro128+ and xoroshiro128** use the same long jump array -- /* 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 uint64_t LONG_JUMP[] = { 0xd2a98b26625eee7b, 0xdddf9b1090aa7ac1 }; -- -- uint64_t s0 = 0; -- uint64_t s1 = 0; -- for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (LONG_JUMP[i] & UINT64_C(1) << b) { -- s0 ^= s[0]; -- s1 ^= s[1]; -- } -- next(); -- } -- -- s[0] = s0; -- s[1] = s1; -- } -- xoroshiro128++ uses a different jump array than the other two 128 bit modes -- -- static const uint64_t LONG_JUMP[] = { 0x360fd5f2cf8d5d99, 0x9c6e6877736c46e3 }; v_longjump state128_array_t := CASE WHEN g_mode IN (c_128plus, c_128starstar) THEN state128_array_t(0 => TO_NUMBER('d2a98b26625eee7b', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('dddf9b1090aa7ac1', 'xxxxxxxxxxxxxxxx')) WHEN g_mode = c_128plusplus THEN state128_array_t(0 => TO_NUMBER('360fd5f2cf8d5d99', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('9c6e6877736c46e3', 'xxxxxxxxxxxxxxxx')) END; v_s0 uint64_t := 0; v_s1 uint64_t := 0; v_dummy uint64_t; BEGIN FOR i IN 0 .. 1 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_longjump(i), shl64(1, b)) != 0 THEN v_s0 := bitxor(v_s0, g_128_state(0)); v_s1 := bitxor(v_s1, g_128_state(1)); END IF; v_dummy := next128; END LOOP; END LOOP; g_128_state(0) := v_s0; g_128_state(1) := v_s1; END long_jump128; ----------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------- FUNCTION next256 RETURN uint64_t IS -- The result for the current state is calculated differently for each mode -- but the steps to increment the state for the next value are the same -- const uint64_t result = s[0] + s[3]; -- xoshiro256+ -- const uint64_t result = rotl(s[0] + s[3], 23) + s[0]; -- xoshiro256++ -- const uint64_t result = rotl(s[1] * 5, 7) * 9; -- xoshiro256** -- uint64_t next(void) { -- const uint64_t result = s[0] + s[3]; -- -- const uint64_t t = s[1] << 17; -- -- s[2] ^= s[0]; -- s[3] ^= s[1]; -- s[1] ^= s[2]; -- s[0] ^= s[3]; -- -- s[2] ^= t; -- -- s[3] = rotl(s[3], 45); -- -- return result; -- } v_result uint64_t; v_t uint64_t := shl64(g_256_state(1), 17); BEGIN v_result := BITAND( CASE g_mode WHEN c_256plus THEN g_256_state(0) + g_256_state(3) WHEN c_256plusplus THEN rotl(BITAND(g_256_state(0) + g_256_state(3), c_64bitmask), 23) + g_256_state(0) WHEN c_256starstar THEN rotl(BITAND(g_256_state(1) * 5, c_64bitmask), 7) * 9 END, c_64bitmask); g_256_state(2) := bitxor(g_256_state(2), g_256_state(0)); g_256_state(3) := bitxor(g_256_state(3), g_256_state(1)); g_256_state(1) := bitxor(g_256_state(1), g_256_state(2)); g_256_state(0) := bitxor(g_256_state(0), g_256_state(3)); g_256_state(2) := bitxor(g_256_state(2), v_t); g_256_state(3) := rotl(g_256_state(3), 45); RETURN v_result; END next256; PROCEDURE jump256 IS -- /* This is the jump function for the generator. It is equivalent -- to 2^128 calls to next(); it can be used to generate 2^128 -- non-overlapping subsequences for parallel computations. */ -- -- void jump(void) { -- static const uint64_t JUMP[] = { 0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c }; -- -- uint64_t s0 = 0; -- uint64_t s1 = 0; -- uint64_t s2 = 0; -- uint64_t s3 = 0; -- for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (JUMP[i] & UINT64_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 CONSTANT state256_array_t := state256_array_t(0 => TO_NUMBER('180ec6d33cfd0aba', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('d5a61266f0c9392c', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('a9582618e03fc9aa', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('39abdc4529b1661c', 'xxxxxxxxxxxxxxxx')) ; v_s0 uint64_t := 0; v_s1 uint64_t := 0; v_s2 uint64_t := 0; v_s3 uint64_t := 0; v_dummy uint64_t; BEGIN FOR i IN 0 .. 3 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_jump(i), shl64(1, b)) != 0 THEN v_s0 := bitxor(v_s0, g_256_state(0)); v_s1 := bitxor(v_s1, g_256_state(1)); v_s2 := bitxor(v_s2, g_256_state(2)); v_s3 := bitxor(v_s3, g_256_state(3)); END IF; v_dummy := next256; END LOOP; END LOOP; g_256_state(0) := v_s0; g_256_state(1) := v_s1; g_256_state(2) := v_s2; g_256_state(3) := v_s3; END jump256; PROCEDURE long_jump256 IS -- /* This is the long-jump function for the generator. It is equivalent to -- 2^192 calls to next(); it can be used to generate 2^64 starting points, -- from each of which jump() will generate 2^64 non-overlapping -- subsequences for parallel distributed computations. */ -- -- void long_jump(void) { -- static const uint64_t LONG_JUMP[] = { 0x76e15d3efefdcbbf, 0xc5004e441c522fb3, 0x77710069854ee241, 0x39109bb02acbe635 }; -- -- uint64_t s0 = 0; -- uint64_t s1 = 0; -- uint64_t s2 = 0; -- uint64_t s3 = 0; -- for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (LONG_JUMP[i] & UINT64_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_longjump CONSTANT state256_array_t := state256_array_t(0 => TO_NUMBER('76e15d3efefdcbbf', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('c5004e441c522fb3', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('77710069854ee241', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('39109bb02acbe635', 'xxxxxxxxxxxxxxxx')) ; v_s0 uint64_t := 0; v_s1 uint64_t := 0; v_s2 uint64_t := 0; v_s3 uint64_t := 0; v_dummy uint64_t; BEGIN FOR i IN 0 .. 3 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_longjump(i), shl64(1, b)) != 0 THEN v_s0 := bitxor(v_s0, g_256_state(0)); v_s1 := bitxor(v_s1, g_256_state(1)); v_s2 := bitxor(v_s2, g_256_state(2)); v_s3 := bitxor(v_s3, g_256_state(3)); END IF; v_dummy := next256; END LOOP; END LOOP; g_256_state(0) := v_s0; g_256_state(1) := v_s1; g_256_state(2) := v_s2; g_256_state(3) := v_s3; END long_jump256; ----------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------- FUNCTION next512 RETURN uint64_t IS -- The result for the current state is calculated differently for each mode -- but the steps to increment the state for the next value are the same -- const uint64_t result = s[0] + s[2]; -- xoshiro512+ -- const uint64_t result = rotl(s[0] + s[2], 17) + s[2]; -- xoshiro512++ -- const uint64_t result = rotl(s[1] * 5, 7) * 9; -- xoshiro512** -- uint64_t next(void) { -- const uint64_t result = s[0] + s[2]; -- -- const uint64_t t = s[1] << 11; -- -- s[2] ^= s[0]; -- s[5] ^= s[1]; -- s[1] ^= s[2]; -- s[7] ^= s[3]; -- s[3] ^= s[4]; -- s[4] ^= s[5]; -- s[0] ^= s[6]; -- s[6] ^= s[7]; -- -- s[6] ^= t; -- -- s[7] = rotl(s[7], 21); -- -- return result; -- } v_result uint64_t := BITAND( CASE g_mode WHEN c_512plus THEN g_512_state(0) + g_512_state(2) WHEN c_512plusplus THEN rotl(BITAND(g_512_state(0) + g_512_state(2), c_64bitmask), 17) + g_512_state(2) WHEN c_512starstar THEN rotl(BITAND(g_512_state(1) * 5, c_64bitmask), 7) * 9 END, c_64bitmask); v_t uint64_t := shl64(g_512_state(1), 11); BEGIN g_512_state(2) := bitxor(g_512_state(2), g_512_state(0)); g_512_state(5) := bitxor(g_512_state(5), g_512_state(1)); g_512_state(1) := bitxor(g_512_state(1), g_512_state(2)); g_512_state(7) := bitxor(g_512_state(7), g_512_state(3)); g_512_state(3) := bitxor(g_512_state(3), g_512_state(4)); g_512_state(4) := bitxor(g_512_state(4), g_512_state(5)); g_512_state(0) := bitxor(g_512_state(0), g_512_state(6)); g_512_state(6) := bitxor(g_512_state(6), g_512_state(7)); g_512_state(6) := bitxor(g_512_state(6), v_t); g_512_state(7) := rotl(g_512_state(7), 21); RETURN v_result; END next512; PROCEDURE jump512 IS -- /* This is the jump function for the generator. It is equivalent -- to 2^256 calls to next(); it can be used to generate 2^256 -- non-overlapping subsequences for parallel computations. */ -- -- void jump(void) { -- static const uint64_t JUMP[] = { 0x33ed89b6e7a353f9, 0x760083d7955323be, 0x2837f2fbb5f22fae, 0x4b8c5674d309511c, 0xb11ac47a7ba28c25, 0xf1be7667092bcc1c, 0x53851efdb6df0aaf, 0x1ebbc8b23eaf25db }; -- -- uint64_t t[sizeof s / sizeof *s]; -- memset(t, 0, sizeof t); -- for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (JUMP[i] & UINT64_C(1) << b) -- for(int w = 0; w < sizeof s / sizeof *s; w++) -- t[w] ^= s[w]; -- next(); -- } -- -- memcpy(s, t, sizeof s); -- } v_jump CONSTANT state512_array_t := state512_array_t(0 => TO_NUMBER('33ed89b6e7a353f9', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('760083d7955323be', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('2837f2fbb5f22fae', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('4b8c5674d309511c', 'xxxxxxxxxxxxxxxx'), 4 => TO_NUMBER('b11ac47a7ba28c25', 'xxxxxxxxxxxxxxxx'), 5 => TO_NUMBER('f1be7667092bcc1c', 'xxxxxxxxxxxxxxxx'), 6 => TO_NUMBER('53851efdb6df0aaf', 'xxxxxxxxxxxxxxxx'), 7 => TO_NUMBER('1ebbc8b23eaf25db', 'xxxxxxxxxxxxxxxx')) ; v_t state512_array_t := state512_array_t(0 => 0, 1 => 0, 2 => 0, 3 => 0, 4 => 0, 5 => 0, 6 => 0, 7 => 0); v_dummy uint64_t; BEGIN FOR i IN 0 .. 7 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_jump(i), shl64(1, b)) != 0 THEN FOR w IN 0 .. 7 LOOP v_t(w) := bitxor(v_t(w), g_512_state(w)); END LOOP; END IF; v_dummy := next512; END LOOP; END LOOP; FOR i IN 0 .. 7 LOOP g_512_state(i) := v_t(i); END LOOP; END jump512; PROCEDURE long_jump512 IS -- /* This is the long-jump function for the generator. It is equivalent to -- 2^384 calls to next(); it can be used to generate 2^128 starting points, -- from each of which jump() will generate 2^128 non-overlapping -- subsequences for parallel distributed computations. */ -- -- void long_jump(void) { -- static const uint64_t LONG_JUMP[] = { 0x11467fef8f921d28, 0xa2a819f2e79c8ea8, 0xa8299fc284b3959a, 0xb4d347340ca63ee1, 0x1cb0940bedbff6ce, 0xd956c5c4fa1f8e17, 0x915e38fd4eda93bc, 0x5b3ccdfa5d7daca5 }; -- -- uint64_t t[sizeof s / sizeof *s]; -- memset(t, 0, sizeof t); -- for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (LONG_JUMP[i] & UINT64_C(1) << b) -- for(int w = 0; w < sizeof s / sizeof *s; w++) -- t[w] ^= s[w]; -- next(); -- } -- -- memcpy(s, t, sizeof s); -- } -- v_longjump CONSTANT state512_array_t := state512_array_t(0 => TO_NUMBER('11467fef8f921d28', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('a2a819f2e79c8ea8', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('a8299fc284b3959a', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('b4d347340ca63ee1', 'xxxxxxxxxxxxxxxx'), 4 => TO_NUMBER('1cb0940bedbff6ce', 'xxxxxxxxxxxxxxxx'), 5 => TO_NUMBER('d956c5c4fa1f8e17', 'xxxxxxxxxxxxxxxx'), 6 => TO_NUMBER('915e38fd4eda93bc', 'xxxxxxxxxxxxxxxx'), 7 => TO_NUMBER('5b3ccdfa5d7daca5', 'xxxxxxxxxxxxxxxx')) ; v_t state512_array_t := state512_array_t(0 => 0, 1 => 0, 2 => 0, 3 => 0, 4 => 0, 5 => 0, 6 => 0, 7 => 0); v_dummy uint64_t; BEGIN FOR i IN 0 .. 7 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_longjump(i), shl64(1, b)) != 0 THEN FOR w IN 0 .. 7 LOOP v_t(w) := bitxor(v_t(w), g_512_state(w)); END LOOP; END IF; v_dummy := next512; END LOOP; END LOOP; FOR i IN 0 .. 7 LOOP g_512_state(i) := v_t(i); END LOOP; END long_jump512; ----------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------- FUNCTION next1024 RETURN uint64_t IS -- The result for the current state is calculated differently for each mode -- but the steps to increment the state for the next value are the same -- const uint64_t result = s0 * 0x9e3779b97f4a7c13; -- xoshiro1024* -- const uint64_t result = rotl(s0 + s15, 23) + s15; -- xoshiro1024++ -- cconst uint64_t result = rotl(s0 * 5, 7) * 9; -- xoshiro1024** -- uint64_t next(void) { -- const int q = p; -- const uint64_t s0 = s[p = (p + 1) & 15]; -- uint64_t s15 = s[q]; -- const uint64_t result = s0 * 0x9e3779b97f4a7c13; -- -- s15 ^= s0; -- s[q] = rotl(s0, 25) ^ s15 ^ (s15 << 27); -- s[p] = rotl(s15, 36); -- -- return result; -- } v_q INTEGER; v_s0 uint64_t; v_s15 uint64_t; v_result uint64_t; BEGIN v_q := g_1024_p; g_1024_p := BITAND(g_1024_p + 1, 15); v_s0 := g_1024_state(g_1024_p); v_s15 := g_1024_state(v_q); v_result := CASE g_mode WHEN c_1024star THEN big_bitand(v_s0 * TO_NUMBER('9e3779b97f4a7c13', 'xxxxxxxxxxxxxxxx'), c_64bitmask) WHEN c_1024plusplus THEN BITAND(rotl(BITAND(v_s0 + v_s15, c_64bitmask), 23) + v_s15, c_64bitmask) WHEN c_1024starstar THEN BITAND(rotl(BITAND(v_s0 * 5, c_64bitmask), 7) * 9, c_64bitmask) END; v_s15 := bitxor(v_s15, v_s0); g_1024_state(v_q) := bitxor(bitxor(rotl(v_s0, 25), v_s15), shl64(v_s15, 27)); g_1024_state(g_1024_p) := rotl(v_s15, 36); RETURN v_result; END next1024; PROCEDURE jump1024 IS -- /* This is the jump function for the generator. It is equivalent -- to 2^512 calls to next(); it can be used to generate 2^512 -- non-overlapping subsequences for parallel computations. */ -- -- void jump() { -- static const uint64_t JUMP[] = { 0x931197d8e3177f17, -- 0xb59422e0b9138c5f, 0xf06a6afb49d668bb, 0xacb8a6412c8a1401, -- 0x12304ec85f0b3468, 0xb7dfe7079209891e, 0x405b7eec77d9eb14, -- 0x34ead68280c44e4a, 0xe0e4ba3e0ac9e366, 0x8f46eda8348905b7, -- 0x328bf4dbad90d6ff, 0xc8fd6fb31c9effc3, 0xe899d452d4b67652, -- 0x45f387286ade3205, 0x03864f454a8920bd, 0xa68fa28725b1b384 }; -- -- uint64_t t[sizeof s / sizeof *s]; -- memset(t, 0, sizeof t); -- for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (JUMP[i] & UINT64_C(1) << b) -- for(int j = 0; j < sizeof s / sizeof *s; j++) -- t[j] ^= s[(j + p) & sizeof s / sizeof *s - 1]; -- next(); -- } -- -- for(int i = 0; i < sizeof s / sizeof *s; i++) { -- s[(i + p) & sizeof s / sizeof *s - 1] = t[i]; -- } -- } v_jump CONSTANT state1024_array_t := state1024_array_t(0 => TO_NUMBER('931197d8e3177f17', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('b59422e0b9138c5f', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('f06a6afb49d668bb', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('acb8a6412c8a1401', 'xxxxxxxxxxxxxxxx'), 4 => TO_NUMBER('12304ec85f0b3468', 'xxxxxxxxxxxxxxxx'), 5 => TO_NUMBER('b7dfe7079209891e', 'xxxxxxxxxxxxxxxx'), 6 => TO_NUMBER('405b7eec77d9eb14', 'xxxxxxxxxxxxxxxx'), 7 => TO_NUMBER('34ead68280c44e4a', 'xxxxxxxxxxxxxxxx'), 8 => TO_NUMBER('e0e4ba3e0ac9e366', 'xxxxxxxxxxxxxxxx'), 9 => TO_NUMBER('8f46eda8348905b7', 'xxxxxxxxxxxxxxxx'), 10 => TO_NUMBER('328bf4dbad90d6ff', 'xxxxxxxxxxxxxxxx'), 11 => TO_NUMBER('c8fd6fb31c9effc3', 'xxxxxxxxxxxxxxxx'), 12 => TO_NUMBER('e899d452d4b67652', 'xxxxxxxxxxxxxxxx'), 13 => TO_NUMBER('45f387286ade3205', 'xxxxxxxxxxxxxxxx'), 14 => TO_NUMBER('03864f454a8920bd', 'xxxxxxxxxxxxxxxx'), 15 => TO_NUMBER('a68fa28725b1b384', 'xxxxxxxxxxxxxxxx')) ; v_t state1024_array_t := state1024_array_t(0 => 0, 1 => 0, 2 => 0, 3 => 0, 4 => 0, 5 => 0, 6 => 0, 7 => 0, 8 => 0, 9 => 0, 10 => 0, 11 => 0, 12 => 0, 13 => 0, 14 => 0, 15 => 0); v_dummy uint64_t; BEGIN FOR i IN 0 .. 15 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_jump(i), shl64(1, b)) != 0 THEN FOR j IN 0 .. 15 LOOP v_t(j) := bitxor(v_t(j), g_1024_state(BITAND(j + g_1024_p, 15))); END LOOP; END IF; v_dummy := next1024; END LOOP; END LOOP; FOR i IN 0 .. 15 LOOP g_1024_state(BITAND(i + g_1024_p, 15)) := v_t(i); END LOOP; END jump1024; PROCEDURE long_jump1024 IS -- /* This is the long-jump function for the generator. It is equivalent to -- 2^768 calls to next(); it can be used to generate 2^256 starting points, -- from each of which jump() will generate 2^256 non-overlapping -- subsequences for parallel distributed computations. */ -- -- void long_jump(void) { -- static const uint64_t LONG_JUMP[] = { 0x7374156360bbf00f, -- 0x4630c2efa3b3c1f6, 0x6654183a892786b1, 0x94f7bfcbfb0f1661, -- 0x27d8243d3d13eb2d, 0x9701730f3dfb300f, 0x2f293baae6f604ad, -- 0xa661831cb60cd8b6, 0x68280c77d9fe008c, 0x50554160f5ba9459, -- 0x2fc20b17ec7b2a9a, 0x49189bbdc8ec9f8f, 0x92a65bca41852cc1, -- 0xf46820dd0509c12a, 0x52b00c35fbf92185, 0x1e5b3b7f589e03c1 }; -- -- uint64_t t[sizeof s / sizeof *s]; -- memset(t, 0, sizeof t); -- for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++) -- for(int b = 0; b < 64; b++) { -- if (LONG_JUMP[i] & UINT64_C(1) << b) -- for(int j = 0; j < sizeof s / sizeof *s; j++) -- t[j] ^= s[(j + p) & sizeof s / sizeof *s - 1]; -- next(); -- } -- -- for(int i = 0; i < sizeof s / sizeof *s; i++) { -- s[(i + p) & sizeof s / sizeof *s - 1] = t[i]; -- } -- } v_longjump CONSTANT state1024_array_t := state1024_array_t(0 => TO_NUMBER('7374156360bbf00f', 'xxxxxxxxxxxxxxxx'), 1 => TO_NUMBER('4630c2efa3b3c1f6', 'xxxxxxxxxxxxxxxx'), 2 => TO_NUMBER('6654183a892786b1', 'xxxxxxxxxxxxxxxx'), 3 => TO_NUMBER('94f7bfcbfb0f1661', 'xxxxxxxxxxxxxxxx'), 4 => TO_NUMBER('27d8243d3d13eb2d', 'xxxxxxxxxxxxxxxx'), 5 => TO_NUMBER('9701730f3dfb300f', 'xxxxxxxxxxxxxxxx'), 6 => TO_NUMBER('2f293baae6f604ad', 'xxxxxxxxxxxxxxxx'), 7 => TO_NUMBER('a661831cb60cd8b6', 'xxxxxxxxxxxxxxxx'), 8 => TO_NUMBER('68280c77d9fe008c', 'xxxxxxxxxxxxxxxx'), 9 => TO_NUMBER('50554160f5ba9459', 'xxxxxxxxxxxxxxxx'), 10 => TO_NUMBER('2fc20b17ec7b2a9a', 'xxxxxxxxxxxxxxxx'), 11 => TO_NUMBER('49189bbdc8ec9f8f', 'xxxxxxxxxxxxxxxx'), 12 => TO_NUMBER('92a65bca41852cc1', 'xxxxxxxxxxxxxxxx'), 13 => TO_NUMBER('f46820dd0509c12a', 'xxxxxxxxxxxxxxxx'), 14 => TO_NUMBER('52b00c35fbf92185', 'xxxxxxxxxxxxxxxx'), 15 => TO_NUMBER('1e5b3b7f589e03c1', 'xxxxxxxxxxxxxxxx')) ; v_t state1024_array_t := state1024_array_t(0 => 0, 1 => 0, 2 => 0, 3 => 0, 4 => 0, 5 => 0, 6 => 0, 7 => 0, 8 => 0, 9 => 0, 10 => 0, 11 => 0, 12 => 0, 13 => 0, 14 => 0, 15 => 0); v_dummy uint64_t; BEGIN FOR i IN 0 .. 15 LOOP FOR b IN 0 .. 63 LOOP IF BITAND(v_longjump(i), shl64(1, b)) != 0 THEN FOR j IN 0 .. 15 LOOP v_t(j) := bitxor(v_t(j), g_1024_state(BITAND(j + g_1024_p, 15))); END LOOP; END IF; v_dummy := next1024; END LOOP; END LOOP; FOR i IN 0 .. 15 LOOP g_1024_state(BITAND(i + g_1024_p, 15)) := v_t(i); END LOOP; END long_jump1024; ----------------------------------------------------------------------------------------------------------------- ----------------------------------------------------------------------------------------------------------------- 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 'xoroshiro128+' THEN c_128plus WHEN 'xoroshiro128++' THEN c_128plusplus WHEN 'xoroshiro128**' THEN c_128starstar WHEN 'xoshiro256+' THEN c_256plus WHEN 'xoshiro256++' THEN c_256plusplus WHEN 'xoshiro256**' THEN c_256starstar WHEN 'xoshiro512+' THEN c_512plus WHEN 'xoshiro512++' THEN c_512plusplus WHEN 'xoshiro512**' THEN c_512starstar WHEN 'xoroshiro1024*' THEN c_1024star WHEN 'xoroshiro1024++' THEN c_1024plusplus WHEN 'xoroshiro1024**' THEN c_1024starstar 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_128plus THEN 'xoroshiro128+' WHEN c_128plusplus THEN 'xoroshiro128++' WHEN c_128starstar THEN 'xoroshiro128**' WHEN c_256plus THEN 'xoshiro256+' WHEN c_256plusplus THEN 'xoshiro256++' WHEN c_256starstar THEN 'xoshiro256**' WHEN c_512plus THEN 'xoshiro512+' WHEN c_512plusplus THEN 'xoshiro512++' WHEN c_512starstar THEN 'xoshiro512**' WHEN c_1024star THEN 'xoroshiro1024*' WHEN c_1024plusplus THEN 'xoroshiro1024++' WHEN c_1024starstar THEN 'xoroshiro1024**' END; END get_mode_name; FUNCTION get_next RETURN uint64_t IS BEGIN RETURN CASE WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN next128 WHEN g_mode IN (c_256plus, c_256plusplus, c_256starstar) THEN next256 WHEN g_mode IN (c_512plus, c_512plusplus, c_512starstar) THEN next512 WHEN g_mode IN (c_1024star, c_1024plusplus, c_1024starstar) THEN next1024 END; END get_next; PROCEDURE jump IS BEGIN CASE WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN jump128; WHEN g_mode IN (c_256plus, c_256plusplus, c_256starstar) THEN jump256; WHEN g_mode IN (c_512plus, c_512plusplus, c_512starstar) THEN jump512; WHEN g_mode IN (c_1024star, c_1024plusplus, c_1024starstar) THEN jump1024; ELSE RAISE invalid_mode; END CASE; END jump; PROCEDURE long_jump IS BEGIN CASE WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN long_jump128; WHEN g_mode IN (c_256plus, c_256plusplus, c_256starstar) THEN long_jump256; WHEN g_mode IN (c_512plus, c_512plusplus, c_512starstar) THEN long_jump512; WHEN g_mode IN (c_1024star, c_1024plusplus, c_1024starstar) THEN long_jump1024; ELSE RAISE invalid_mode; END CASE; END long_jump; PROCEDURE set_state(p_seed IN seed_tab) IS BEGIN CASE WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN g_128_state := state128_array_t(); FOR i IN 0 .. 1 LOOP g_128_state(i) := p_seed(i + 1); END LOOP; WHEN g_mode IN (c_256plus, c_256plusplus, c_256starstar) THEN g_256_state := state256_array_t(); FOR i IN 0 .. 3 LOOP g_256_state(i) := p_seed(i + 1); END LOOP; WHEN g_mode IN (c_512plus, c_512plusplus, c_512starstar) THEN g_512_state := state512_array_t(); FOR i IN 0 .. 7 LOOP g_512_state(i) := p_seed(i + 1); END LOOP; WHEN g_mode IN (c_1024star, c_1024plusplus, c_1024starstar) THEN g_1024_p := 0; g_1024_state := state1024_array_t(); FOR i IN 0 .. 15 LOOP g_1024_state(i) := p_seed(i + 1); END LOOP; ELSE RAISE invalid_mode; END CASE; END set_state; PROCEDURE splitmix64_set_state(p_seed uint64_t) IS BEGIN g_splitmix64_state := p_seed; END splitmix64_set_state; FUNCTION splitmix64_next RETURN uint64_t IS /* Based on the splitmix64.c code by Sebastiano Vigna https://xorshift.di.unimi.it/splitmix64.c static uint64_t x; uint64_t next() { uint64_t z = (x += 0x9e3779b97f4a7c15); z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9; z = (z ^ (z >> 27)) * 0x94d049bb133111eb; return z ^ (z >> 31); } This will not be nearly as fast as the original as it uses division instead of bit shifting. However, the functionality of it should be consistent */ v_z NUMBER; BEGIN g_splitmix64_state := BITAND(g_splitmix64_state + TO_NUMBER('9e3779b97f4a7c15', 'xxxxxxxxxxxxxxxx'), c_64bitmask); v_z := g_splitmix64_state; v_z := bitxor(v_z, shr64(v_z, 30)); v_z := v_z * TO_NUMBER('bf58476d1ce4e5b9', 'xxxxxxxxxxxxxxxx'); v_z := big_bitand(v_z, c_64bitmask); v_z := bitxor(v_z, shr64(v_z, 27)); v_z := v_z * TO_NUMBER('94d049bb133111eb', 'xxxxxxxxxxxxxxxx'); v_z := big_bitand(v_z, c_64bitmask); RETURN BITAND(bitxor(v_z, shr64(v_z, 31)), c_64bitmask); END splitmix64_next; END;
The BIG_BITAND function I described in detail previously is used within the shl64, next1024, and splitmix64_next functions to ensure large (128-bit) values are handled correctly.
I hope you find the code interesting and useful. I have one more article planned in this series – covering the the 32-bit variations within this generator family. While they could have been combined into a single package; I felt it was a better solution to keep them separate as they are unlikely to be used together in practical application. In most cases I expect users will want results of a single single bit-size, 64 or 32. Combining them wouldn’t serve much purpose.
Questions and comments, as always, are welcome.