This is how I implemented it. The script contains comments, test vectors and debug instructions. The actual code is quite simple.
I hope that this is useful to someone.
Code: Select all
:global getRandom64 do={
# Parameters:
# none
# Return value:
# A 64 bit pseudo random number obtained from the xorshift64* prng.
# The number lies in the range [1, 2**64). Note that it will NEVER be 0.
# The number is returned as a ROS array of integers { t=top_32_bits; b=bottom_32_bits }
#
# xorshift64* is defined in http://xorshift.di.unimi.it/
# It is the smallest state 64bit prng that is recommended especially for memory limited devices.
# As well as passing all the statistical tests (except MatrixRank test) it has the property that even on the first iteration it
# produces numbers with each bit =1 has 0.5 probability (see
# chapter 7.2 "Escaping Zeroland" in http://vigna.di.unimi.it/ftp/papers/xorshift.pdf )
# End chapter 5 (page 14) "suggests that for these generators it is better to extract the lower bits,
# rather than the high bits, when just a subsequence is needed."
#
# The algorithm is (from the above references):
# uint64_t x; /* The state must be seeded with a nonzero value. */
# uint64_t next() {
# x ^= x >> 12; // a
# x ^= x << 25; // b
# x ^= x >> 27; // c
# return x * 2685821657736338717LL;
# }
# http://www.javascripter.net/math/calculators/primefactorscalculator.htm
#
# ROS does not have 64 bit unsigned integers so we need to emulate using a pair of
# 32bit unsigned integers which are contained in the lower 32 bits of two ROS 64 bit signed integers.
# Care must be taken when doing addition and multiplication to avoid hitting the sign bit in the 64bit ROS variable.
# Fortunately, the xorshift64* algorithm only requires multiplication by the constant 2685821657736338717LL which has
# a hex representation of 2545f491 4f6cdd1d where the space indicates the 32 bit boundary. This means that the largest multiplication
# result within a "32bit half" is 4f6cdd1d * (2**32 -1 ) < 4f6cdd1d00000000. Thus overflow will never occur.
# Note that additions are done only using 32-bit halves thus overflow can never occur.
# The following test vectors were obtained from an independent implementation
# Initial seed (t, b) xorshift64* result
# 0,0 0
# 0,1 0x47e4ce4b896cdd1d
# 0,2 0x8fc99c9712d9ba3a
# 1,1 0x84e9495afd4c80bd
# 0x21212121,0x32323232 0xff3ddfa115d6198d
# 0xfefefefe,0xcacacaca 0xc67bafe0b4bc30cf
#
# In this implementation a 64 bit unsigned variable x is represented as the pair xt ("x top") and xb ("x bottom"). The
# most significant 32 bits of x are in xt whilst the least significant are in xb
:global Random64Statet
:global Random64Stateb
# Uncomment following two lines to debug
#:set Random64Statet [:tonum $1]
#:set Random64Stateb [:tonum $2]
# Define the multiplier 2685821657736338717LL:
:local ct 0x2545F491
:local cb 0x4F6CDD1D
# Define the multiplier used for calculating modulo 2**32
:local mask32 0xFFFFFFFF
# Do x ^= x >> 12
:set Random64Stateb ($Random64Stateb ^ ((($Random64Stateb >> 12) | ($Random64Statet << 20)) & $mask32))
:set Random64Statet ($Random64Statet ^ ($Random64Statet >> 12))
# Do x ^= x << 25
:set Random64Statet ($Random64Statet ^ ((($Random64Statet << 25) | ($Random64Stateb >> 7)) & $mask32))
:set Random64Stateb ($Random64Stateb ^ (($Random64Stateb << 25) & $mask32))
# Do x ^= x >> 27
:set Random64Stateb ($Random64Stateb ^ ((($Random64Stateb >> 27) | ($Random64Statet << 5)) & $mask32))
:set Random64Statet ($Random64Statet ^ ($Random64Statet >> 27))
# Calculate result = x * 2685821657736338717LL
# Use the formulae (we work mod 2**64):
# x * c = (xt * 2**32 + xb) * (ct * 2**32 + cb)
# = xt * ct * 2**64 + (xt * cb + xb * ct) * 2**32 + xb * cb
# = (xt * cb + xb * ct) * 2**32 + xb * cb # because we are working mod 2**64
# Note that the top half of xb * cb must be carried to the top half of the result
# and the top half of (xt * cb + xb * ct) can be discarded.
:local resultt (((($Random64Statet * $cb) & $mask32) + (($Random64Stateb * $ct) & $mask32) + (($Random64Stateb * $cb) >> 32)) & $mask32)
:local resultb (($Random64Stateb * $cb) & $mask32)
# Uncomment following line to debug
#:global convertToHex
#:put "$[$convertToHex $resultt] $[$convertToHex $resultb]"
:return { t=$resultt; b=$resultb}
}