An implementation of a xorshift128+ RNG in shell code

xorshift128plus

#!/bin/bash
#set -x

# Function to generate a reliable seed for whatever method requires one
# Because 'date +%s' isn't entirely portable, we try other methods to get the epoch
Fn_getSeed() {
  # First we check if /dev/urandom is available.
  # We used to have a method_urandom.  /dev/urandom can generate numbers fast
  # But selecting numbers in ranges etc made for a fairly slow method
  if [ -c /dev/urandom ] && command -v od >/dev/null 2>&1; then
    # Get a string of bytes from /dev/urandom using od
    od -N 4 -A n -t uL /dev/urandom | tr -d '\n' | awk '{$1=$1+0};1'
  # Otherwise we try to get another seed
  elif date '+%s' >/dev/null 2>&1; then
    date '+%s'
  # Portable workaround based on http://www.etalabs.net/sh_tricks.html
  # We take extra steps to try to prevent accidental octal interpretation
  else
    local secsVar=$(TZ=GMT0 date +%S)
    local minsVar=$(TZ=GMT0 date +%M)
    local hourVar=$(TZ=GMT0 date +%H)
    local dayVar=$(TZ=GMT0 date +%j | sed 's/^0*//')
    local yrOffset=$(( $(TZ=GMT0 date +%Y) - 1600 ))
    local yearVar=$(( (yrOffset * 365 + yrOffset / 4 - yrOffset / 100 + yrOffset / 400 + dayVar - 135140) * 86400 ))

    printf '%s\n' "$(( yearVar + (${hourVar#0} * 3600) + (${minsVar#0} * 60) + ${secsVar#0} ))"
  fi
}

nCount="${1:-1}"
nMin="${2:-1}"
seed0=4          # RFC 1149.5
seed1="${4:-$(Fn_getSeed)}"

# Figure out our default maxRand, using 'getconf'
if [ "$(getconf LONG_BIT)" -eq 64 ]; then
  # 2^63-1
  maxRand=9223372036854775807
elif [ "$(getconf LONG_BIT)" -eq 32 ]; then
  # 2^31-1
  maxRand=2147483647
else
  # 2^15-1
  maxRand=32767
fi

nMax="${3:-maxRand}"

nRange=$(( nMax - nMin  + 1 ))

int0="${seed0}"
int1="${seed1}"
while (( nCount > 0 )); do
  # This is required for modulo de-biasing
  # We figure out the number of problematic integers to skip
  # in order to bring modulo back into uniform alignment
  # This is significantly faster than naive rejection sampling
  #skipInts=$(( ((maxRand % nMax) + 1) % nMax ))
  skipInts=$(( (maxRand - nRange) % nRange ))

  # xorshift128+ with preselected triples
  int1=$(( int1 ^ int1 << 23 ))
  int1=$(( int1 ^ int1 >> 17 ))
  int1=$(( int1 ^ int0 ))
  int1=$(( int1 ^ int0 >> 26 ))
  seed1="${int1}"
 
  # If our generated int is larger than the number of problematic ints
  # Then we can modulo it safely, otherwise drop it and generate again
  if (( (seed0 + seed1) > skipInts )); then
    #printf '%u\n' "$(( ((seed0 + seed1) % nMax) + 1 ))"
    printf '%u\n' "$(( ((seed0 + seed1) % nRange) + nMin))"
    nCount=$(( nCount - 1 ))
  fi
  
done