Close

PL/SQL package for 64-bit xoshiro/xoroshiro pseudorandom number generators

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:

  1. Choose the algorithm mode (the package defaults to xoshiro256++)
  2. Seed the state array
  3. 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.

Leave a Reply