whrandom.sh - POSIX compliant random number generator using Wichmann-Hill method

##################################################################################################################################
# whrandom v2.31 - POSIX compliant random number generator using old or new improved Wichmann-Hill method.                       #
# Handshake:                                                                                                                     #
# $WH_MIN contains the minimum value for the random number (default: 0).                                                         #
# $WH_MAX contains the maximum value for the random number. The limit is 4 billion. If rejection sampling is used the process    #
#  will slow down significantly if the max value exceeds 999999999 (default: 4000000000).                                        #
# $WH_LOOP tells how many random numbers to be generated. Set it to "U" for an unlimited random number generation (default: 1)   #
# $WH_DEL contains the delimiter used to separate multiple values (default: \n (new line)).                                      #
# $WH_ALGO tells which algorithm will be used - WH_OLD, WH_NEW or WH_NEW64() (default: WH_NEW).                                  #
# $WH_SCALE sets the scaling mode if $WH_MIN and/or $WH_MAX is used. Leave empty to trigger the slower, but better rejection     #
#  sampling method. Set it to anything to trigger a simpler, but slightly faster rounding method (default: fast).                #
# Usage:                                                                                                                         #
# Simply call WH_RANDOM() to generate one or more random numbers. These numbers will be written to stdout.                       #
# The old Wichmann-Hill algorrithm requires 3, the new one 4 seeds each between 1 and 30000. By default, whrandom automatically  #
# creates 4 unique seeds using the millisecond counter of the system clock. Alternatively, you can call WH_RANDOM() with up to 4 #
# integers. They will be checked and exchanged with automatically generated ones if they don't fit the needs. Note that if you   #
# use 4 fixed seeds whrandom will always create an identical chain of random numbers as long as you use the same seeds.          #
##################################################################################################################################

WH_MIN=0
WH_MAX=4000000000
WH_LOOP=1
WH_DEL="\n"
WH_ALGO=WH_NEW
WH_SCALE=fast

WH_CHECK(){ # checks if a value is given, if this value is a number and if it's between 1 and 30000. If so it returns that very
#             number and if not it generates a suitable one using the system clock's nanosecond counter.
# IN:  $1 should contain an integer between 1 and 30000.
# OUT: Prints the fitting number to stdout.
    if [ -n "${1##*[!0-9]*}" ] && [ $1 -gt 0 ] && [ $1 -le 30000 ]; then
        printf $1
    else # fetch 5 digits of nanoseconds; removes leading zeros & adds 1 to prevent zero; if value is >30000 divide by 4
        WH_TMP=$(date +%5N)
        [ $((WH_TMP=${WH_TMP#${WH_TMP%%[1-9]*}}+1)) -gt 30000 ] && WH_TMP=$(($WH_TMP/4))
        printf $WH_TMP
    fi
}

WH_OLD(){ # generates a random number using the older 16bit integer based algorithm.
# IN:  $WH_S1, $WH_S2 & $WH_S3 contain the values (seeds) needed by the algorithm for computation.
# OUT: $WH_RND contains the random number to be handed over to the core function for further operations.
    [ $((WH_S1=171*($WH_S1%177)-2*($WH_S1/177))) -lt 0 ] && WH_S1=$(($WH_S1+30269))
    [ $((WH_S2=172*($WH_S2%176)-35*($WH_S2/176))) -lt 0 ] && WH_S2=$(($WH_S2+30307))
    [ $((WH_S3=170*($WH_S3%178)-63*($WH_S3/178))) -lt 0 ] && WH_S3=$(($WH_S3+30323))
    WH_RND=$(($WH_S1*1000000000/30269+$WH_S2*1000000000/30307+$WH_S3*1000000000/30323))
}

WH_NEW(){ # generates a random number using the newer 32bit integer based algorithm.
# IN:  $WH_S1, $WH_S2, $WH_S3 & $WH_S4 contain the values (seeds) needed by the algorithm for computation.
# OUT: $WH_RND contains the random number to be handed over to the core function for further operations.
    [ $((WH_S1=11600*($WH_S1%185127)-10379*($WH_S1/185127))) -lt 0 ] && WH_S1=$(($WH_S1+2147483579))
    [ $((WH_S2=47003*($WH_S2%45688)-10479*($WH_S2/45688))) -lt 0 ] && WH_S2=$(($WH_S2+2147483543))
    [ $((WH_S3=23000*($WH_S3%93368)-19423*($WH_S3/93368))) -lt 0 ] && WH_S3=$(($WH_S3+2147483423))
    [ $((WH_S4=33000*($WH_S4%65075)-8123*($WH_S4/65075))) -lt 0 ] && WH_S4=$(($WH_S4+2147483123))
    WH_RND=$(($WH_S1*1000000000/2147483579+$WH_S2*1000000000/2147483543+$WH_S3*1000000000/2147483423+$WH_S4*1000000000/2147483123))
}

WH_NEW64(){ # same result as WH_NEW(), but uses 64bit integer based algorithm with less steps & thus being slightly faster.
    WH_S1=$((11600*$WH_S1%2147483579))
    WH_S2=$((47003*$WH_S2%2147483543))
    WH_S3=$((23000*$WH_S3%2147483423))
    WH_S4=$((33000*$WH_S4%2147483123))
    WH_RND=$(($WH_S1*1000000000/2147483579+$WH_S2*1000000000/2147483543+$WH_S3*1000000000/2147483423+$WH_S4*1000000000/2147483123))
}

WH_RANDOM(){ # core function: prepares seeds and generates the output.
# IN:  $1, $2, $3 and $4 can contain integers between 1 and 30000 to be used as seeds. These values will be checked and - if not
#       suitable - replaced with auto generated ones.
#      $WH_LOOP tells how many numbers will be generated. "U" will generate a never ending chain of random numbers.
#      $WH_RND contains the random number which was handed over by the algorithm for further operations.
#      $WH_MIN contains the minimum value which will not be underrun by the random number.
#      $WH_MAX contains the maximum value which will not be exceeded by the random number. 0 leaves the result unaltered.
#      $WH_DEL contains the delimiter which is placed after each value.
#      $WH_SCALE tells which method is used to scale the random number. If left empty, rejection sampling is used which is slower, but
#       with better diffusion & no bias. If set to anything, it uses a slightly faster rounding method with a little worse diffusion.
# OUT: $WH_SEEDS contains the original seeds used for the 1st iteration.
#      $WH_RNG is the range from $WH_MIN to $WH_MAX.
#      $WH_LMT is the largest multiple of $WH_RNG which is still < or = 1000000000.
#      $WH_LO contains the lower block (everything below 1000000000 if $WH_MAX exceeds 999999999).
#      $WH_HI contains the block index to decide which range above 999999999 is used in addition to $WH_LO.
#      $WH_RND contains the current random number
#      Also prints the random numbers with delimiter appended to stdout.
    WH_TMP=1
    [ -z $WH_MAX ] || [ $WH_MAX -le 0 ] || [ $WH_MAX -gt 4000000000 ] && WH_MAX=4000000000
    [ $WH_MIN -ge $WH_MAX ] && return 1
    set -- $(WH_CHECK "$1") $(WH_CHECK "$2") $(WH_CHECK "$3") $(WH_CHECK "$4")
    WH_S1=$1 WH_S2=$2 WH_S3=$3 WH_S4=$4
    WH_SEEDS="$WH_S1 $WH_S2 $WH_S3 $WH_S4"
    while [ $WH_LOOP = U ] || [ $((WH_TMP=$WH_TMP+1)) -le $(($WH_LOOP+1)) ]; do
        $WH_ALGO
        if [ -z $WH_SCALE ]; then
            if [ $((WH_RNG=$WH_MAX-$WH_MIN+1)) -gt 1000000000 ]; then
                WH_LO=$(($WH_RND%1000000000))
                $WH_ALGO
                WH_HI=$(($WH_RND%(($WH_RNG+999999999)/1000000000)))
                [ $((WH_RND=$WH_HI*1000000000+$WH_LO)) -lt $WH_RNG ] && printf $(($WH_RND+$WH_MIN))"$WH_DEL" || WH_TMP=$(($WH_TMP-1))
            else
                WH_LMT=$((1000000000/$WH_RNG*$WH_RNG))
                [ $((WH_RND=$WH_RND%1000000000)) -lt $WH_LMT ] && printf $(($WH_RND%$WH_RNG+$WH_MIN))"$WH_DEL" || WH_TMP=$(($WH_TMP-1))
            fi
        elif [ -z $WH_MAX ] || [ $WH_MAX -ge 4000000000 ]; then
            [ $WH_RND -ge $WH_MIN ] && printf $WH_RND"$WH_DEL" || WH_TMP=$(($WH_TMP-1))
        else
            printf $(((($WH_RND%1000000000*($WH_MAX-$WH_MIN)+555555555)/1000000000)+$WH_MIN))"$WH_DEL"
        fi
    done
}