Close

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

Wrapping up my series on the, xoshiro/xoroshiro algorithms, in this article I present a package to return 32-bit values using the following variants.

  • xoroshiro64*
  • xoroshiro64**
  • xoshiro128+
  • xoshiro128++
  • xoshiro128**

The numbers in the names refer to the size of the state array that each method operates on. The 64s operate on two 32-bit values, while the 128s operate on four 32-bit values.

Using each of them follows the same pattern:

  1. Choose the algorithm mode (the package defaults to xoshiro64*)
  2. Seed the state array
  3. Retrieve the next number from the chosen algorithm as many times as needed

For example…

SQL> BEGIN
  2      xoshiro32.set_mode('xoroshiro64**');
  3      xoshiro32.set_state(xoshiro32.seed_tab(12345, 67890));
  4
  5      FOR i IN 0 .. 9
  6      LOOP
  7          DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.get_next, 10));
  8      END LOOP;
  9  END;
 10  /
0 3157960260
1 4142509522
2 1831851427
3  506054173
4 2910589752
5 1819521659
6 3282141937
7 2257682835
8 2133372007
9 3757018772

While all of the 64-bit variants supported jump and long_jump functions, only those with 128-bit state arrays support them in the 32-bit variants. in the state by simulating many (from 2^64 to 2^768 depending on the mode and jump size) calls of the get_next function. As noted by Blackman and Vigna, these are typically used to generate non-overlapping sequences for parallel processing. In single threaded use there would be no need for them.

The 64-bit variants included splitmix64, which I ported from Vigna’s original c implementation. No corresponding method was provided for the 32-bit variants, so I went back to the original SplitMix paper by Steele, Lea, and Flood. Mirroring the logic of splitmix64, I also used a 32-bit version of the MurmurHash3 by Appleby to calculate the Golden Gamma.

SQL> BEGIN
  2      xoshiro32.splitmix32_set_state(12345);
  3
  4      FOR i IN 0 .. 9
  5      LOOP
  6          DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.splitmix32_next, 10));
  7      END LOOP;
  8  END;
  9  /
0 1200724404
1  818072533
2  996137225
3 2397394836
4 4079075752
5 2274189806
6 2795887828
7 4161515127
8 3291005408
9  722528451

Just as with splitmix64, using splitmix32 directly to generate pseudorandom numbers is not the typical usage . Rather, splitmix32 is best used as a means of seeding the state arrays of the other algorithms. Below, I set the splitmix32 state, and then create a seed of 4 values for the xoshiro128** generator.

SQL> DECLARE
  2      v_seed   xoshiro32.seed_tab := xoshiro32.seed_tab();
  3  BEGIN
  4      xoshiro32.splitmix32_set_state(12345);
  5      v_seed.EXTEND(4);
  6
  7      FOR i IN 1 .. 4
  8      LOOP
  9          v_seed(i) := xoshiro32.splitmix32_next;
 10      END LOOP;
 11
 12      xoshiro32.set_mode('xoshiro128**');
 13      xoshiro32.set_state(v_seed);
 14
 15      FOR i IN 0 .. 9
 16      LOOP
 17          DBMS_OUTPUT.put_line(i || ' ' || LPAD(xoshiro32.get_next, 10));
 18      END LOOP;
 19  END;
 20  /
0  518667457
1  440444462
2 4232892992
3 3757857622
4 3939018813
5 1334683535
6 3795058715
7 2092637810
8 2829112157
9  779180383

Following is the complete package for all of the 32-bit algorithms.

CREATE OR REPLACE PACKAGE xoshiro32
IS
    -- Based on xoshiro/xoroshiro family of pseudorandom number generators
    --  by by David Blackman and Sebastiano Vigna
    -- original  https://prng.di.unimi.it

    --                    .///.
    --                   (0 o)
    ---------------0000--(_)--0000---------------
    --
    --  Sean D. Stuber
    --  sean.stuber@gmail.com
    --
    --             oooO      Oooo
    --------------(   )-----(   )---------------
    --             \ (       ) /
    --              \_)     (_/

    -- It would be nice to define an explicit RANGE 0..4294967295 for the subtype
    -- but as of 23ai, only subtypes of PLS_INTEGER may have a range.
    SUBTYPE uint32_t IS NUMBER(10, 0);

    TYPE seed_tab IS TABLE OF uint32_t;

    SUBTYPE bits_t IS SIMPLE_INTEGER RANGE 0 .. 31;

    SUBTYPE mode_t IS SIMPLE_INTEGER RANGE 13 .. 17;

    -- These do not support jump and long_jump operations
    c_64star        CONSTANT mode_t := 13; -- xoroshiro64*
    c_64starstar    CONSTANT mode_t := 14; -- xoroshiro64**

    -- These do support jump and long_jump operations    
    c_128plus       CONSTANT mode_t := 15; -- xoshiro128+
    c_128plusplus   CONSTANT mode_t := 16; -- xoshiro128++
    c_128starstar   CONSTANT mode_t := 17; -- xoshiro128**

    --Exceptions
    invalid_mode             EXCEPTION;

    -- Select a mode by numeric value or by name
    PROCEDURE set_mode(p_mode IN mode_t);

    PROCEDURE set_mode(p_mode IN VARCHAR2);

    -- Return the current mode number or name
    FUNCTION get_mode
        RETURN mode_t;

    FUNCTION get_mode_name
        RETURN VARCHAR2;

    -- Given a set of seed values, initialize the state array for the current mode size
    -- Different modes for a particular size can share the same state array.
    PROCEDURE set_state(p_seed seed_tab);

    -- Return the next pseudorandom generator value for the current mode
    FUNCTION get_next
        RETURN uint32_t;

    -- Simulate many calls to "get_next".
    -- Typically used to create non-overlapping sequences from the same initial state.
    -- xoroshiro64* and xoroshiro64** do not support jump or long_jump operations
    PROCEDURE jump;

    -- Simulate even more calls to "get_next".
    --Also used to create non-overlapping sequences from the same initial state.
    -- xoroshiro64* and xoroshiro64** do not support jump or long_jump operations
    PROCEDURE long_jump;

    -- SplitMix32 isn't part of the xoshiro/xoroshiro family
    -- However using a function like it to fill the state arrays instead of a simple
    -- integer is recommended by the original authors.
    PROCEDURE splitmix32_set_state(p_seed uint32_t);

    FUNCTION splitmix32_next
        RETURN uint32_t;
END;
/

As a comparison reference for myself and readers, I included snippets of the original c code within the comments of the package body for the various procedures and functions of each mode.

CREATE OR REPLACE PACKAGE BODY xoshiro32
IS
    c_32bitmask   CONSTANT INTEGER := POWER(2, 32) - 1;

    SUBTYPE array_64index_t IS PLS_INTEGER RANGE 0 .. 1;

    SUBTYPE array_128index_t IS PLS_INTEGER RANGE 0 .. 3;

    TYPE state64_array_t IS TABLE OF uint32_t
        INDEX BY array_64index_t;

    TYPE state128_array_t IS TABLE OF uint32_t
        INDEX BY array_128index_t;

    g_mode                 mode_t := c_64star;

    g_splitmix32_state     uint32_t;

    g_64_state             state64_array_t;
    g_128_state            state128_array_t;

    FUNCTION bitor32(p_a IN uint32_t, p_b IN uint32_t)
        RETURN uint32_t
    IS
    BEGIN
        RETURN BITAND(p_a + p_b - BITAND(p_a, p_b), c_32bitmask);
    END bitor32;

    FUNCTION bitxor32(p_a IN uint32_t, p_b IN uint32_t)
        RETURN uint32_t
    IS
    BEGIN
        RETURN BITAND((p_a + p_b) - (BITAND(p_a, p_b) * 2), c_32bitmask);
    END bitxor32;

    FUNCTION shl32(p_n IN uint32_t, p_bits IN bits_t)
        RETURN uint32_t
    IS
    BEGIN
        RETURN CASE p_bits WHEN 0 THEN p_n ELSE BITAND(p_n * (2 ** p_bits), c_32bitmask) END;
    END shl32;

    FUNCTION shr32(p_n IN uint32_t, p_bits IN bits_t)
        RETURN uint32_t
    IS
    BEGIN
        RETURN CASE p_bits WHEN 0 THEN p_n ELSE BITAND(FLOOR(p_n / (2 ** p_bits)), c_32bitmask) END;
    END shr32;

    FUNCTION rotl32(p_x uint32_t, p_k IN bits_t)
        RETURN uint32_t
    IS
    --  static inline uint32_t rotl(const uint32_t x, int k) {
    --   return (x << k) | (x >> (32 - k));
    --  }
    BEGIN
        RETURN CASE WHEN p_k = 0 THEN p_x ELSE bitor32(shl32(p_x, p_k), shr32(p_x, 32 - p_k)) END;
    END rotl32;

    -----------------------------------------------------------------------------------------------------------------
    -----------------------------------------------------------------------------------------------------------------

    FUNCTION next64
        RETURN uint32_t
    IS
        --  uint32_t next(void) {
        --   const uint32_t s0 = s[0];
        --   uint32_t s1 = s[1];

        --   const uint32_t result = s0 * 0x9E3779BB;                  -- xoroshiro64*
        --   const uint32_t result = rotl(s0 * 0x9E3779BB, 5) * 5;     -- xoroshiro64**
        --
        --   s1 ^= s0;
        --   s[0] = rotl(s0, 26) ^ s1 ^ (s1 << 9); // a, b
        --   s[1] = rotl(s1, 13); // c
        --
        --   return result;
        --  }

        v_s0       uint32_t := g_64_state(0);
        v_s1       uint32_t := g_64_state(1);
        v_result   uint32_t
            := BITAND(
                   CASE g_mode
                       WHEN c_64star
                       THEN
                           v_s0 * TO_NUMBER('9E3779BB', 'xxxxxxxx')
                       WHEN c_64starstar
                       THEN
                           rotl32(BITAND(v_s0 * TO_NUMBER('9E3779BB', 'xxxxxxxx'), c_32bitmask), 5) * 5
                   END,
                   c_32bitmask);
    BEGIN
        v_s1 := bitxor32(v_s1, v_s0);
        g_64_state(0) := bitxor32(bitxor32(rotl32(v_s0, 26), v_s1), shl32(v_s1, 9));
        g_64_state(1) := rotl32(v_s1, 13);
        RETURN v_result;
    END next64;

    FUNCTION next128
        RETURN uint32_t
    IS
        --  uint32_t next(void) {
        --   const uint32_t result = s[0] + s[3];             -- xoshiro128+
        --   const uint32_t result = s[0] + s[3];             -- xoshiro128++
        --   const uint32_t result = rotl(s[1] * 5, 7) * 9;   -- xoshiro128**
        --
        --   const uint32_t t = s[1] << 9;
        --
        --   s[2] ^= s[0];
        --   s[3] ^= s[1];
        --   s[1] ^= s[2];
        --   s[0] ^= s[3];
        --
        --   s[2] ^= t;
        --
        --   s[3] = rotl(s[3], 11);
        --
        --   return result;
        --  }

        v_result   uint32_t
            := BITAND(
                   CASE g_mode
                       WHEN c_128plus
                       THEN
                           g_128_state(0) + g_128_state(3)
                       WHEN c_128plusplus
                       THEN
                           rotl32(BITAND(g_128_state(0) + g_128_state(3), c_32bitmask), 7) + g_128_state(0)
                       WHEN c_128starstar
                       THEN
                           rotl32(BITAND(g_128_state(1) * 5, c_32bitmask), 7) * 9
                   END,
                   c_32bitmask);

        v_t        uint32_t := shl32(g_128_state(1), 9);
    BEGIN
        g_128_state(2) := bitxor32(g_128_state(2), g_128_state(0));
        g_128_state(3) := bitxor32(g_128_state(3), g_128_state(1));
        g_128_state(1) := bitxor32(g_128_state(1), g_128_state(2));
        g_128_state(0) := bitxor32(g_128_state(0), g_128_state(3));

        g_128_state(2) := bitxor32(g_128_state(2), v_t);

        g_128_state(3) := rotl32(g_128_state(3), 11);
        RETURN v_result;
    END next128;

    PROCEDURE jump
    IS
        --  /* This is the jump function for the generator. It is equivalent
        --     to 2^64 calls to next(); it can be used to generate 2^64
        --     non-overlapping subsequences for parallel computations. */
        --
        --  void jump(void) {
        --   static const uint32_t JUMP[] = { 0x8764000b, 0xf542d2d3, 0x6fa035c3, 0x77f2db5b };
        --
        --   uint32_t s0 = 0;
        --   uint32_t s1 = 0;
        --   uint32_t s2 = 0;
        --   uint32_t s3 = 0;
        --   for(int i = 0; i < sizeof JUMP / sizeof *JUMP; i++)
        --    for(int b = 0; b < 32; b++) {
        --     if (JUMP[i] & UINT32_C(1) << b) {
        --      s0 ^= s[0];
        --      s1 ^= s[1];
        --      s2 ^= s[2];
        --      s3 ^= s[3];
        --     }
        --     next();
        --    }
        --
        --   s[0] = s0;
        --   s[1] = s1;
        --   s[2] = s2;
        --   s[3] = s3;
        --  }
        v_jump    state128_array_t
            := state128_array_t(0   => TO_NUMBER('8764000b', 'xxxxxxxx'),
                                1   => TO_NUMBER('f542d2d3', 'xxxxxxxx'),
                                2   => TO_NUMBER('6fa035c3', 'xxxxxxxx'),
                                3   => TO_NUMBER('77f2db5b', 'xxxxxxxx'));

        v_s0      uint32_t := 0;
        v_s1      uint32_t := 0;
        v_s2      uint32_t := 0;
        v_s3      uint32_t := 0;
        v_dummy   uint32_t;
    BEGIN
        CASE
            WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
            THEN
                FOR i IN 0 .. 3
                LOOP
                    FOR b IN 0 .. 31
                    LOOP
                        IF BITAND(v_jump(i), shl32(1, b)) != 0
                        THEN
                            v_s0 := bitxor32(v_s0, g_128_state(0));
                            v_s1 := bitxor32(v_s1, g_128_state(1));
                            v_s2 := bitxor32(v_s2, g_128_state(2));
                            v_s3 := bitxor32(v_s3, g_128_state(3));
                        END IF;

                        v_dummy := next128;
                    END LOOP;
                END LOOP;

                g_128_state(0) := v_s0;
                g_128_state(1) := v_s1;
                g_128_state(2) := v_s2;
                g_128_state(3) := v_s3;
            ELSE
                -- 64-bit state modes don't support jump and long_jump operations
                RAISE invalid_mode;
        END CASE;
    END jump;

    PROCEDURE long_jump
    IS
        --  /* This is the long-jump function for the generator. It is equivalent to
        --     2^96 calls to next(); it can be used to generate 2^32 starting points,
        --     from each of which jump() will generate 2^32 non-overlapping
        --     subsequences for parallel distributed computations. */
        --
        --  void long_jump(void) {
        --   static const uint32_t LONG_JUMP[] = { 0xb523952e, 0x0b6f099f, 0xccf5a0ef, 0x1c580662 };
        --
        --   uint32_t s0 = 0;
        --   uint32_t s1 = 0;
        --   uint32_t s2 = 0;
        --   uint32_t s3 = 0;
        --   for(int i = 0; i < sizeof LONG_JUMP / sizeof *LONG_JUMP; i++)
        --    for(int b = 0; b < 32; b++) {
        --     if (LONG_JUMP[i] & UINT32_C(1) << b) {
        --      s0 ^= s[0];
        --      s1 ^= s[1];
        --      s2 ^= s[2];
        --      s3 ^= s[3];
        --     }
        --     next();
        --    }
        --
        --   s[0] = s0;
        --   s[1] = s1;
        --   s[2] = s2;
        --   s[3] = s3;
        --  }

        v_jump    state128_array_t
            := state128_array_t(0   => TO_NUMBER('b523952e', 'xxxxxxxx'),
                                1   => TO_NUMBER('0b6f099f', 'xxxxxxxx'),
                                2   => TO_NUMBER('ccf5a0ef', 'xxxxxxxx'),
                                3   => TO_NUMBER('1c580662', 'xxxxxxxx'));

        v_s0      uint32_t := 0;
        v_s1      uint32_t := 0;
        v_s2      uint32_t := 0;
        v_s3      uint32_t := 0;
        v_dummy   uint32_t;
    BEGIN
        CASE
            WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
            THEN
                FOR i IN 0 .. 3
                LOOP
                    FOR b IN 0 .. 31
                    LOOP
                        IF BITAND(v_jump(i), shl32(1, b)) != 0
                        THEN
                            v_s0 := bitxor32(v_s0, g_128_state(0));
                            v_s1 := bitxor32(v_s1, g_128_state(1));
                            v_s2 := bitxor32(v_s2, g_128_state(2));
                            v_s3 := bitxor32(v_s3, g_128_state(3));
                        END IF;

                        v_dummy := next128;
                    END LOOP;
                END LOOP;

                g_128_state(0) := v_s0;
                g_128_state(1) := v_s1;
                g_128_state(2) := v_s2;
                g_128_state(3) := v_s3;
            ELSE
                -- 64-bit state modes don't support jump and long_jump operations
                RAISE invalid_mode;
        END CASE;
    END long_jump;

    -----------------------------------------------------------------------------------------------------------------
    -----------------------------------------------------------------------------------------------------------------
    PROCEDURE set_mode(p_mode IN mode_t)
    IS
    BEGIN
        g_mode := p_mode;
    END set_mode;

    PROCEDURE set_mode(p_mode IN VARCHAR2)
    IS
    BEGIN
        set_mode(
            CASE LOWER(p_mode)
                WHEN 'xoroshiro64*' THEN c_64star
                WHEN 'xoroshiro64**' THEN c_64starstar
                WHEN 'xoshiro128+' THEN c_128plus
                WHEN 'xoshiro128++' THEN c_128plusplus
                WHEN 'xoshiro128**' THEN c_128starstar
            END);
    END set_mode;

    FUNCTION get_mode
        RETURN mode_t
    IS
    BEGIN
        RETURN g_mode;
    END get_mode;

    FUNCTION get_mode_name
        RETURN VARCHAR2
    IS
    BEGIN
        RETURN CASE g_mode
                   WHEN c_64star THEN 'xoroshiro64*'
                   WHEN c_64starstar THEN 'xoroshiro64**'
                   WHEN c_128plus THEN 'xoshiro128+'
                   WHEN c_128plusplus THEN 'xoshiro128++'
                   WHEN c_128starstar THEN 'xoshiro128**'
               END;
    END get_mode_name;

    FUNCTION get_next
        RETURN uint32_t
    IS
    BEGIN
        RETURN CASE
                   WHEN g_mode IN (c_64star, c_64starstar) THEN next64
                   WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar) THEN next128
               END;
    END get_next;

    PROCEDURE set_state(p_seed IN seed_tab)
    IS
    BEGIN
        CASE
            WHEN g_mode IN (c_64star, c_64starstar)
            THEN
                g_64_state := state64_array_t();

                FOR i IN 0 .. 1
                LOOP
                    g_64_state(i) := p_seed(i + 1);
                END LOOP;
            WHEN g_mode IN (c_128plus, c_128plusplus, c_128starstar)
            THEN
                g_128_state := state128_array_t();

                FOR i IN 0 .. 3
                LOOP
                    g_128_state(i) := p_seed(i + 1);
                END LOOP;
            ELSE
                RAISE invalid_mode;
        END CASE;
    END set_state;

    -----------------------------------------------------------------------------------------------------------------
    -----------------------------------------------------------------------------------------------------------------

    PROCEDURE splitmix32_set_state(p_seed uint32_t)
    IS
    BEGIN
        g_splitmix32_state := p_seed;
    END splitmix32_set_state;

    FUNCTION splitmix32_next
        RETURN uint32_t
    IS
        -- splitmix32 is implemented using MurmurHash3
        -- with the extra step of adding the "Golden Gamma" 
        -- This modification is base on the SplitMix algorithm described by
        -- Guy L. Steele Jr., Doug Lea, and Christine H. Flood
        -- https://gee.cs.oswego.edu/dl/papers/oopsla14.pdf 
        --
        -- The MurmurHash3 algorithm by Austin Appleby can be found here.
        -- https://github.com/aappleby/smhasher/blob/master/src/MurmurHash3.cpp
        --

        -- The Golden Gamma is the closest odd number to 2^32/golden_ratio
        -- Where the golden ratio is (1+sqrt(5))/2
        c_golden_gamma   CONSTANT INTEGER := 2654435769; -- 0x9e3779b9

        v_z                       NUMBER;
    BEGIN
        --  In Appleby's code, the state is passed in as a parameter.
        --  In this package, the state is maintained as a package variable.
        g_splitmix32_state := BITAND(g_splitmix32_state + c_golden_gamma, c_32bitmask);
        v_z := g_splitmix32_state;

        v_z := bitxor32(v_z, shr32(v_z, 16));
        v_z := v_z * TO_NUMBER('85ebca6b', 'xxxxxxxx');
        v_z := BITAND(v_z, c_32bitmask);

        v_z := bitxor32(v_z, shr32(v_z, 13));
        v_z := v_z * TO_NUMBER('c2b2ae35', 'xxxxxxxx');
        v_z := BITAND(v_z, c_32bitmask);

        RETURN BITAND(bitxor32(v_z, shr32(v_z, 16)), c_32bitmask);
    END splitmix32_next;
END;
/

Because the seed numbers and intermediate results are smaller, these algorithm variants can use the native BITAND function and don’t requite the BIG_BITAND function for extremely large values.

This concludes the series and I hope you find the code interesting and useful. As mentioned in my 64-bit article, I intentionally split the 32-bit and 64-bit algorithms into distinct packages on the assumption they would be unlikely to be used together. But, if you find you do want them to be combined, I defined the mode constants such that they are mutually exclusive thus could be combined without fear of ambiguity from overlapping values.

Questions and comments, as always, are welcome.

Leave a Reply