[PATCH 1/1] zsh/random

[Clinton Bunch <cdb_zsh@xxxxxxxxxxx>]

diff --git a/Completion/Zsh/Type/_module_math_func 
b/Completion/Zsh/Type/_module_math_func
index 5044bdf4c..c997bde8b 100644
--- a/Completion/Zsh/Type/_module_math_func
+++ b/Completion/Zsh/Type/_module_math_func
@@ -2,7 +2,7 @@

 local mod
 local -a funcs alts
-local -a modules=( example mathfunc system )
+local -a modules=( example mathfunc system random)

 for mod in $modules; do
   funcs=( ${${${(f)"$(zmodload -Fl zsh/$mod 2>/dev/null)"}:#^+f:*}##+f:} )
diff --git a/Doc/Makefile.in b/Doc/Makefile.in
index d9be182e9..1197e74aa 100644
--- a/Doc/Makefile.in
+++ b/Doc/Makefile.in
@@ -68,7 +68,7 @@ Zsh/mod_hlgroup.yo Zsh/mod_langinfo.yo \
 Zsh/mod_ksh93.yo Zsh/mod_mapfile.yo Zsh/mod_mathfunc.yo \
 Zsh/mod_nearcolor.yo Zsh/mod_newuser.yo \
 Zsh/mod_parameter.yo Zsh/mod_pcre.yo Zsh/mod_private.yo \
-Zsh/mod_regex.yo Zsh/mod_sched.yo Zsh/mod_socket.yo \
+Zsh/mod_regex.yo Zsh/mod_random.yo Zsh/mod_sched.yo Zsh/mod_socket.yo \
 Zsh/mod_stat.yo  Zsh/mod_system.yo Zsh/mod_tcp.yo \
 Zsh/mod_termcap.yo Zsh/mod_terminfo.yo \
 Zsh/mod_watch.yo \
diff --git a/Doc/Zsh/mod_random.yo b/Doc/Zsh/mod_random.yo
new file mode 100644
index 000000000..17f6e5ff6
--- /dev/null
+++ b/Doc/Zsh/mod_random.yo
@@ -0,0 +1,57 @@
+COMMENT(!MOD!zsh/random
+Some High-quality randomness parameters and functions.
+!MOD!)
+The tt(zsh/random) module gets random data from the kernel random pool. 
If no
+kernel random pool can be found, the module will not load.
+
+subsect(Parameters)
+
+startitem()
+vindex(SRANDOM)
+item(tt(SRANDOM)) (
+A random positive 32-bit integer between 0 and 4,294,967,295. This 
parameter
+is read-only. The name was chosen for compatibility with Bash and to
+distinguish it from tt(RANDOM) which has a documented repeatable behavior.
+)
+enditem()
+
+subsect(Math Functions)
+
+startitem()
+item(tt(zrand_float+LPAR()RPAR())) (
+Returns a random floating point number between 0 and 1.
+)
+enditem()
+
+startitem()
+item(tt(zrand_int)+LPAR()tt(upper), tt(lower), tt(inclusive)RPAR()) (
+Returns a random integer between tt(lower) and tt(upper). All 
parameters are
+optional.  If none are specified it is equivalent to
+tt(SRANDOM).
+
+tt(upper) is the upper bound on the resultant number and defaults to
+4,294,967,295.
+
+tt(lower) is the lower bound and defaults to 0.
+
+The defaults of these two arguments are also the maximum and minimum to 
which
+either can be set.
+
+tt(inclusive) is a flag that controls whether the result is ever equal to
+tt(upper).  By default it is not. If this argument is set to a non-zero 
value
+then it may be.
+
+This is to facilitate a construct like tt($a[zrand_int($#a)+1]) rather
+than tt($a[zrand_int+LPAR()$#a-1+RPAR()+1]).
+For example, if $#a is 16, you would use tt(zrand_int+LPAR()16RPAR()) 
which has
+16 possible return values 0-15.  Because the function can return zero, 
in order
+to use it as an array index from 1-16 you need to add one.  It would
+be an array index range error for it to also potentially return 16 
($#a). You
+could, however, use the construct tt(zrand_int+LPAR()16,1,1+RPAR()) 
instead of
+adding 1 to achieve the same result, but it is more verbose.
+
+Most statistics algorithms seem to also expect 0 to tt(upper)-1, so 
this was
+deemed the most commonly desired case and chosen as the default.
+)
+enditem()
+
diff --git a/Src/Modules/random.c b/Src/Modules/random.c
new file mode 100644
index 000000000..b45483afe
--- /dev/null
+++ b/Src/Modules/random.c
@@ -0,0 +1,324 @@
+/*
+ * random.c - module to access kernel random sources.
+ *
+ * This file is part of zsh, the Z shell.
+ *
+ * Copyright (c) 2022 Clinton Bunch
+ * All rights reserved.
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and to distribute modified versions of this software for any
+ * purpose, provided that the above copyright notice and the following
+ * two paragraphs appear in all copies of this software.
+ *
+ * In no event shall Clinton Bunch or the Zsh Development Group be liable
+ * to any party for direct, indirect, special, incidental, or consequential
+ * damages arising out of the use of this software and its documentation,
+ * even if Clinton Bunch and the Zsh Development Group have been advised of
+ * the possibility of such damage.
+ *
+ * Clinton Bunch and the Zsh Development Group specifically disclaim any
+ * warranties, including, but not limited to, the implied warranties of
+ * merchantability and fitness for a particular purpose.  The software
+ * provided hereunder is on an "as is" basis, and Clinton Bunch and the
+ * Zsh Development Group have no obligation to provide maintenance,
+ * support, updates, enhancements, or modifications.
+ *
+ */
+
+#include "random.mdh"
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <stdio.h>
+
+#include <stdbool.h>
+#include <stdint.h>
+
+#ifdef HAVE_SYS_RANDOM_H
+#include <sys/random.h>
+#endif
+
+/* Simplify select URANDOM specific code */
+#if !defined(HAVE_ARC4RANDOM_BUF) && !defined(HAVE_GETRANDOM)
+#define USE_URANDOM
+#endif
+
+/* buffer to pre-load integers for SRANDOM to lessen the context 
switches */
+uint32_t rand_buff[8];
+static int buf_cnt = -1;
+
+#ifdef USE_URANDOM
+/* File descriptor for /dev/urandom */
+int randfd = -1;
+#endif /* USE_URANDOM */
+
+static zlong get_srandom(UNUSED(Param p));
+
+/**/
+ssize_t
+getrandom_buffer(void *buf, size_t len)
+{
+    ssize_t ret;
+    size_t  val     = 0;
+    uint8_t *bufptr = buf;
+
+    do {
+       errno = 0;
+#ifdef HAVE_ARC4RANDOM_BUF
+       arc4random_buf(buf,len);
+       ret = len;
+#elif defined(HAVE_GETRANDOM)
+       ret=getrandom(bufptr,(len - val),0);
+#else
+       ret=read(randfd,bufptr,(len - val));
+#endif
+       if (ret < 0) {
+           if (errno != EINTR || errflag || retflag || breaks || 
contflag) {
+               zwarn("Unable to get random data: %e.", errno);
+               return -1;
+           }
+       }
+       bufptr += ret;
+       val    += ret;
+    } while (ret < len);
+    return ret;
+}
+
+/*
+ * Generate count number of random 32-bit integers between 0 and max-1
+ * Got this algorithm from
+ * https://lemire.me/blog/2016/06/30/fast-random-shuffling/
+ * adapting the public domain code.
+ *
+ */
+
+/**/
+void
+get_bound_random_buffer(uint32_t *buffer, size_t count, uint32_t max)
+{
+    uint64_t multi_result;
+    uint32_t threshold;
+    uint32_t leftover;
+
+    size_t i; /* loop counter */
+
+    getrandom_buffer((void*) buffer, count*sizeof(uint32_t));
+    if (max == UINT32_MAX)
+       return;
+
+    for(i=0;i<count;i++) {
+       multi_result = ((uint64_t) buffer[i]) * (uint64_t) max;
+       leftover = (uint32_t) multi_result;
+
+       /*
+        * The following if statement should (according to Google's Gemini)
+        * only be executed with a probability of 1/2**32 or 2.33e-10
+        */
+       if(leftover < max) {
+           threshold= -max % max;
+           while (leftover < threshold) {
+               uint32_t j=get_srandom(NULL);
+               multi_result=(uint64_t) j * (uint64_t) max;
+               leftover= (uint32_t) multi_result;
+           }
+       }
+       buffer[i]=multi_result >> 32;
+    }
+}
+
+/*
+ * Provides for the SRANDOM parameter and returns an unsigned 32-bit random
+ * integer.
+ */
+
+/**/
+static zlong
+get_srandom(UNUSED(Param pm)) {
+
+    if(buf_cnt <= 0) {
+       getrandom_buffer((void*) rand_buff,sizeof(rand_buff));
+       buf_cnt=sizeof(rand_buff)/sizeof(rand_buff[0]);
+    }
+    return rand_buff[--buf_cnt];
+}
+
+/*
+ * Implements math function zrand_int takes 0 to 3 arguments an upper 
bound,
+ * a lower bound and a flag as to whether the range is inclusive or 
not.  The
+ * default is exclusive.  If neither upper or lower is specified this is no
+ * different than SRANDOM.
+ */
+
+/**/
+static mnumber
+math_zrand_int(UNUSED(char *name), int argc, mnumber *argv, UNUSED(int id))
+{
+    mnumber ret;
+    uint32_t i;
+    zlong lower=0, upper=UINT32_MAX,incl=0, diff;
+
+    ret.type = MN_INTEGER;
+
+    switch (argc) {
+       case 0: ret.u.l=get_srandom(NULL);
+               return ret;
+               break;
+       case 3: incl = (argv[2].u.l != 0)?1:0;
+       case 2: lower = argv[1].u.l;
+       case 1: upper = argv[0].u.l;
+       default: diff = upper-lower+incl;
+    }
+
+    if (lower < 0 || lower >= UINT32_MAX) {
+       zwarn("Lower bound (%z) out of range: 0-4294967295",lower);
+    } else if (upper < lower) {
+       zwarn("Upper bound (%z) must be greater than Lower Bound 
(%z)",upper,lower);
+    } else if (upper < 0 || upper >= UINT32_MAX) {
+       zwarn("Upper bound (%z) out of range: 0-4294967295",upper);
+    }
+
+    if ( diff == 0 ) {
+       ret.u.l=upper; /* still not convinced this shouldn't be an error. */
+    } else {
+       get_bound_random_buffer(&i,1,(uint32_t) diff);
+       ret.u.l=i+lower;
+    }
+    return ret;
+}
+
+/*
+ * Implements the mathfunc zrand_float and returns a random floating-point
+ * number between 0 and 1.
+ *
+ */
+
+/**/
+static mnumber
+math_zrand_float(UNUSED(char *name), UNUSED(int argc), UNUSED(mnumber 
*argv),
+            UNUSED(int id))
+{
+    mnumber ret;
+    double r;
+
+    r = random_real();
+    if (r < 0) {
+       zwarnnam(name, "Failed to get sufficient random data.");
+    }
+    ret.type = MN_FLOAT;
+    ret.u.d = r;
+
+    return ret;
+}
+
+static const struct gsu_integer srandom_gsu =
+{ get_srandom, nullintsetfn, stdunsetfn };
+
+static struct paramdef patab[] = {
+    {"SRANDOM", PM_INTEGER | PM_READONLY_SPECIAL | PM_HIDEVAL, NULL,
+           &srandom_gsu, NULL, NULL, NULL},
+};
+
+static struct mathfunc mftab[] = {
+    NUMMATHFUNC("zrand_float", math_zrand_float, 0, 0, 0),
+    NUMMATHFUNC("zrand_int", math_zrand_int, 0, 3, 0),
+};
+
+static struct features module_features = {
+    NULL, 0,
+    NULL, 0,
+    mftab, sizeof(mftab)/sizeof(*mftab),
+    patab, sizeof(patab)/sizeof(*patab),
+    0
+};
+
+/**/
+int
+setup_(UNUSED(Module m))
+{
+#ifdef USE_URANDOM
+    /* Check for the existence of /dev/urandom */
+
+    struct stat st;
+
+    errno=0;
+    if (stat("/dev/urandom",&st) < 0) {
+       zwarn("Error getting kernel random pool: %e.", errno);
+       return 1;
+    }
+
+    errno=0;
+    if (!(S_ISCHR(st.st_mode)) ) {
+       zwarn("Error getting kernel random pool: %e.", errno);
+       return 1;
+    }
+#endif /* USE_URANDOM */
+    return 0;
+}
+
+/**/
+int
+features_(Module m, char ***features)
+{
+    *features = featuresarray(m, &module_features);
+    return 0;
+}
+
+/**/
+int
+enables_(Module m, int **enables)
+{
+    return handlefeatures(m, &module_features, enables);
+}
+
+/**/
+int
+boot_(Module m)
+{
+#ifdef USE_URANDOM
+    /*
+     * Open /dev/urandom.  Here because of a weird issue on HP-UX 11.31
+     * When opening in setup_ open returned 0.  In boot_, it returns
+     * an unused file descriptor. Decided against ifdef HPUX as it works
+     * here just as well for other platforms.
+     *
+     */
+
+    int tmpfd=-1;
+
+    errno=0;
+    if ((tmpfd = open("/dev/urandom", O_RDONLY)) < 0) {
+       zwarn("Could not access kernel random pool: %e.",errno);
+       return 1;
+    }
+    randfd = movefd(tmpfd);
+    addmodulefd(randfd,FDT_MODULE);
+    if (randfd < 0) {
+       zwarn("Could not access kernel random pool.");
+       return 1;
+    }
+#endif /* USE_URANDOM */
+    return 0;
+}
+
+/**/
+int
+cleanup_(Module m)
+{
+    return setfeatureenables(m, &module_features, NULL);
+}
+
+/**/
+int
+finish_(UNUSED(Module m))
+{
+#ifdef USE_URANDOM
+    if (randfd >= 0)
+       zclose(randfd);
+#endif /* USE_URANDOM */
+    return 0;
+}
+
diff --git a/Src/Modules/random.mdd b/Src/Modules/random.mdd
new file mode 100644
index 000000000..7a75f29ff
--- /dev/null
+++ b/Src/Modules/random.mdd
@@ -0,0 +1,7 @@
+name=zsh/random
+link=either
+load=yes
+
+autofeatures="p:SRANDOM f:zrand_float f:zrand_int"
+
+objects="random.o random_real.o"
diff --git a/Src/Modules/random_real.c b/Src/Modules/random_real.c
new file mode 100644
index 000000000..900360414
--- /dev/null
+++ b/Src/Modules/random_real.c
@@ -0,0 +1,213 @@
+/*  This file contains code under different copyrights separated by */
+/* ====@@@@@=== */
+
+/*
+ * random_real.c - module to access kernel random sources.
+ *
+ * This file is part of zsh, the Z shell.
+ *
+ * Copyright (c) 2022 Clinton Bunch
+ * All rights reserved.
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and to distribute modified versions of this software for any
+ * purpose, provided that the above copyright notice and the following
+ * two paragraphs appear in all copies of this software.
+ *
+ * In no event shall Clinton Bunch or the Zsh Development Group be liable
+ * to any party for direct, indirect, special, incidental, or consequential
+ * damages arising out of the use of this software and its documentation,
+ * even if Clinton Bunch and the Zsh Development Group have been advised of
+ * the possibility of such damage.
+ *
+ * Clinton Bunch and the Zsh Development Group specifically disclaim any
+ * warranties, including, but not limited to, the implied warranties of
+ * merchantability and fitness for a particular purpose.  The software
+ * provided hereunder is on an "as is" basis, and Clinton Bunch and the
+ * Zsh Development Group have no obligation to provide maintenance,
+ * support, updates, enhancements, or modifications.
+ *
+ */
+
+#include "random.mdh"
+
+#include <math.h>
+#include <stdint.h>
+#include <errno.h>
+
+
+/* Count the number of leading zeros, hopefully in gcc/clang by HW
+ * instruction */
+#if defined(__GNUC__) || defined(__clang__)
+#define clz64(x) __builtin_clzll(x)
+#else
+#define clz64(x) _zclz64(x)
+
+/**/
+int
+_zclz64(uint64_t x) {
+    int n = 0;
+
+    if (x == 0)
+       return 64;
+
+    if (!(x & 0xFFFFFFFF00000000ull)) {
+       n+=32;
+       x<<=32;
+    }
+    if (!(x & 0xFFFF000000000000ull)) {
+       n+=16;
+       x<<=16;
+    }
+    if (!(x & 0xFF00000000000000ull)) {
+       n+=8;
+       x<<=8;
+    }
+    if (!(x & 0xF000000000000000ull)) {
+       n+=4;
+       x<<=4;
+    }
+    if (!(x & 0xC000000000000000ull)) {
+       n+=2;
+       x<<=1;
+    }
+    if (!(x & 0x8000000000000000ull)) {
+       n+=1;
+    }
+    return n;
+}
+#endif /* __GNU_C__ or __clang__ */
+
+/**/
+uint64_t
+random_64bit(void) {
+    uint64_t r;
+
+    if(getrandom_buffer(&r,sizeof(r)) < 0) {
+       zwarn("zsh/random: Can't get sufficient random data.");
+       return 1; /* 0 will cause loop */
+    }
+
+    return r;
+}
+
+/* ====@@@@@=== */
+/* Following code is under the below copyright, despite changes for error
+ * handling and removing GCCisms */
+
+/*-
+ * Copyright (c) 2014 Taylor R. Campbell
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
+ */
+
+/*
+ * Uniform random floats: How to generate a double-precision
+ * floating-point numbers in [0, 1] uniformly at random given a uniform
+ * random source of bits.
+ *
+ * See <http://mumble.net/~campbell/2014/04/28/uniform-random-float>
+ * for explanation.
+ *
+ * Updated 2015-02-22 to replace ldexp(x, <constant>) by x * ldexp(1,
+ * <constant>), since glibc and NetBSD libm both have slow software
+ * bit-twiddling implementations of ldexp, but GCC can constant-fold
+ * the latter.
+ */
+
+/*
+ * random_real: Generate a stream of bits uniformly at random and
+ * interpret it as the fractional part of the binary expansion of a
+ * number in [0, 1], 0.00001010011111010100...; then round it.
+ */
+
+/**/
+double
+random_real(void)
+{
+       int exponent = 0;
+       uint64_t significand = 0;
+       uint64_t r = 0;
+       unsigned shift;
+
+       /*
+        * Read zeros into the exponent until we hit a one; the rest
+        * will go into the significand.
+        */
+       while (significand == 0) {
+               exponent -= 64;
+
+               /* Get random_64bit and check for error */
+               errno = 0;
+               significand = random_64bit();
+               if (errno)
+                   return -1;
+               /*
+                * If the exponent falls below -1074 = emin + 1 - p,
+                * the exponent of the smallest subnormal, we are
+                * guaranteed the result will be rounded to zero.  This
+                * case is so unlikely it will happen in realistic
+                * terms only if random_64bit is broken.
+                */
+               if (exponent < -1074)
+                       return 0;
+       }
+
+       /*
+        * There is a 1 somewhere in significand, not necessarily in
+        * the most significant position.  If there are leading zeros,
+        * shift them into the exponent and refill the less-significant
+        * bits of the significand.  Can't predict one way or another
+        * whether there are leading zeros: there's a fifty-fifty
+        * chance, if random_64bit is uniformly distributed.
+        */
+       shift = clz64(significand);
+       if (shift != 0) {
+
+               errno = 0;
+               r = random_64bit();
+               if (errno)
+                   return -1;
+
+               exponent -= shift;
+               significand <<= shift;
+               significand |= (r >> (64 - shift));
+       }
+
+       /*
+        * Set the sticky bit, since there is almost surely another 1
+        * in the bit stream.  Otherwise, we might round what looks
+        * like a tie to even when, almost surely, were we to look
+        * further in the bit stream, there would be a 1 breaking the
+        * tie.
+        */
+       significand |= 1;
+
+       /*
+        * Finally, convert to double (rounding) and scale by
+        * 2^exponent.
+        */
+       return ldexp((double)significand, exponent);
+}
+/* ====@@@@@=== */
diff --git a/Test/V15random.ztst b/Test/V15random.ztst
new file mode 100644
index 000000000..889324380
--- /dev/null
+++ b/Test/V15random.ztst
@@ -0,0 +1,120 @@
+%prep
+
+  if ! zmodload zsh/mathfunc 2>/dev/null; then
+    ZTST_unimplemented="can't load the zsh/mathfunc module for testing"
+  fi
+
+  print -ru $ZTST_fd 'This test may take two seconds...'
+
+  function calc_chi2() {
+# Calculate Chi Squared
+       integer n=$1
+       for ((k=1;k<=n;k++))
+       do
+               ((delta=${observed[$k]}-expected))
+               #echo "$k=$delta \c"
+               ((chi2+=delta**2/expected))
+       done
+       #echo
+  }
+  function ligf() {
+       float -E S=$1 Z=$2
+       float -gE RET
+       if (( Z < 0.0 )); then
+               RET=0.0
+       fi
+       float -E Sc=$((1.0/S))
+       ((Sc *= Z**S ))
+       ((Sc *= exp(-1*Z) ))
+       float -E Sum=1.0 Nom=1.0 Denom=1.0
+       integer I
+       # 200 iterations seems like enough for an approximation of an 
infinate
+       # series
+       for (( I=0;I<200;I++))
+       do
+               ((Nom *= Z))
+               ((S++))
+               ((Denom *= S))
+               ((Sum += (Nom/Denom) ))
+       done
+       RET=$((Sum * Sc))
+  }
+#echo $observed
+  function calc_p() {
+       #echo "chi2=$chi2"
+       float chi2=0 cdf=0
+        float RET
+       calc_chi2 "$n"
+       ligf $(((n-1)*0.5)) $((chi2*0.5))
+       ((cdf=$RET/gamma((n-1)*0.5) ))
+       ((p=1-cdf))
+  }
+%test
+  fail=0
+  for ((o=0;o<100;o++))
+  do
+      n=20
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( zrand_int(n,1,1)  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:random integer $samples samples between 1-$n
+
+  fail=0
+  unset observed
+  for ((o=0;o<100;o++))
+  do
+      n=20
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( ceil(zrand_float()*n)  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:random float $samples samples
+
+  fail=0
+  unset observed
+  for ((o=0;o<100;o++))
+  do
+      n=16
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( (SRANDOM & 0x0F)+1  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:SRANDOM mod 16 $samples samples
+
diff --git a/configure.ac b/configure.ac
index 78621042d..35f641cde 100644
--- a/configure.ac
+++ b/configure.ac
@@ -636,6 +636,7 @@ fi
 AC_CHECK_HEADERS(sys/time.h sys/times.h sys/select.h termcap.h termio.h \
                 termios.h sys/param.h sys/filio.h string.h memory.h \
                 limits.h fcntl.h libc.h sys/utsname.h sys/resource.h \
+                 sys/random.h \
                 locale.h errno.h stdio.h stdarg.h varargs.h stdlib.h \
                 unistd.h sys/capability.h \
                 utmp.h utmpx.h sys/types.h pwd.h grp.h poll.h sys/mman.h \
@@ -1292,6 +1293,7 @@ AC_CHECK_FUNCS(strftime strptime mktime timelocal \
               cygwin_conv_path \
               nanosleep \
               srand_deterministic \
+               getrandom arc4random_buf \
               setutxent getutxent endutxent getutent)
 AC_FUNC_STRCOLL

diff --git a/Completion/Zsh/Type/_module_math_func b/Completion/Zsh/Type/_module_math_func
index 5044bdf4c..c997bde8b 100644
--- a/Completion/Zsh/Type/_module_math_func
+++ b/Completion/Zsh/Type/_module_math_func
@@ -2,7 +2,7 @@
 
 local mod
 local -a funcs alts
-local -a modules=( example mathfunc system )
+local -a modules=( example mathfunc system random)
 
 for mod in $modules; do
   funcs=( ${${${(f)"$(zmodload -Fl zsh/$mod 2>/dev/null)"}:#^+f:*}##+f:} )
diff --git a/Doc/Makefile.in b/Doc/Makefile.in
index d9be182e9..1197e74aa 100644
--- a/Doc/Makefile.in
+++ b/Doc/Makefile.in
@@ -68,7 +68,7 @@ Zsh/mod_hlgroup.yo Zsh/mod_langinfo.yo \
 Zsh/mod_ksh93.yo Zsh/mod_mapfile.yo Zsh/mod_mathfunc.yo \
 Zsh/mod_nearcolor.yo Zsh/mod_newuser.yo \
 Zsh/mod_parameter.yo Zsh/mod_pcre.yo Zsh/mod_private.yo \
-Zsh/mod_regex.yo Zsh/mod_sched.yo Zsh/mod_socket.yo \
+Zsh/mod_regex.yo Zsh/mod_random.yo Zsh/mod_sched.yo Zsh/mod_socket.yo \
 Zsh/mod_stat.yo  Zsh/mod_system.yo Zsh/mod_tcp.yo \
 Zsh/mod_termcap.yo Zsh/mod_terminfo.yo \
 Zsh/mod_watch.yo \
diff --git a/Doc/Zsh/mod_random.yo b/Doc/Zsh/mod_random.yo
new file mode 100644
index 000000000..17f6e5ff6
--- /dev/null
+++ b/Doc/Zsh/mod_random.yo
@@ -0,0 +1,57 @@
+COMMENT(!MOD!zsh/random
+Some High-quality randomness parameters and functions.
+!MOD!)
+The tt(zsh/random) module gets random data from the kernel random pool. If no
+kernel random pool can be found, the module will not load.
+
+subsect(Parameters)
+
+startitem()
+vindex(SRANDOM)
+item(tt(SRANDOM)) (
+A random positive 32-bit integer between 0 and 4,294,967,295.  This parameter
+is read-only. The name was chosen for compatibility with Bash and to
+distinguish it from tt(RANDOM) which has a documented repeatable behavior.
+)
+enditem()
+
+subsect(Math Functions)
+
+startitem()
+item(tt(zrand_float+LPAR()RPAR())) (
+Returns a random floating point number between 0 and 1.
+)
+enditem()
+
+startitem()
+item(tt(zrand_int)+LPAR()tt(upper), tt(lower), tt(inclusive)RPAR()) (
+Returns a random integer between tt(lower) and tt(upper). All parameters are
+optional.  If none are specified it is equivalent to
+tt(SRANDOM).
+
+tt(upper) is the upper bound on the resultant number and defaults to
+4,294,967,295.
+
+tt(lower) is the lower bound and defaults to 0.
+
+The defaults of these two arguments are also the maximum and minimum to which
+either can be set.
+
+tt(inclusive) is a flag that controls whether the result is ever equal to
+tt(upper).  By default it is not. If this argument is set to a non-zero value
+then it may be.
+
+This is to facilitate a construct like tt($a[zrand_int($#a)+1]) rather
+than tt($a[zrand_int+LPAR()$#a-1+RPAR()+1]).
+For example, if $#a is 16, you would use tt(zrand_int+LPAR()16RPAR()) which has
+16 possible return values 0-15.  Because the function can return zero, in order
+to use it as an array index from 1-16 you need to add one.  It would
+be an array index range error for it to also potentially return 16 ($#a). You
+could, however, use the construct tt(zrand_int+LPAR()16,1,1+RPAR()) instead of
+adding 1 to achieve the same result, but it is more verbose.
+
+Most statistics algorithms seem to also expect 0 to tt(upper)-1, so this was
+deemed the most commonly desired case and chosen as the default.
+)
+enditem()
+
diff --git a/Src/Modules/random.c b/Src/Modules/random.c
new file mode 100644
index 000000000..b45483afe
--- /dev/null
+++ b/Src/Modules/random.c
@@ -0,0 +1,324 @@
+/*
+ * random.c - module to access kernel random sources.
+ *
+ * This file is part of zsh, the Z shell.
+ *
+ * Copyright (c) 2022 Clinton Bunch
+ * All rights reserved.
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and to distribute modified versions of this software for any
+ * purpose, provided that the above copyright notice and the following
+ * two paragraphs appear in all copies of this software.
+ *
+ * In no event shall Clinton Bunch or the Zsh Development Group be liable
+ * to any party for direct, indirect, special, incidental, or consequential
+ * damages arising out of the use of this software and its documentation,
+ * even if Clinton Bunch and the Zsh Development Group have been advised of
+ * the possibility of such damage.
+ *
+ * Clinton Bunch and the Zsh Development Group specifically disclaim any
+ * warranties, including, but not limited to, the implied warranties of
+ * merchantability and fitness for a particular purpose.  The software
+ * provided hereunder is on an "as is" basis, and Clinton Bunch and the
+ * Zsh Development Group have no obligation to provide maintenance,
+ * support, updates, enhancements, or modifications.
+ *
+ */
+
+#include "random.mdh"
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <stdio.h>
+
+#include <stdbool.h>
+#include <stdint.h>
+
+#ifdef HAVE_SYS_RANDOM_H
+#include <sys/random.h>
+#endif
+
+/* Simplify select URANDOM specific code */
+#if !defined(HAVE_ARC4RANDOM_BUF) && !defined(HAVE_GETRANDOM)
+#define USE_URANDOM
+#endif
+
+/* buffer to pre-load integers for SRANDOM to lessen the context switches */
+uint32_t rand_buff[8];
+static int buf_cnt = -1;
+
+#ifdef USE_URANDOM
+/* File descriptor for /dev/urandom */
+int randfd = -1;
+#endif /* USE_URANDOM */
+
+static zlong get_srandom(UNUSED(Param p));
+
+/**/
+ssize_t
+getrandom_buffer(void *buf, size_t len)
+{
+    ssize_t ret;
+    size_t  val     = 0;
+    uint8_t *bufptr = buf;
+
+    do {
+	errno = 0;
+#ifdef HAVE_ARC4RANDOM_BUF
+	arc4random_buf(buf,len);
+	ret = len;
+#elif defined(HAVE_GETRANDOM)
+	ret=getrandom(bufptr,(len - val),0);
+#else
+	ret=read(randfd,bufptr,(len - val));
+#endif
+	if (ret < 0) {
+	    if (errno != EINTR || errflag || retflag || breaks || contflag) {
+		zwarn("Unable to get random data: %e.", errno);
+		return -1;
+	    }
+	}
+	bufptr += ret;
+	val    += ret;
+    } while (ret < len);
+    return ret;
+}
+
+/*
+ * Generate count number of random 32-bit integers between 0 and max-1 
+ * Got this algorithm from
+ * https://lemire.me/blog/2016/06/30/fast-random-shuffling/
+ * adapting the public domain code.
+ *
+ */
+
+/**/
+void
+get_bound_random_buffer(uint32_t *buffer, size_t count, uint32_t max)
+{
+    uint64_t multi_result;
+    uint32_t threshold;
+    uint32_t leftover;
+
+    size_t i; /* loop counter */
+
+    getrandom_buffer((void*) buffer, count*sizeof(uint32_t));
+    if (max == UINT32_MAX) 
+	return;
+
+    for(i=0;i<count;i++) {
+	multi_result = ((uint64_t) buffer[i]) * (uint64_t) max;
+	leftover = (uint32_t) multi_result;
+
+	/*
+	 * The following if statement should (according to Google's Gemini)
+	 * only be executed with a probability of 1/2**32 or 2.33e-10
+	 */
+	if(leftover < max) {
+	    threshold= -max % max;
+	    while (leftover < threshold) {
+		uint32_t j=get_srandom(NULL);
+		multi_result=(uint64_t) j * (uint64_t) max;
+		leftover= (uint32_t) multi_result;
+	    }
+	}
+	buffer[i]=multi_result >> 32;
+    }
+}
+
+/*
+ * Provides for the SRANDOM parameter and returns an unsigned 32-bit random
+ * integer.
+ */
+
+/**/
+static zlong
+get_srandom(UNUSED(Param pm)) {
+
+    if(buf_cnt <= 0) {
+	getrandom_buffer((void*) rand_buff,sizeof(rand_buff));
+	buf_cnt=sizeof(rand_buff)/sizeof(rand_buff[0]);
+    }
+    return rand_buff[--buf_cnt];
+}
+
+/*
+ * Implements math function zrand_int takes 0 to 3 arguments an upper bound,
+ * a lower bound and a flag as to whether the range is inclusive or not.  The
+ * default is exclusive.  If neither upper or lower is specified this is no
+ * different than SRANDOM.
+ */
+
+/**/
+static mnumber
+math_zrand_int(UNUSED(char *name), int argc, mnumber *argv, UNUSED(int id))
+{
+    mnumber ret;
+    uint32_t i;
+    zlong lower=0, upper=UINT32_MAX,incl=0, diff;
+
+    ret.type = MN_INTEGER;
+
+    switch (argc) {
+	case 0: ret.u.l=get_srandom(NULL);
+		return ret;
+		break;
+	case 3: incl = (argv[2].u.l != 0)?1:0;
+	case 2: lower = argv[1].u.l;
+	case 1: upper = argv[0].u.l;
+	default: diff = upper-lower+incl;
+    }
+
+    if (lower < 0 || lower >= UINT32_MAX) {
+	zwarn("Lower bound (%z) out of range: 0-4294967295",lower);
+    } else if (upper < lower) {
+	zwarn("Upper bound (%z) must be greater than Lower Bound (%z)",upper,lower);
+    } else if (upper < 0 || upper >= UINT32_MAX) {
+	zwarn("Upper bound (%z) out of range: 0-4294967295",upper);
+    }
+
+    if ( diff == 0 ) {
+	ret.u.l=upper; /* still not convinced this shouldn't be an error. */
+    } else {
+	get_bound_random_buffer(&i,1,(uint32_t) diff);
+	ret.u.l=i+lower;
+    }
+    return ret;
+}
+
+/*
+ * Implements the mathfunc zrand_float and returns a random floating-point
+ * number between 0 and 1.  
+ *
+ */
+
+/**/
+static mnumber
+math_zrand_float(UNUSED(char *name), UNUSED(int argc), UNUSED(mnumber *argv),
+	     UNUSED(int id))
+{
+    mnumber ret;
+    double r;
+
+    r = random_real();
+    if (r < 0) {
+	zwarnnam(name, "Failed to get sufficient random data.");
+    }
+    ret.type = MN_FLOAT;
+    ret.u.d = r;
+
+    return ret;
+}
+
+static const struct gsu_integer srandom_gsu =
+{ get_srandom, nullintsetfn, stdunsetfn };
+
+static struct paramdef patab[] = {
+    {"SRANDOM", PM_INTEGER | PM_READONLY_SPECIAL | PM_HIDEVAL, NULL,
+	    &srandom_gsu, NULL, NULL, NULL},
+};
+
+static struct mathfunc mftab[] = {
+    NUMMATHFUNC("zrand_float", math_zrand_float, 0, 0, 0),
+    NUMMATHFUNC("zrand_int", math_zrand_int, 0, 3, 0),
+};
+
+static struct features module_features = {
+    NULL, 0,
+    NULL, 0,
+    mftab, sizeof(mftab)/sizeof(*mftab),
+    patab, sizeof(patab)/sizeof(*patab),
+    0
+};
+
+/**/
+int
+setup_(UNUSED(Module m))
+{
+#ifdef USE_URANDOM
+    /* Check for the existence of /dev/urandom */
+
+    struct stat st;
+
+    errno=0;
+    if (stat("/dev/urandom",&st) < 0) {
+	zwarn("Error getting kernel random pool: %e.", errno);
+	return 1;
+    }
+
+    errno=0;
+    if (!(S_ISCHR(st.st_mode)) ) {
+	zwarn("Error getting kernel random pool: %e.", errno);
+	return 1;
+    }
+#endif /* USE_URANDOM */
+    return 0;
+}
+
+/**/
+int
+features_(Module m, char ***features)
+{
+    *features = featuresarray(m, &module_features);
+    return 0;
+}
+
+/**/
+int
+enables_(Module m, int **enables)
+{
+    return handlefeatures(m, &module_features, enables);
+}
+
+/**/
+int
+boot_(Module m)
+{
+#ifdef USE_URANDOM
+    /*
+     * Open /dev/urandom.  Here because of a weird issue on HP-UX 11.31
+     * When opening in setup_ open returned 0.  In boot_, it returns
+     * an unused file descriptor. Decided against ifdef HPUX as it works
+     * here just as well for other platforms.
+     *
+     */
+
+    int tmpfd=-1;
+
+    errno=0;
+    if ((tmpfd = open("/dev/urandom", O_RDONLY)) < 0) {
+	zwarn("Could not access kernel random pool: %e.",errno);
+	return 1;
+    }
+    randfd = movefd(tmpfd);
+    addmodulefd(randfd,FDT_MODULE);
+    if (randfd < 0) {
+	zwarn("Could not access kernel random pool.");
+	return 1;
+    }
+#endif /* USE_URANDOM */
+    return 0;
+}
+
+/**/
+int
+cleanup_(Module m)
+{
+    return setfeatureenables(m, &module_features, NULL);
+}
+
+/**/
+int
+finish_(UNUSED(Module m))
+{
+#ifdef USE_URANDOM
+    if (randfd >= 0)
+	zclose(randfd);
+#endif /* USE_URANDOM */
+    return 0;
+}
+
diff --git a/Src/Modules/random.mdd b/Src/Modules/random.mdd
new file mode 100644
index 000000000..7a75f29ff
--- /dev/null
+++ b/Src/Modules/random.mdd
@@ -0,0 +1,7 @@
+name=zsh/random
+link=either
+load=yes
+
+autofeatures="p:SRANDOM f:zrand_float f:zrand_int"
+
+objects="random.o random_real.o"
diff --git a/Src/Modules/random_real.c b/Src/Modules/random_real.c
new file mode 100644
index 000000000..900360414
--- /dev/null
+++ b/Src/Modules/random_real.c
@@ -0,0 +1,213 @@
+/*  This file contains code under different copyrights separated by */
+/* ====@@@@@=== */
+
+/*
+ * random_real.c - module to access kernel random sources.
+ *
+ * This file is part of zsh, the Z shell.
+ *
+ * Copyright (c) 2022 Clinton Bunch
+ * All rights reserved.
+ *
+ * Permission is hereby granted, without written agreement and without
+ * license or royalty fees, to use, copy, modify, and distribute this
+ * software and to distribute modified versions of this software for any
+ * purpose, provided that the above copyright notice and the following
+ * two paragraphs appear in all copies of this software.
+ *
+ * In no event shall Clinton Bunch or the Zsh Development Group be liable
+ * to any party for direct, indirect, special, incidental, or consequential
+ * damages arising out of the use of this software and its documentation,
+ * even if Clinton Bunch and the Zsh Development Group have been advised of
+ * the possibility of such damage.
+ *
+ * Clinton Bunch and the Zsh Development Group specifically disclaim any
+ * warranties, including, but not limited to, the implied warranties of
+ * merchantability and fitness for a particular purpose.  The software
+ * provided hereunder is on an "as is" basis, and Clinton Bunch and the
+ * Zsh Development Group have no obligation to provide maintenance,
+ * support, updates, enhancements, or modifications.
+ *
+ */
+
+#include "random.mdh"
+
+#include <math.h>
+#include <stdint.h>
+#include <errno.h>
+
+
+/* Count the number of leading zeros, hopefully in gcc/clang by HW
+ * instruction */
+#if defined(__GNUC__) || defined(__clang__)
+#define clz64(x) __builtin_clzll(x)
+#else
+#define clz64(x) _zclz64(x)
+
+/**/
+int
+_zclz64(uint64_t x) {
+    int n = 0;
+
+    if (x == 0)
+	return 64;
+
+    if (!(x & 0xFFFFFFFF00000000ull)) {
+	n+=32;
+	x<<=32;
+    }
+    if (!(x & 0xFFFF000000000000ull)) {
+	n+=16;
+	x<<=16;
+    }
+    if (!(x & 0xFF00000000000000ull)) {
+	n+=8;
+	x<<=8;
+    }
+    if (!(x & 0xF000000000000000ull)) {
+	n+=4;
+	x<<=4;
+    }
+    if (!(x & 0xC000000000000000ull)) {
+	n+=2;
+	x<<=1;
+    }
+    if (!(x & 0x8000000000000000ull)) {
+	n+=1;
+    }
+    return n;
+}
+#endif /* __GNU_C__ or __clang__ */
+
+/**/
+uint64_t
+random_64bit(void) {
+    uint64_t r;
+
+    if(getrandom_buffer(&r,sizeof(r)) < 0) {
+	zwarn("zsh/random: Can't get sufficient random data.");
+	return 1; /* 0 will cause loop */
+    }
+
+    return r;
+}
+
+/* ====@@@@@=== */
+/* Following code is under the below copyright, despite changes for error
+ * handling and removing GCCisms */
+
+/*-
+ * Copyright (c) 2014 Taylor R. Campbell
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. 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.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR 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 AUTHOR 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.
+ */
+
+/*
+ * Uniform random floats: How to generate a double-precision
+ * floating-point numbers in [0, 1] uniformly at random given a uniform
+ * random source of bits.
+ *
+ * See <http://mumble.net/~campbell/2014/04/28/uniform-random-float>
+ * for explanation.
+ *
+ * Updated 2015-02-22 to replace ldexp(x, <constant>) by x * ldexp(1,
+ * <constant>), since glibc and NetBSD libm both have slow software
+ * bit-twiddling implementations of ldexp, but GCC can constant-fold
+ * the latter.
+ */
+
+/*
+ * random_real: Generate a stream of bits uniformly at random and
+ * interpret it as the fractional part of the binary expansion of a
+ * number in [0, 1], 0.00001010011111010100...; then round it.
+ */
+
+/**/
+double
+random_real(void)
+{
+	int exponent = 0;
+	uint64_t significand = 0;
+	uint64_t r = 0;
+	unsigned shift;
+
+	/*
+	 * Read zeros into the exponent until we hit a one; the rest
+	 * will go into the significand.
+	 */
+	while (significand == 0) {
+		exponent -= 64;
+
+		/* Get random_64bit and check for error */
+		errno = 0;
+		significand = random_64bit();
+		if (errno)
+		    return -1;
+		/*
+		 * If the exponent falls below -1074 = emin + 1 - p,
+		 * the exponent of the smallest subnormal, we are
+		 * guaranteed the result will be rounded to zero.  This
+		 * case is so unlikely it will happen in realistic
+		 * terms only if random_64bit is broken.
+		 */
+		if (exponent < -1074)
+			return 0;
+	}
+
+	/*
+	 * There is a 1 somewhere in significand, not necessarily in
+	 * the most significant position.  If there are leading zeros,
+	 * shift them into the exponent and refill the less-significant
+	 * bits of the significand.  Can't predict one way or another
+	 * whether there are leading zeros: there's a fifty-fifty
+	 * chance, if random_64bit is uniformly distributed.
+	 */
+	shift = clz64(significand);
+	if (shift != 0) {
+
+		errno = 0;
+		r = random_64bit();
+		if (errno)
+		    return -1;
+
+		exponent -= shift;
+		significand <<= shift;
+		significand |= (r >> (64 - shift));
+	}
+
+	/*
+	 * Set the sticky bit, since there is almost surely another 1
+	 * in the bit stream.  Otherwise, we might round what looks
+	 * like a tie to even when, almost surely, were we to look
+	 * further in the bit stream, there would be a 1 breaking the
+	 * tie.
+	 */
+	significand |= 1;
+
+	/*
+	 * Finally, convert to double (rounding) and scale by
+	 * 2^exponent.
+	 */
+	return ldexp((double)significand, exponent);
+}
+/* ====@@@@@=== */
diff --git a/Test/V15random.ztst b/Test/V15random.ztst
new file mode 100644
index 000000000..889324380
--- /dev/null
+++ b/Test/V15random.ztst
@@ -0,0 +1,120 @@
+%prep
+
+  if ! zmodload zsh/mathfunc 2>/dev/null; then
+    ZTST_unimplemented="can't load the zsh/mathfunc module for testing"
+  fi
+
+  print -ru $ZTST_fd 'This test may take two seconds...'
+
+  function calc_chi2() {
+# Calculate Chi Squared
+	integer n=$1
+	for ((k=1;k<=n;k++))
+	do
+		((delta=${observed[$k]}-expected))
+		#echo "$k=$delta \c"
+		((chi2+=delta**2/expected))
+	done
+	#echo
+  }
+  function ligf() {
+	float -E S=$1 Z=$2
+	float -gE RET
+	if (( Z < 0.0 )); then
+		RET=0.0
+	fi
+	float -E Sc=$((1.0/S))
+	((Sc *= Z**S ))
+	((Sc *= exp(-1*Z) ))
+	float -E Sum=1.0 Nom=1.0 Denom=1.0
+	integer I
+	# 200 iterations seems like enough for an approximation of an infinate
+	# series
+	for (( I=0;I<200;I++))
+	do
+		((Nom *= Z))
+		((S++))
+		((Denom *= S))
+		((Sum += (Nom/Denom) ))
+	done
+	RET=$((Sum * Sc))
+  }
+#echo $observed
+  function calc_p() {
+	#echo "chi2=$chi2"
+	float chi2=0 cdf=0
+        float RET
+	calc_chi2 "$n"
+	ligf $(((n-1)*0.5)) $((chi2*0.5))
+	((cdf=$RET/gamma((n-1)*0.5) ))
+	((p=1-cdf))
+  }
+%test
+  fail=0
+  for ((o=0;o<100;o++))
+  do
+      n=20
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( zrand_int(n,1,1)  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:random integer $samples samples between 1-$n
+
+  fail=0
+  unset observed
+  for ((o=0;o<100;o++))
+  do
+      n=20
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( ceil(zrand_float()*n)  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:random float $samples samples
+
+  fail=0
+  unset observed
+  for ((o=0;o<100;o++))
+  do
+      n=16
+      samples=1000
+      integer j=0
+      float -F expected=$((samples*1.0/n))
+      float -F p
+      typeset -a observed=([$n]=0)
+      for ((i=1;i<=samples;i++))
+      do
+        j=$(( (SRANDOM & 0x0F)+1  ))
+        ((observed[$j]++))
+      done
+      calc_p
+      if ((p<0.05)); then
+         ((fail++))
+      fi
+  done
+  ((fail <= 10))
+0q:SRANDOM mod 16 $samples samples
+
diff --git a/configure.ac b/configure.ac
index 78621042d..35f641cde 100644
--- a/configure.ac
+++ b/configure.ac
@@ -636,6 +636,7 @@ fi
 AC_CHECK_HEADERS(sys/time.h sys/times.h sys/select.h termcap.h termio.h \
 		 termios.h sys/param.h sys/filio.h string.h memory.h \
 		 limits.h fcntl.h libc.h sys/utsname.h sys/resource.h \
+                 sys/random.h \
 		 locale.h errno.h stdio.h stdarg.h varargs.h stdlib.h \
 		 unistd.h sys/capability.h \
 		 utmp.h utmpx.h sys/types.h pwd.h grp.h poll.h sys/mman.h \
@@ -1292,6 +1293,7 @@ AC_CHECK_FUNCS(strftime strptime mktime timelocal \
 	       cygwin_conv_path \
 	       nanosleep \
 	       srand_deterministic \
+               getrandom arc4random_buf \
 	       setutxent getutxent endutxent getutent)
 AC_FUNC_STRCOLL