--- a/NEWS	Tue Dec 18 16:43:05 2018 -0800
+++ b/NEWS	Tue Dec 18 22:04:23 2018 -0800
@@ -8,6 +8,13 @@
     aforementioned functions by properly overloading the "size"
     function.
 
+ ** The function randi has been recoded to produce an unbiased (all
+    results are equally likely) sample of integers.  This may produce
+    different results in existing code.  If it is necessary to reproduce
+    the exact random integer sequence as in previous versions use
+    
+      ri = imin + floor ((imax - imin + 1) * rand ());
+
  ** A new core function movfun will apply a function to a sliding
     window of arbitrary size on a dataset and accumulate the results.
     Many common cases have been implemented using the naming

--- a/scripts/general/randi.m	Tue Dec 18 16:43:05 2018 -0800
+++ b/scripts/general/randi.m	Tue Dec 18 22:04:23 2018 -0800
@@ -79,13 +79,6 @@
     endif
   endif
 
-  if (nargin > 1 && ischar (varargin{end}))
-    rclass = varargin{end};
-    varargin(end) = [];
-  else
-    rclass = "double";
-  endif
-
   ## Limit set by use of class double in rand(): Any consecutive integer in the
   ## range [-flintmax(), flintmax()] can be represented by a double.
   if ((abs (imax) >= flintmax ()) || (abs (imin) >= flintmax ()))
@@ -95,7 +88,35 @@
     error ("randi: integer range must be smaller than flintmax()-1");
   endif
 
-  ri = imin + floor ((imax - imin + 1) * rand (varargin{:}));
+  if (nargin > 1 && ischar (varargin{end}))
+    rclass = varargin{end};
+    varargin(end) = [];
+    nargin = nargin - 1;
+  else
+    rclass = "double";
+  endif
+
+  ## Expand dimension argument to at least 2-D for reshape
+  if (nargin == 1)
+    varargin = {1, 1};
+  elseif (nargin == 2 && isscalar (varargin{1}))
+    varargin(2) = varargin(1);
+  endif
+    
+  ## Rejection Algorithm to guarantee unbiased results.  See bug #54619.
+  rng = (imax - imin) + 1;              # requested range
+  N = prod ([varargin{:}]);             # number of requested elements
+  K = floor ((flintmax () + 1) / rng);  # number of primary integers ...
+                                        # mapped to single output
+  p = (K*rng) / (flintmax () + 1);      # expected proportion of used primaries
+  
+  do  
+    M = ceil (N/p + 10*sqrt (N/p - N)); # number of requested primary integers
+    r_prim = floor (rand (M,1) * (flintmax () + 1));
+    r_prim = r_prim(r_prim < K*rng);
+  until (numel (r_prim) >= N)           # should practically always be true
+  
+  ri = imin + floor (reshape (r_prim(1:N), varargin{:}) / K);
 
   if (! strcmp (rclass, "double"))
     if (strfind (rclass, "int"))