CREATE OR REPLACE PACKAGE sds.sfmt IS -- PL/SQL Implementation of the SIMD oriented Fast Mersenne Twister(SFMT) -- Conversion to PL/SQL by Sean D. Stuber -- .///. -- (0 o) ---------------0000--(_)--0000--------------- -- -- Sean D. Stuber -- sean.stuber@gmail.com -- -- oooO Oooo --------------( )-----( )--------------- -- \ ( ) / -- \_) (_/ --Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima --University. All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- * Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- * Redistributions in binary form must reproduce the above -- copyright notice, this list of conditions and the following -- disclaimer in the documentation and/or other materials provided -- with the distribution. -- * Neither the name of the Hiroshima University nor the names of -- its contributors may be used to endorse or promote products -- derived from this software without specific prior written -- permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS --"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT --LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR --A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT --OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, --SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT --LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, --DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY --THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT --(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE --OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- -- PL/SQL doesn't have a 32-bit or 64-bit unsigned integer type -- Nor is there a way to declare a subtype with that constraint -- So we just have to fake it with regular INTEGER type -- Within the package, the code ensures no values execeed the bit limit. -- Inputs will be validated procedurally before accepting. SUBTYPE t_uint32 IS INTEGER; SUBTYPE t_uint64 IS INTEGER; TYPE t_uint32_array IS TABLE OF t_uint32 INDEX BY BINARY_INTEGER; TYPE t_uint64_array IS TABLE OF t_uint64 INDEX BY BINARY_INTEGER; -- The random number generator constructs an internal array of 32 bit values -- 64bit values are returned by combining 2 of these values -- If the offset into the array is an odd number this exception will be raised -- This exception should only occur if both 32bit and 64bit generation is -- used for a single initialization. To avoid this error, avoid using both -- types of generation, or be sure to invoke 32bit generation in pairs -- between 64bit calls. exc_invalid_64bit_offset EXCEPTION; -- Before the number generation can begin it must be initialized -- with a seed number or array of seed values -- Failure to do so will raise this exception exc_not_initialized EXCEPTION; -- If initialization value isn't a valid unsigned 32bit integer -- then this exception will be raised exc_invalid_uint32 EXCEPTION; -- If an invalid Mersenne Exponent is specified then -- this exception will be raised -- valid MEXP values are: (607,1279,2281,4253,11213,19937,44497,86243,132049,216091) exc_invalid_mexp EXCEPTION; -- Determines the period of the pseudo-random number generator -- Setting this defines several other parameters -- The original c-code defined these at compile time -- p_mexp must be one of the following -- 607, 1279, 2281, 4253, 11213, 19937, 44497, 86243, 132049, 216091 -- Changing the MEXP requires a reinitialization of the generator -- The default MEXP is initialized to 19937 during package initialization PROCEDURE set_mersenne_exponent(p_mexp IN INTEGER DEFAULT 19937); -- Simple initialization -- The MEXP (Mersenne Exponent) must be set prior to calling this PROCEDURE init_gen_rand(p_seed IN t_uint32); -- The initialization array must have 0..p_key_length-1 elements populated -- if key length is null, or unspecified the length will be p_init_key.count -- So, length isn't strictly necessary but the original api is maintained -- The MEXP (Mersenne Exponent) must be set prior to calling this PROCEDURE init_by_array( p_init_key IN OUT NOCOPY t_uint32_array, p_key_length INTEGER DEFAULT NULL ); -- Return a single unsigned 32-bit integer pseudo-random number -- subsequent calls will return new values -- The generator must be initialized prior to calling this FUNCTION gen_rand32 RETURN t_uint32; -- Return a single unsigned 64-bit integer pseudo-random number -- subsequent calls will return new values -- The generator must be initialized prior to calling this FUNCTION gen_rand64 RETURN t_uint64; -- Fills a 0-based array with pseudo-random unsigned 32-bit integers -- The generator must be initialized prior to calling this PROCEDURE fill_array32( p_array IN OUT NOCOPY t_uint32_array, p_size IN INTEGER ); -- Fills a 0-based array with pseudo-random unsigned 64-bit integers -- The generator must be initialized prior to calling this PROCEDURE fill_array64( p_array IN OUT NOCOPY t_uint64_array, p_size IN INTEGER ); /* These real versions are due to Isaku Wada */ -- generates a random number on [0,1]-real-interval FUNCTION genrand_real1 RETURN NUMBER; -- generates a random number on [0,1)-real-interval FUNCTION genrand_real2 RETURN NUMBER; -- generates a random number on (0,1)-real-interval FUNCTION genrand_real3 RETURN NUMBER; -- generates a random number on [0,1) with 53-bit resolution FUNCTION genrand_res53 RETURN NUMBER; -- generates a random number on [0,1) with 53-bit resolution -- value is generated from from two 32 bit integers FUNCTION genrand_res53_mix RETURN NUMBER; -- Replicate the checking processes of the original for whatever -- MEXP you want. Normally the results will be printed to -- the DBMS_OUTPUT buffer. -- If the directory and file name are provided the output will -- also be written via UTL_FILE to the specified file and location. PROCEDURE check32( p_mexp IN INTEGER DEFAULT 19937, p_file_dir IN VARCHAR2 DEFAULT NULL, p_file_name IN VARCHAR2 DEFAULT NULL ); PROCEDURE check64( p_mexp IN INTEGER DEFAULT 19937, p_file_dir IN VARCHAR2 DEFAULT NULL, p_file_name IN VARCHAR2 DEFAULT NULL ); -- Run speed tests PROCEDURE speed32(p_count IN INTEGER, p_block IN INTEGER); PROCEDURE speed64(p_count IN INTEGER, p_block IN INTEGER); END sfmt; / CREATE OR REPLACE PACKAGE BODY sds.sfmt IS -- PL/SQL Implementation of the SIMD oriented Fast Mersenne Twister(SFMT) -- Conversion to PL/SQL by Sean D. Stuber -- .///. -- (0 o) ---------------0000--(_)--0000--------------- -- -- Sean D. Stuber -- sean.stuber@gmail.com -- -- oooO Oooo --------------( )-----( )--------------- -- \ ( ) / -- \_) (_/ --Copyright (c) 2006,2007 Mutsuo Saito, Makoto Matsumoto and Hiroshima --University. All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- * Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- * Redistributions in binary form must reproduce the above -- copyright notice, this list of conditions and the following -- disclaimer in the documentation and/or other materials provided -- with the distribution. -- * Neither the name of the Hiroshima University nor the names of -- its contributors may be used to endorse or promote products -- derived from this software without specific prior written -- permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS --"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT --LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR --A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT --OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, --SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT --LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, --DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY --THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT --(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE --OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. c_8hexdigits CONSTANT VARCHAR2(8) := RPAD('x', 8, 'x'); c_16hexdigits CONSTANT VARCHAR2(16) := RPAD('x', 16, 'x'); c_32bit_mask CONSTANT INTEGER := TO_NUMBER(RPAD('f', 8, 'f'), c_8hexdigits) ; c_64bit_mask CONSTANT INTEGER := TO_NUMBER(RPAD('f', 16, 'f'), c_16hexdigits) ; ---- Pseudo-constants ---- These define the period of the generator ---- They are initialized once before processing ---- and remain constant thereafter ---- In the original c-code they were defined at compile time g_mexp INTEGER; g_pos1 INTEGER; g_sl1 INTEGER; g_sl2 INTEGER; g_sr1 INTEGER; g_sr2 INTEGER; g_msk1 INTEGER; g_msk2 INTEGER; g_msk3 INTEGER; g_msk4 INTEGER; g_parity1 INTEGER; g_parity2 INTEGER; g_parity3 INTEGER; g_parity4 INTEGER; g_idstr VARCHAR2(100); g_n INTEGER; g_n32 INTEGER; g_n64 INTEGER; g_period_set BOOLEAN := FALSE; ----------------------------------------------------------------------------------- TYPE t_w128 IS RECORD(u t_uint32_array); -- Ideally this would be a varray, but we can't declare it -- until the MEXP is defined so we'll use the more flexible -- associative array TYPE t_w128_array IS TABLE OF t_w128 INDEX BY PLS_INTEGER; g_sfmt t_w128_array; g_idx INTEGER; g_initialized BOOLEAN := FALSE; g_parity t_uint32_array; -- PL/SQL doesn't have syntax corresponding to the pointer -- tricks that c does with multi-typed pointers -- So, SET/GET routines to the nested array simulate them PROCEDURE psfmt32set(p_index IN INTEGER, p_value IN INTEGER) IS BEGIN g_sfmt(FLOOR(p_index / 4)).u(MOD(p_index, 4)) := p_value; END psfmt32set; FUNCTION psfmt32get(p_index IN INTEGER) RETURN INTEGER IS BEGIN RETURN g_sfmt(FLOOR(p_index / 4)).u(MOD(p_index, 4)); END psfmt32get; -- combines the functions of the GET and SET -- used where the value within the array is adjusted by some other -- value with 32bit roll over -- rather than simply assigning a new 32 bit value -- Statements like this: psfmt32[i] += r; PROCEDURE psfmt32delta(p_index IN INTEGER, p_delta IN INTEGER) IS x INTEGER; BEGIN x := g_sfmt(FLOOR(p_index / 4)).u(MOD(p_index, 4)) + p_delta; IF x < POWER(2, 32) THEN g_sfmt(FLOOR(p_index / 4)).u(MOD(p_index, 4)) := x; ELSE -- Rollover values greater than 32bits g_sfmt(FLOOR(p_index / 4)).u(MOD(p_index, 4)) := x - POWER(2, 32); END IF; END psfmt32delta; -- Determines the period of the pseudo-random number generator -- Setting this defines several other parameters -- The original c-code defined these at compile time -- p_mexp must be one of the following -- 607, 1279, 2281, 4253, 11213, 19937, 44497, 86243, 132049, 216091 PROCEDURE set_mersenne_exponent(p_mexp IN INTEGER DEFAULT 19937) IS BEGIN IF p_mexp IN (607, 1279, 2281, 4253, 11213, 19937, 44497, 86243, 132049, 216091) THEN -- If we change the exponent -- then the generator must be reinitalized g_initialized := FALSE; g_mexp := p_mexp; g_n := FLOOR(g_mexp / 128) + 1; g_n32 := g_n * 4; g_n64 := g_n * 2; IF p_mexp = 607 THEN g_pos1 := 2; g_sl1 := 15; g_sl2 := 3; g_sr1 := 13; g_sr2 := 3; g_msk1 := TO_NUMBER('fdff37ff', c_8hexdigits); g_msk2 := TO_NUMBER('ef7f3f7d', c_8hexdigits); g_msk3 := TO_NUMBER('ff777b7d', c_8hexdigits); g_msk4 := TO_NUMBER('7ff7fb2f', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('00000000', c_8hexdigits); g_parity4 := TO_NUMBER('5986f054', c_8hexdigits); g_idstr := 'SFMT-607:2-15-3-13-3:fdff37ff-ef7f3f7d-ff777b7d-7ff7fb2f'; ELSIF p_mexp = 1279 THEN g_pos1 := 7; g_sl1 := 14; g_sl2 := 3; g_sr1 := 5; g_sr2 := 1; g_msk1 := TO_NUMBER('f7fefffd', c_8hexdigits); g_msk2 := TO_NUMBER('7fefcfff', c_8hexdigits); g_msk3 := TO_NUMBER('aff3ef3f', c_8hexdigits); g_msk4 := TO_NUMBER('b5ffff7f', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('00000000', c_8hexdigits); g_parity4 := TO_NUMBER('20000000', c_8hexdigits); g_idstr := 'SFMT-1279:7-14-3-5-1:f7fefffd-7fefcfff-aff3ef3f-b5ffff7f'; ELSIF p_mexp = 2281 THEN g_pos1 := 12; g_sl1 := 19; g_sl2 := 1; g_sr1 := 5; g_sr2 := 1; g_msk1 := TO_NUMBER('bff7ffbf', c_8hexdigits); g_msk2 := TO_NUMBER('fdfffffe', c_8hexdigits); g_msk3 := TO_NUMBER('f7ffef7f', c_8hexdigits); g_msk4 := TO_NUMBER('f2f7cbbf', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('00000000', c_8hexdigits); g_parity4 := TO_NUMBER('41dfa600', c_8hexdigits); g_idstr := 'SFMT-2281:12-19-1-5-1:bff7ffbf-fdfffffe-f7ffef7f-f2f7cbbf'; ELSIF p_mexp = 4253 THEN g_pos1 := 17; g_sl1 := 20; g_sl2 := 1; g_sr1 := 7; g_sr2 := 1; g_msk1 := TO_NUMBER('9f7bffff', c_8hexdigits); g_msk2 := TO_NUMBER('9fffff5f', c_8hexdigits); g_msk3 := TO_NUMBER('3efffffb', c_8hexdigits); g_msk4 := TO_NUMBER('fffff7bb', c_8hexdigits); g_parity1 := TO_NUMBER('a8000001', c_8hexdigits); g_parity2 := TO_NUMBER('af5390a3', c_8hexdigits); g_parity3 := TO_NUMBER('b740b3f8', c_8hexdigits); g_parity4 := TO_NUMBER('6c11486d', c_8hexdigits); g_idstr := 'SFMT-4253:17-20-1-7-1:9f7bffff-9fffff5f-3efffffb-fffff7bb'; ELSIF p_mexp = 11213 THEN g_pos1 := 68; g_sl1 := 14; g_sl2 := 3; g_sr1 := 7; g_sr2 := 3; g_msk1 := TO_NUMBER('effff7fb', c_8hexdigits); g_msk2 := TO_NUMBER('ffffffef', c_8hexdigits); g_msk3 := TO_NUMBER('dfdfbfff', c_8hexdigits); g_msk4 := TO_NUMBER('7fffdbfd', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('e8148000', c_8hexdigits); g_parity4 := TO_NUMBER('d0c7afa3', c_8hexdigits); g_idstr := 'SFMT-11213:68-14-3-7-3:effff7fb-ffffffef-dfdfbfff-7fffdbfd'; ELSIF p_mexp = 19937 THEN g_pos1 := 122; g_sl1 := 18; g_sl2 := 1; g_sr1 := 11; g_sr2 := 1; g_msk1 := TO_NUMBER('dfffffef', c_8hexdigits); g_msk2 := TO_NUMBER('ddfecb7f', c_8hexdigits); g_msk3 := TO_NUMBER('bffaffff', c_8hexdigits); g_msk4 := TO_NUMBER('bffffff6', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('00000000', c_8hexdigits); g_parity4 := TO_NUMBER('13c9e684', c_8hexdigits); g_idstr := 'SFMT-19937:122-18-1-11-1:dfffffef-ddfecb7f-bffaffff-bffffff6'; ELSIF p_mexp = 44497 THEN g_pos1 := 330; g_sl1 := 5; g_sl2 := 3; g_sr1 := 9; g_sr2 := 3; g_msk1 := TO_NUMBER('effffffb', c_8hexdigits); g_msk2 := TO_NUMBER('dfbebfff', c_8hexdigits); g_msk3 := TO_NUMBER('bfbf7bef', c_8hexdigits); g_msk4 := TO_NUMBER('9ffd7bff', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('a3ac4000', c_8hexdigits); g_parity4 := TO_NUMBER('ecc1327a', c_8hexdigits); g_idstr := 'SFMT-44497:330-5-3-9-3:effffffb-dfbebfff-bfbf7bef-9ffd7bff'; ELSIF p_mexp = 86243 THEN g_pos1 := 366; g_sl1 := 6; g_sl2 := 7; g_sr1 := 19; g_sr2 := 1; g_msk1 := TO_NUMBER('fdbffbff', c_8hexdigits); g_msk2 := TO_NUMBER('bff7ff3f', c_8hexdigits); g_msk3 := TO_NUMBER('fd77efff', c_8hexdigits); g_msk4 := TO_NUMBER('bf9ff3ff', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('00000000', c_8hexdigits); g_parity4 := TO_NUMBER('e9528d85', c_8hexdigits); g_idstr := 'SFMT-86243:366-6-7-19-1:fdbffbff-bff7ff3f-fd77efff-bf9ff3ff'; ELSIF p_mexp = 132049 THEN g_pos1 := 110; g_sl1 := 19; g_sl2 := 1; g_sr1 := 21; g_sr2 := 1; g_msk1 := TO_NUMBER('ffffbb5f', c_8hexdigits); g_msk2 := TO_NUMBER('fb6ebf95', c_8hexdigits); g_msk3 := TO_NUMBER('fffefffa', c_8hexdigits); g_msk4 := TO_NUMBER('cff77fff', c_8hexdigits); g_parity1 := TO_NUMBER('00000001', c_8hexdigits); g_parity2 := TO_NUMBER('00000000', c_8hexdigits); g_parity3 := TO_NUMBER('cb520000', c_8hexdigits); g_parity4 := TO_NUMBER('c7e91c7d', c_8hexdigits); g_idstr := 'SFMT-132049:110-19-1-21-1:ffffbb5f-fb6ebf95-fffefffa-cff77fff'; ELSIF p_mexp = 216091 THEN g_pos1 := 627; g_sl1 := 11; g_sl2 := 3; g_sr1 := 10; g_sr2 := 1; g_msk1 := TO_NUMBER('bff7bff7', c_8hexdigits); g_msk2 := TO_NUMBER('bfffffff', c_8hexdigits); g_msk3 := TO_NUMBER('bffffa7f', c_8hexdigits); g_msk4 := TO_NUMBER('ffddfbfb', c_8hexdigits); g_parity1 := TO_NUMBER('f8000001', c_8hexdigits); g_parity2 := TO_NUMBER('89e80709', c_8hexdigits); g_parity3 := TO_NUMBER('3bd2b64b', c_8hexdigits); g_parity4 := TO_NUMBER('0c64b1e4', c_8hexdigits); g_idstr := 'SFMT-216091:627-11-3-10-1:bff7bff7-bfffffff-bffffa7f-ffddfbfb'; END IF; g_parity(0) := g_parity1; g_parity(1) := g_parity2; g_parity(2) := g_parity3; g_parity(3) := g_parity4; ELSE RAISE exc_invalid_mexp; END IF; END set_mersenne_exponent; -- Return an integer that is a result of OR-ing each bit of a and b FUNCTION bitor(p_a IN INTEGER, p_b IN INTEGER) RETURN INTEGER IS BEGIN RETURN (p_a + p_b) - BITAND(p_a, p_b); END bitor; -- Return an integer that is a result of XOR-ing each bit of a and b -- enforce 32bit result FUNCTION bitxor(p_a IN INTEGER, p_b IN INTEGER) RETURN INTEGER IS BEGIN RETURN BITAND((p_a + p_b) - BITAND(p_a, p_b) * 2, c_32bit_mask); END bitxor; -- Shift bits right FUNCTION bitshr(p_bits IN INTEGER, p_shift IN INTEGER) RETURN INTEGER IS BEGIN RETURN FLOOR(p_bits / (2 ** p_shift)); END bitshr; -- Shift bits left, and enforce 32bit or 64bit results FUNCTION bitshl( p_bits IN INTEGER, p_shift IN INTEGER, p_length IN INTEGER DEFAULT 32 ) RETURN INTEGER IS BEGIN IF p_length = 32 THEN RETURN BITAND(p_bits * (2 ** p_shift), c_32bit_mask); ELSE RETURN BITAND(p_bits * (2 ** p_shift), c_64bit_mask); END IF; END bitshl; PROCEDURE rshift128( p_out IN OUT NOCOPY t_w128, p_in IN t_w128, p_shift IN INTEGER ) IS th INTEGER; tl INTEGER; oh INTEGER; ol INTEGER; BEGIN th := bitor(bitshl(p_in.u(3), 32, 64), p_in.u(2)); tl := bitor(bitshl(p_in.u(1), 32, 64), p_in.u(0)); oh := bitshr(th, p_shift * 8); ol := bitshr(tl, p_shift * 8); ol := bitor(ol, bitshl(th, 64 - p_shift * 8, 64)); p_out.u(1) := BITAND(bitshr(ol, 32), c_32bit_mask); p_out.u(0) := BITAND(ol, c_32bit_mask); p_out.u(3) := BITAND(bitshr(oh, 32), c_32bit_mask); p_out.u(2) := BITAND(oh, c_32bit_mask); END rshift128; PROCEDURE lshift128( p_out IN OUT NOCOPY t_w128, p_in IN t_w128, p_shift IN INTEGER ) IS th INTEGER; tl INTEGER; oh INTEGER; ol INTEGER; BEGIN th := bitor(bitshl(p_in.u(3), 32, 64), p_in.u(2)); tl := bitor(bitshl(p_in.u(1), 32, 64), p_in.u(0)); oh := bitshl(th, p_shift * 8, 64); ol := bitshl(tl, p_shift * 8, 64); oh := bitor(oh, bitshr(tl, 64 - p_shift * 8)); p_out.u(1) := BITAND(bitshr(ol, 32), c_32bit_mask); p_out.u(0) := BITAND(ol, c_32bit_mask); p_out.u(3) := BITAND(bitshr(oh, 32), c_32bit_mask); p_out.u(2) := BITAND(oh, c_32bit_mask); END lshift128; PROCEDURE do_recursion( p_r IN OUT NOCOPY t_w128, p_a IN t_w128, p_b IN t_w128, p_c IN INTEGER, p_d IN INTEGER ) IS x t_w128; y t_w128; BEGIN lshift128(x, p_a, g_sl2); rshift128(y, g_sfmt(p_c), g_sr2); p_r.u(0) := bitxor( bitxor( bitxor( bitxor(p_a.u(0), x.u(0)), BITAND(bitshr(p_b.u(0), g_sr1), g_msk1) ), y.u(0) ), bitshl(g_sfmt(p_d).u(0), g_sl1) ); p_r.u(1) := bitxor( bitxor( bitxor( bitxor(p_a.u(1), x.u(1)), BITAND(bitshr(p_b.u(1), g_sr1), g_msk2) ), y.u(1) ), bitshl(g_sfmt(p_d).u(1), g_sl1) ); p_r.u(2) := bitxor( bitxor( bitxor( bitxor(p_a.u(2), x.u(2)), BITAND(bitshr(p_b.u(2), g_sr1), g_msk3) ), y.u(2) ), bitshl(g_sfmt(p_d).u(2), g_sl1) ); p_r.u(3) := bitxor( bitxor( bitxor( bitxor(p_a.u(3), x.u(3)), BITAND(bitshr(p_b.u(3), g_sr1), g_msk4) ), y.u(3) ), bitshl(g_sfmt(p_d).u(3), g_sl1) ); END do_recursion; PROCEDURE gen_rand_all IS i INTEGER; r1 INTEGER; r2 INTEGER; BEGIN r1 := g_n - 2; r2 := g_n - 1; i := 0; WHILE i < g_n - g_pos1 LOOP do_recursion( g_sfmt(i), g_sfmt(i), g_sfmt(i + g_pos1), r1, r2 ); r1 := r2; r2 := i; i := i + 1; END LOOP; WHILE i < g_n LOOP do_recursion( g_sfmt(i), g_sfmt(i), g_sfmt(i + g_pos1 - g_n), r1, r2 ); r1 := r2; r2 := i; i := i + 1; END LOOP; END gen_rand_all; FUNCTION func1(p_num IN INTEGER) RETURN INTEGER IS x INTEGER := BITAND(p_num, c_32bit_mask); BEGIN RETURN BITAND(bitxor(x, bitshr(x, 27)) * 1664525, c_32bit_mask); END func1; FUNCTION func2(p_num IN INTEGER) RETURN INTEGER IS BEGIN RETURN BITAND( bitxor(p_num, bitshr(p_num, 27)) * 1566083941, c_32bit_mask ); END func2; PROCEDURE period_certification IS v_inner INTEGER := 0; v_work INTEGER; v_shift INTEGER; BEGIN FOR i IN 0 .. 3 LOOP v_inner := bitxor(v_inner, BITAND(psfmt32get(i), g_parity(i))); END LOOP; v_shift := 16; WHILE v_shift > 0 LOOP v_inner := bitxor(v_inner, bitshr(v_inner, v_shift)); v_shift := FLOOR(v_shift / 2); END LOOP; v_inner := BITAND(v_inner, 1); IF (v_inner = 1) THEN RETURN; END IF; FOR i IN 0 .. 3 LOOP v_work := 1; FOR j IN 0 .. 31 LOOP IF (BITAND(v_work, g_parity(i)) != 0) THEN psfmt32set(i, bitxor(psfmt32get(i), v_work)); RETURN; END IF; v_work := BITAND(v_work * 2, c_32bit_mask); END LOOP; END LOOP; END period_certification; FUNCTION get_idstring RETURN VARCHAR2 IS BEGIN RETURN g_idstr; END get_idstring; FUNCTION gen_rand32 RETURN t_uint32 IS r INTEGER; BEGIN IF NOT g_initialized THEN RAISE exc_not_initialized; END IF; IF (g_idx >= g_n32) THEN gen_rand_all; g_idx := 0; END IF; r := psfmt32get(g_idx); g_idx := g_idx + 1; RETURN r; END gen_rand32; FUNCTION gen_rand64 RETURN t_uint64 IS v_lo INTEGER; v_hi INTEGER; BEGIN IF NOT g_initialized THEN RAISE exc_not_initialized; END IF; IF MOD(g_idx, 2) != 0 THEN RAISE exc_invalid_64bit_offset; END IF; IF (g_idx >= g_n32) THEN gen_rand_all; g_idx := 0; END IF; v_lo := psfmt32get(g_idx); v_hi := psfmt32get(g_idx + 1); g_idx := g_idx + 2; RETURN bitor(bitshl(v_hi, 32, 64), v_lo); END gen_rand64; -- Array filling in the original c-code relied on pointer manipulation for -- rapid array population. No corresponding facility exists for pl/sql -- so the array is populated by simple iteration. -- As a side effect of this, some of the extra assertions on the -- array limitations aren't necessary in this implementation PROCEDURE fill_array32( p_array IN OUT NOCOPY t_uint32_array, p_size IN INTEGER ) IS v_temp_w128 t_w128_array; BEGIN IF NOT g_initialized THEN RAISE exc_not_initialized; END IF; p_array.delete; FOR i IN 0 .. p_size - 1 LOOP p_array(i) := gen_rand32; END LOOP; END fill_array32; PROCEDURE fill_array64( p_array IN OUT NOCOPY t_uint64_array, p_size IN INTEGER ) IS v_temp_w128 t_w128_array; v_x INTEGER; v_y INTEGER; BEGIN IF NOT g_initialized THEN RAISE exc_not_initialized; END IF; p_array.delete; FOR i IN 0 .. p_size - 1 LOOP v_x := gen_rand32; v_y := gen_rand32; p_array(i) := bitor(bitshl(v_y, 32, 64), v_x); -- END LOOP; END fill_array64; PROCEDURE init_gen_rand(p_seed IN t_uint32) IS BEGIN -- Validate the seed value is within the -- unsigned 32bit integer range 0..2^32-1 IF p_seed < 0 OR p_seed >= POWER(2, 32) OR p_seed IS NULL THEN RAISE exc_invalid_uint32; END IF; g_sfmt.delete; psfmt32set(0, p_seed); FOR i IN 1 .. g_n32 - 1 LOOP psfmt32set( i, BITAND( 1812433253 * bitxor( psfmt32get(i - 1), bitshr(psfmt32get(i - 1), 30) ) + i, c_32bit_mask ) ); END LOOP; g_idx := g_n32; period_certification; g_initialized := TRUE; END init_gen_rand; PROCEDURE init_by_array( p_init_key IN OUT NOCOPY t_uint32_array, p_key_length INTEGER DEFAULT NULL ) IS i INTEGER; j INTEGER; v_count INTEGER; v_r INTEGER; v_lag INTEGER; v_mid INTEGER; v_size INTEGER := g_n * 4; v_length INTEGER := NVL(p_key_length, p_init_key.COUNT); BEGIN -- Validate the initialization array is fully populated -- on the required range. Also confirm each element -- is within the 0-2^32-1 unsigned 32-bit integer range. FOR i IN 0 .. v_length - 1 LOOP IF p_init_key(i) < 0 OR p_init_key(i) >= POWER(2, 32) OR p_init_key(i) IS NULL THEN RAISE exc_invalid_uint32; END IF; END LOOP; g_sfmt.delete; IF v_size >= 623 THEN v_lag := 11; ELSIF v_size >= 68 THEN v_lag := 7; ELSIF v_size >= 39 THEN v_lag := 5; ELSE v_lag := 3; END IF; v_mid := FLOOR((v_size - v_lag) / 2); -- quasi-memset FOR p IN 0 .. g_n - 1 LOOP FOR q IN 0 .. 4 LOOP g_sfmt(p).u(q) := TO_NUMBER('8b8b8b8b', c_8hexdigits); END LOOP; END LOOP; IF v_length + 1 > g_n32 THEN v_count := v_length + 1; ELSE v_count := g_n32; END IF; v_r := func1( bitxor( bitxor(psfmt32get(0), psfmt32get(v_mid)), psfmt32get(g_n32 - 1) ) ); psfmt32delta(v_mid, v_r); v_r := v_r + v_length; psfmt32delta(v_mid + v_lag, v_r); psfmt32set(0, v_r); v_count := v_count - 1; i := 1; j := 0; WHILE j < v_count AND j < v_length LOOP v_r := func1( bitxor( bitxor( psfmt32get(i), psfmt32get(MOD(i + v_mid, g_n32)) ), psfmt32get(MOD(i + g_n32 - 1, g_n32)) ) ); psfmt32delta(MOD(i + v_mid, g_n32), v_r); v_r := v_r + p_init_key(j) + i; psfmt32delta(MOD(i + v_mid + v_lag, g_n32), v_r); psfmt32set(i, v_r); i := MOD(i + 1, g_n32); j := j + 1; END LOOP; WHILE j < v_count LOOP v_r := func1( bitxor( bitxor( psfmt32get(i), psfmt32get(MOD(i + v_mid, g_n32)) ), psfmt32get(MOD(i + g_n32 - 1, g_n32)) ) ); psfmt32delta(MOD(i + v_mid, g_n32), v_r); v_r := v_r + i; psfmt32delta(MOD(i + v_mid + v_lag, g_n32), v_r); psfmt32set(i, v_r); i := MOD(i + 1, g_n32); j := j + 1; END LOOP; j := 0; WHILE j < g_n32 LOOP v_r := func2( BITAND( BITAND( psfmt32get(i) + psfmt32get(MOD(i + v_mid, g_n32)), c_32bit_mask ) + psfmt32get(MOD(i + g_n32 - 1, g_n32)), c_32bit_mask ) ); psfmt32set( MOD(i + v_mid, g_n32), bitxor(psfmt32get(MOD(i + v_mid, g_n32)), v_r) ); v_r := v_r - i; psfmt32set( MOD(i + v_mid + v_lag, g_n32), bitxor(psfmt32get(MOD(i + v_mid + v_lag, g_n32)), v_r) ); psfmt32set(i, v_r); i := MOD(i + 1, g_n32); j := j + 1; END LOOP; g_idx := g_n32; period_certification; g_initialized := TRUE; END init_by_array; -- generates a random number on [0,1]-real-interval FUNCTION genrand_real1 RETURN NUMBER IS BEGIN RETURN gen_rand32 / 4294967295; END genrand_real1; -- generates a random number on [0,1)-real-interval FUNCTION genrand_real2 RETURN NUMBER IS BEGIN RETURN gen_rand32 / 4294967296; END genrand_real2; -- generates a random number on (0,1)-real-interval FUNCTION genrand_real3 RETURN NUMBER IS BEGIN RETURN (gen_rand32 + 0.5) / 4294967296; END genrand_real3; FUNCTION genrand_res53 RETURN NUMBER IS BEGIN RETURN gen_rand64() / 18446744073709551616; END genrand_res53; FUNCTION genrand_res53_mix RETURN NUMBER IS BEGIN -- Functionally this is the same as genrand_res53 -- replicated here for to mimic origional c api RETURN bitor(gen_rand32, bitshl(gen_rand32, 32, 64)) / 18446744073709551616; END genrand_res53_mix; PROCEDURE check32( p_mexp IN INTEGER DEFAULT 19937, p_file_dir IN VARCHAR2 DEFAULT NULL, p_file_name IN VARCHAR2 DEFAULT NULL ) IS v_out1 t_uint32_array; v_out2 t_uint32_array; v_init t_uint32_array; v_str VARCHAR2(1000); v_file UTL_FILE.file_type; v_r INTEGER; PROCEDURE p(v_str IN VARCHAR2, p_newline IN BOOLEAN DEFAULT FALSE) IS BEGIN IF UTL_FILE.is_open(v_file) THEN IF p_newline THEN UTL_FILE.put_line(v_file, v_str); ELSE UTL_FILE.put(v_file, v_str); END IF; END IF; IF p_newline THEN DBMS_OUTPUT.put_line(v_str); ELSE DBMS_OUTPUT.put(v_str); END IF; END; BEGIN set_mersenne_exponent(p_mexp); IF p_file_dir IS NOT NULL AND p_file_name IS NOT NULL THEN v_file := UTL_FILE.fopen( p_file_dir, p_file_name, 'w', 1000 ); END IF; p(get_idstring, TRUE); p('32 bit generated randoms', TRUE); p('init_gen_rand__________', TRUE); init_gen_rand(1234); fill_array32(v_out1, 10000); fill_array32(v_out2, 10000); init_gen_rand(1234); FOR i IN 0 .. 9999 LOOP IF i < 1000 THEN p(LPAD(v_out1(i), 10) || ' '); IF MOD(i, 5) = 4 THEN p(NULL, TRUE); END IF; END IF; v_r := gen_rand32; IF v_r != v_out1(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out1:' || TO_CHAR(v_out1(i), c_8hexdigits) || ' gen:' || TO_CHAR(v_r, c_8hexdigits), TRUE ); RETURN; END IF; END LOOP; FOR i IN 0 .. 699 LOOP v_r := gen_rand32; IF v_r != v_out2(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out2:' || TO_CHAR(v_out2(i), c_8hexdigits) || ' gen:' || TO_CHAR(v_r, c_8hexdigits), TRUE ); RETURN; END IF; END LOOP; p(CHR(10) || 'init_by_array__________', TRUE); v_init(0) := TO_NUMBER('1234', 'xxxx'); v_init(1) := TO_NUMBER('5678', 'xxxx'); v_init(2) := TO_NUMBER('9abc', 'xxxx'); v_init(3) := TO_NUMBER('def0', 'xxxx'); init_by_array(v_init, 4); fill_array32(v_out1, 10000); fill_array32(v_out2, 10000); init_by_array(v_init, 4); FOR i IN 0 .. 9999 LOOP IF i < 1000 THEN p(LPAD(v_out1(i), 10) || ' '); IF MOD(i, 5) = 4 THEN p(NULL, TRUE); END IF; END IF; v_r := gen_rand32; IF v_r != v_out1(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out1:' || TO_CHAR(v_out1(i), c_8hexdigits) || ' gen:' || TO_CHAR(v_r, c_8hexdigits), TRUE ); RETURN; END IF; END LOOP; FOR i IN 0 .. 699 LOOP v_r := gen_rand32; IF v_r != v_out2(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out2:' || TO_CHAR(v_out2(i), c_8hexdigits) || ' gen:' || TO_CHAR(v_r, c_8hexdigits), TRUE ); RETURN; END IF; END LOOP; IF UTL_FILE.is_open(v_file) THEN UTL_FILE.fclose(v_file); END IF; END check32; PROCEDURE check64( p_mexp IN INTEGER DEFAULT 19937, p_file_dir IN VARCHAR2 DEFAULT NULL, p_file_name IN VARCHAR2 DEFAULT NULL ) IS v_out1 t_uint64_array; v_out2 t_uint64_array; v_init t_uint32_array; v_str VARCHAR2(1000); v_file UTL_FILE.file_type; v_r INTEGER; PROCEDURE p(v_str IN VARCHAR2, p_newline IN BOOLEAN DEFAULT FALSE) IS BEGIN IF UTL_FILE.is_open(v_file) THEN IF p_newline THEN UTL_FILE.put_line(v_file, v_str); ELSE UTL_FILE.put(v_file, v_str); END IF; END IF; IF p_newline THEN DBMS_OUTPUT.put_line(v_str); ELSE DBMS_OUTPUT.put(v_str); END IF; END; BEGIN set_mersenne_exponent(p_mexp); IF p_file_dir IS NOT NULL AND p_file_name IS NOT NULL THEN v_file := UTL_FILE.fopen( p_file_dir, p_file_name, 'w', 1000 ); END IF; p(get_idstring, TRUE); p('64 bit generated randoms', TRUE); p('init_gen_rand__________', TRUE); init_gen_rand(4321); fill_array64(v_out1, 5000); fill_array64(v_out2, 5000); init_gen_rand(4321); FOR i IN 0 .. 4999 LOOP IF i < 1000 THEN p(LPAD(v_out1(i), 20) || ' '); IF MOD(i, 3) = 2 THEN p(NULL, TRUE); END IF; END IF; v_r := gen_rand64; IF v_r != v_out1(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out1:PRI' || TO_CHAR(v_out1(i), c_16hexdigits) || ' gen:PRI' || TO_CHAR(v_r, c_16hexdigits), TRUE ); RETURN; END IF; END LOOP; FOR i IN 0 .. 699 LOOP v_r := gen_rand64; IF v_r != v_out2(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out2:PRI' || TO_CHAR(v_out2(i), c_16hexdigits) || ' gen:PRI' || TO_CHAR(v_r, c_16hexdigits), TRUE ); RETURN; END IF; END LOOP; p(CHR(10) || 'init_by_array__________', TRUE); v_init(0) := 5; v_init(1) := 4; v_init(2) := 3; v_init(3) := 2; v_init(4) := 1; init_by_array(v_init, 5); fill_array64(v_out1, 5000); fill_array64(v_out2, 5000); init_by_array(v_init, 5); FOR i IN 0 .. 4999 LOOP IF i < 1000 THEN p(LPAD(v_out1(i), 20) || ' '); IF MOD(i, 3) = 2 THEN p(NULL, TRUE); END IF; END IF; v_r := gen_rand64; IF v_r != v_out1(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out1:' || TO_CHAR(v_out1(i), c_16hexdigits) || ' gen:' || TO_CHAR(v_r, c_16hexdigits), TRUE ); RETURN; END IF; END LOOP; FOR i IN 0 .. 699 LOOP v_r := gen_rand64; IF v_r != v_out2(i) THEN p( CHR(10) || 'mismatch at ' || i || ' v_out2:PRI' || TO_CHAR(v_out2(i), c_16hexdigits) || ' gen:PRI' || TO_CHAR(v_r, c_16hexdigits), TRUE ); RETURN; END IF; END LOOP; IF UTL_FILE.is_open(v_file) THEN UTL_FILE.fclose(v_file); END IF; END check64; PROCEDURE speed32(p_count IN INTEGER, p_block IN INTEGER) IS v_start TIMESTAMP; v_end TIMESTAMP; v_min INTERVAL DAY TO SECOND := NUMTODSINTERVAL(99, 'DAY'); v_array t_uint32_array; v_rand t_uint32; FUNCTION interval2ms(p_interval IN INTERVAL DAY TO SECOND) RETURN INTEGER IS BEGIN RETURN 1000 * ( EXTRACT(DAY FROM p_interval) * 86400 + EXTRACT(HOUR FROM p_interval) * 3600 + EXTRACT(MINUTE FROM p_interval) * 60 + EXTRACT(SECOND FROM p_interval)); END; BEGIN init_gen_rand(1234); FOR i IN 1 .. 10 LOOP v_start := SYSTIMESTAMP; -- Iterate 10 times and report the minimum timing FOR j IN 1 .. p_count LOOP fill_array32(v_array, p_block); END LOOP; v_end := SYSTIMESTAMP; IF v_end - v_start < v_min THEN v_min := v_end - v_start; END IF; END LOOP; DBMS_OUTPUT.put_line( '32 bit BLOCK:' || interval2ms(v_min) || 'ms for ' || p_block * p_count || ' randoms generation' ); init_gen_rand(1234); FOR i IN 1 .. 10 LOOP v_start := SYSTIMESTAMP; -- Iterate 10 times and report the minimum timing FOR j IN 1 .. p_count * p_block LOOP v_rand := gen_rand32; END LOOP; v_end := SYSTIMESTAMP; IF v_end - v_start < v_min THEN v_min := v_end - v_start; END IF; END LOOP; DBMS_OUTPUT.put_line( '32 bit SEQUE:' || interval2ms(v_min) || 'ms for ' || p_block * p_count || ' randoms generation' ); END speed32; PROCEDURE speed64(p_count IN INTEGER, p_block IN INTEGER) IS v_start TIMESTAMP; v_end TIMESTAMP; v_min INTERVAL DAY TO SECOND := NUMTODSINTERVAL(99, 'DAY'); v_array t_uint64_array; v_rand t_uint64; FUNCTION interval2ms(p_interval IN INTERVAL DAY TO SECOND) RETURN INTEGER IS BEGIN RETURN 1000 * ( EXTRACT(DAY FROM p_interval) * 86400 + EXTRACT(HOUR FROM p_interval) * 3600 + EXTRACT(MINUTE FROM p_interval) * 60 + EXTRACT(SECOND FROM p_interval)); END; BEGIN init_gen_rand(1234); FOR i IN 1 .. 10 LOOP v_start := SYSTIMESTAMP; -- Iterate 10 times and report the minimum timing FOR j IN 1 .. p_count LOOP fill_array64(v_array, p_block); END LOOP; v_end := SYSTIMESTAMP; IF v_end - v_start < v_min THEN v_min := v_end - v_start; END IF; END LOOP; DBMS_OUTPUT.put_line( '64 bit BLOCK:' || interval2ms(v_min) || 'ms for ' || p_block * p_count || ' randoms generation' ); init_gen_rand(1234); FOR i IN 1 .. 10 LOOP v_start := SYSTIMESTAMP; -- Iterate 10 times and report the minimum timing FOR j IN 1 .. p_count * p_block LOOP v_rand := gen_rand64; END LOOP; v_end := SYSTIMESTAMP; IF v_end - v_start < v_min THEN v_min := v_end - v_start; END IF; END LOOP; DBMS_OUTPUT.put_line( '64 bit SEQUE:' || interval2ms(v_min) || 'ms for ' || p_block * p_count || ' randoms generation' ); END speed64; BEGIN -- Initialize with default Mersenne Exponent parameters set_mersenne_exponent; END; /