250 lines
6.1 KiB
Plaintext
250 lines
6.1 KiB
Plaintext
note
|
|
description: "A Linear congruential random number generator"
|
|
author: "Colin LeMahieu"
|
|
date: "$Date$"
|
|
revision: "$Revision$"
|
|
quote: "To limit the press is to insult a nation; to prohibit reading of certain books is to declare the inhabitants to be either fools or knaves. - Claude-Adrien Helvetius"
|
|
|
|
class
|
|
LINEAR_CONGRUENTIAL_RNG
|
|
|
|
inherit
|
|
RANDOM_NUMBER_GENERATOR
|
|
LIMB_MANIPULATION
|
|
SPECIAL_UTILITY
|
|
SPECIAL_ARITHMETIC
|
|
rename
|
|
add as add_special,
|
|
sub as sub_special,
|
|
mul as mul_special
|
|
export
|
|
{NONE}
|
|
all
|
|
end
|
|
INTEGER_X_DIVISION
|
|
|
|
create
|
|
make
|
|
|
|
feature
|
|
|
|
make (size: INTEGER)
|
|
require
|
|
size >= 0
|
|
size <= 128
|
|
local
|
|
a: INTEGER_X
|
|
m2exp: INTEGER
|
|
number: STRING_8
|
|
do
|
|
inspect size
|
|
when 0..16 then
|
|
m2exp := 32
|
|
number := "29cf535"
|
|
when 17 then
|
|
m2exp := 34
|
|
number := "A3D73AD"
|
|
when 18 then
|
|
m2exp := 36
|
|
number := "28f825c5"
|
|
when 19 then
|
|
m2exp := 38
|
|
number := "a3dd4cdd"
|
|
when 20 then
|
|
m2exp := 40
|
|
number := "28f5da175"
|
|
when 21..28 then
|
|
m2exp := 56
|
|
number := "AA7D735234C0DD"
|
|
when 29..32 then
|
|
m2exp := 64
|
|
number := "BAECD515DAF0B49D"
|
|
when 33..50 then
|
|
m2exp := 100
|
|
number := "292787EBD3329AD7E7575E2FD"
|
|
when 51..64 then
|
|
m2exp := 128
|
|
number := "48A74F367FA7B5C8ACBB36901308FA85"
|
|
when 65..78 then
|
|
m2exp := 156
|
|
number := "78A7FDDDC43611B527C3F1D760F36E5D7FC7C45"
|
|
when 79..98 then
|
|
m2exp := 196
|
|
number := "41BA2E104EE34C66B3520CE706A56498DE6D44721E5E24F5"
|
|
when 99..100 then
|
|
m2exp := 200
|
|
number := "4E5A24C38B981EAFE84CD9D0BEC48E83911362C114F30072C5"
|
|
when 101..238 then
|
|
m2exp := 256
|
|
number := "AF66BA932AAF58A071FD8F0742A99A0C76982D648509973DB802303128A14CB5"
|
|
end
|
|
create a.make_from_string_base (number, 16)
|
|
randinit_lc_2exp (a, 1, m2exp)
|
|
end
|
|
|
|
feature
|
|
|
|
randinit_lc_2exp (a: READABLE_INTEGER_X; c: NATURAL_32; m2exp: INTEGER)
|
|
local
|
|
seedn: INTEGER
|
|
do
|
|
seedn := bits_to_limbs (m2exp)
|
|
create seed.make (m2exp, a, seedn, c)
|
|
end
|
|
|
|
seed: RAND_LC_STRUCT
|
|
|
|
randseed (seed_a: READABLE_INTEGER_X)
|
|
local
|
|
seedz: READABLE_INTEGER_X
|
|
seedn: INTEGER
|
|
do
|
|
seedz := seed.seed
|
|
seedn := bits_to_limbs (seed.m2exp)
|
|
fdiv_r_2exp (seedz, seed_a, seed.m2exp)
|
|
seedz.item.fill_with (0, seedz.count, seedz.count - (seedn - seedz.count))
|
|
seedz.count := seedn
|
|
end
|
|
|
|
randget (target: SPECIAL [NATURAL_32]; target_offset: INTEGER; nbits: INTEGER)
|
|
local
|
|
rbitpos: INTEGER
|
|
chunk_nbits: INTEGER
|
|
tp: SPECIAL [NATURAL_32]
|
|
tp_offset: INTEGER
|
|
tn: INTEGER
|
|
r2p: SPECIAL [NATURAL_32]
|
|
r2p_offset: INTEGER
|
|
rcy: NATURAL_32
|
|
junk: INTEGER
|
|
last_nbits: INTEGER
|
|
savelimb: NATURAL_32
|
|
carry: CELL [NATURAL_32]
|
|
do
|
|
create carry.put (0)
|
|
chunk_nbits := seed.m2exp // 2
|
|
tn := bits_to_limbs (chunk_nbits)
|
|
create tp.make_filled (0, tn)
|
|
from
|
|
rbitpos := 0
|
|
until
|
|
rbitpos + chunk_nbits > nbits
|
|
loop
|
|
r2p := target
|
|
r2p_offset := target_offset + rbitpos // limb_bits
|
|
if rbitpos \\ limb_bits /= 0 then
|
|
junk := lc (tp, tp_offset)
|
|
savelimb := r2p [r2p_offset]
|
|
lshift (r2p, r2p_offset, tp, tp_offset, tn, rbitpos \\ limb_bits, carry)
|
|
rcy := carry.item
|
|
r2p [r2p_offset] := r2p [r2p_offset].bit_or (savelimb)
|
|
if (chunk_nbits \\ limb_bits + rbitpos \\ limb_bits) > limb_bits then
|
|
r2p [r2p_offset + tn] := rcy
|
|
end
|
|
else
|
|
junk := lc (r2p, r2p_offset)
|
|
end
|
|
rbitpos := rbitpos + chunk_nbits
|
|
end
|
|
if rbitpos /= nbits then
|
|
r2p := target
|
|
r2p_offset := target_offset + rbitpos // limb_bits
|
|
last_nbits := nbits - rbitpos
|
|
tn := bits_to_limbs (last_nbits)
|
|
junk := lc (tp, tp_offset)
|
|
if rbitpos \\ limb_bits /= 0 then
|
|
savelimb := r2p [r2p_offset]
|
|
lshift (r2p, r2p_offset, tp, tp_offset, tn, rbitpos \\ limb_bits, carry)
|
|
rcy := carry.item
|
|
r2p [r2p_offset] := r2p [r2p_offset].bit_or (savelimb)
|
|
if rbitpos + tn * limb_bits - rbitpos \\ limb_bits < nbits then
|
|
r2p [r2p_offset + tn] := rcy
|
|
end
|
|
else
|
|
r2p.copy_data (tp, tp_offset, r2p_offset, tn)
|
|
end
|
|
if nbits \\ limb_bits /= 0 then
|
|
target [target_offset + nbits // limb_bits] := target [target_offset + nbits // limb_bits].bit_and (((0).to_natural_32.bit_not |<< (nbits \\ limb_bits)).bit_not)
|
|
end
|
|
end
|
|
end
|
|
|
|
lc (target: SPECIAL [NATURAL_32]; target_offset: INTEGER): INTEGER
|
|
local
|
|
tp: SPECIAL [NATURAL_32]
|
|
tp_offset: INTEGER
|
|
seedp: SPECIAL [NATURAL_32]
|
|
seedp_offset: INTEGER
|
|
ap: SPECIAL [NATURAL_32]
|
|
ap_offset: INTEGER
|
|
ta: INTEGER
|
|
tn: INTEGER
|
|
seedn: INTEGER
|
|
an: INTEGER
|
|
m2exp: INTEGER
|
|
bits: INTEGER
|
|
cy: INTEGER
|
|
xn: INTEGER
|
|
tmp: INTEGER
|
|
i: INTEGER
|
|
x: NATURAL_32
|
|
limb: NATURAL_32
|
|
count: INTEGER
|
|
junk: NATURAL_32
|
|
carry: CELL [NATURAL_32]
|
|
do
|
|
create carry.put (0)
|
|
m2exp := seed.m2exp
|
|
seedp := seed.seed.item
|
|
seedn := seed.seed.count
|
|
ap := seed.a.item
|
|
an := seed.a.count
|
|
ta := an + seedn + 1
|
|
tn := bits_to_limbs (m2exp)
|
|
if ta <= tn then
|
|
tmp := an + seedn
|
|
ta := tn + 1
|
|
end
|
|
create tp.make_filled (0, ta)
|
|
mul_special (tp, tp_offset, seedp, seedp_offset, seedn, ap, ap_offset, an, carry)
|
|
junk := carry.item
|
|
i := seed.cn
|
|
if i /= 0 then
|
|
add_n (tp, tp_offset, tp, tp_offset, seed.cp, 0, i, carry)
|
|
if carry.item /= 0 then
|
|
from
|
|
cy := 0
|
|
limb := 1
|
|
until
|
|
cy /= 0 or limb = 0
|
|
loop
|
|
if i >= tn then
|
|
cy := 1
|
|
else
|
|
x := tp [tp_offset + i]
|
|
limb := x + 1
|
|
tp [tp_offset + i] := limb
|
|
i := i + 1
|
|
end
|
|
end
|
|
end
|
|
end
|
|
tp [tp_offset + m2exp // limb_bits] := tp [tp_offset + m2exp // limb_bits].bit_and ((integer_to_limb (1) |<< m2exp \\ integer_to_limb (limb_bits)) - 1)
|
|
seed.seed.item.copy_data (tp, tp_offset, 0, tn)
|
|
bits := m2exp // 2
|
|
xn := bits // limb_bits
|
|
tn := tn - xn
|
|
if tn > 0 then
|
|
count := bits \\ limb_bits
|
|
if count /= 0 then
|
|
rshift (tp, tp_offset, tp, tp_offset + xn, tn, count, carry)
|
|
junk := carry.item
|
|
target.copy_data (tp, tp_offset, target_offset, xn + 1)
|
|
else
|
|
target.copy_data (tp, tp_offset + xn, target_offset, tn)
|
|
end
|
|
end
|
|
Result := (m2exp + 1) // 2
|
|
end
|
|
end
|