Перевод публикуется с полученного мною личного разрешения автора, Sean D. Stuber.
Оригинальная статья: "PL/SQL package for 64-bit xoshiro/xoroshiro pseudorandom number generators" , 01.03.2025
Приведённый ниже пакет поддерживает все алгоритмы xoshiro/xoroshiro, возвращающие 64-битные значения.
- xoroshiro128+
- xoroshiro128++
- xoroshiro128**
- xoshiro256+
- xoshiro256++
- xoshiro256**
- xoshiro512+
- xoshiro512++
- xoshiro512**
- xoroshiro1024*
- xoroshiro1024++
- xoroshiro1024**
Цифры в названиях указывают на размер массива состояний, с которым работает каждый метод. Использование каждого из них осуществляется по одной и той же схеме:
- Выберите режим алгоритма (по умолчанию в пакете используется xoshiro256++)
- Заполнение массива состояний
- Извлекайте следующее число из выбранного алгоритма столько раз, сколько необходимо
Например
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При необходимости вы можете вызвать функции jump или long_jump, чтобы перескочить вперёд в состоянии, имитируя множество (от 2^64 до 2^768 в зависимости от режима и размера перескока) вызовов функции get_next. Как отмечают Блэкман и Винья, они обычно используются для создания непересекающихся последовательностей для параллельной обработки. При использовании в одном потоке в них нет необходимости.
В дополнение к ним я также включил функцию splitmix64, которая сама по себе является генератором псевдослучайных чисел и может использоваться аналогичным образом. Задайте начальное значение, а затем вызовите splitmix64_next, чтобы получить следующее число.
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Однако это не типичное использование. Скорее, splitmix64 предлагается в качестве средства инициализации массивов состояний других алгоритмов. Ниже я устанавливаю состояние splitmix64, а затем создаю начальное значение из 4 значений для генератора xoshiro256**.
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Ниже приведен полный пакет для всех 64-разрядных алгоритмов.
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
-- This email address is being protected from spambots. You need JavaScript enabled to view it.
--
-- 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++, xoroshi-ro1024**
-- 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 val-ue.
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;Функция BIG_BITAND, которую я подробно описал ранее, используется в функциях shl64, next1024 и splitmix64_next для корректной обработки больших (128-битных) значений.
Надеюсь, вам понравится этот код и он будет вам полезен. В этой серии я планирую ещё одну статью, посвящённую 32-битным вариациям этого семейства генераторов. Хотя их можно было бы объединить в один пакет, я решил, что лучше оставить их отдельно, так как они вряд ли будут использоваться вместе в практических целях. В большинстве случаев, я думаю, пользователям понадобятся результаты с одним размером битов, 64 или 32. Их объединение не принесёт особой пользы.
Пакет PL/SQL для 32-разрядных генераторов псевдослучайных чисел xoshiro/xoroshiro
PL/SQL package for 32-bit xoshiro/xoroshiro pseudorandom number generators
Завершая серию статей об алгоритмах xoshiro/xoroshiro, в этой статье я представляю пакет для возврата 32-битных значений с использованием следующих вариантов.
- xoroshiro64*
- xoroshiro64**
- xoshiro128+
- xoshiro128++
- xoshiro128**
Цифры в названиях указывают на размер массива состояний, с которым работает каждый метод. 64-битные методы работают с двумя 32-битными значениями, а 128-битные — с четырьмя 32-битными значениями.
Использование каждого из них осуществляется по одной и той же схеме:
- Выберите режим алгоритма (в пакете по умолчанию используется значение xoshiro64*).
- Заполнение массива состояний
- Извлекайте следующее число из выбранного алгоритма столько раз, сколько необходимо
Например
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В то время как все 64-битные варианты поддерживают функции jump и long_jump, только те, у которых есть 128-битные массивы состояний, поддерживают их в 32-битных вариантах. В состоянии, имитируя множество (от 2^64 до 2^768 в зависимости от режима и размера перехода) вызовов функции get_next. Как отмечают Блэкман и Винья, они обычно используются для создания непересекающихся последовательностей для параллельной обработки. При однопоточном использовании в них не было бы необходимости.
64-битные варианты включали splitmix64, который я портировал из оригинальной реализации Vigna на языке C. Для 32-битных вариантов не было предусмотрено соответствующего метода, поэтому я вернулся к оригинальной статье Стила, Ли и Флуда. Следуя логике splitmix64, я также использовал 32-битную версию MurmurHash3 от Эпплби.
|
Fast Splittable Pseudorandom Number Generators https://gee.cs.oswego.edu/dl/papers/oopsla14.pdf Guy L. Steele Jr. Oracle Labs Doug Lea SUNY Oswego Christine H. Flood Red Hat Inc |
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Как и в случае с splitmix64, прямое использование splitmix32 для генерации псевдослучайных чисел не является типичным применением. Скорее, splitmix32 лучше использовать как средство инициализации массивов состояний других алгоритмов. Ниже я устанавливаю состояние splitmix32, а затем создаю начальное значение из 4 значений для генератора xoshiro128**.
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Ниже представлен полный пакет для всех 32-битных алгоритмов.
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
-- This email address is being protected from spambots. You need JavaScript enabled to view it.
--
-- 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;
-- 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 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;
/
Поскольку начальные числа и промежуточные результаты меньше, эти варианты алгоритма могут использовать встроенную функцию BITAND и не требуют применения функции BIG_BITAND для очень больших значений.
На этом я завершаю серию статей и надеюсь, что вы найдёте код интересным и полезным. Как я упоминал в своей статье о 64-битных вычислениях, я намеренно разделил 32-битные и 64-битные алгоритмы на отдельные пакеты, предполагая, что они вряд ли будут использоваться вместе. Но если вы решите их объединить, я определил константы режима таким образом, чтобы они были взаимоисключающими и их можно было комбинировать, не опасаясь неоднозначности из-за совпадающих значений.
Вопросы и комментарии, как всегда, приветствуются.