--- a/scripts/statistics/distributions/trnd.m
+++ b/scripts/statistics/distributions/trnd.m
@@ -1,3 +1,4 @@
+## Copyright (C) 2011 Rik Wehbring
 ## Copyright (C) 1995-2011 Kurt Hornik
 ##
 ## This file is part of Octave.
@@ -17,74 +18,100 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {} trnd (@var{n}, @var{r}, @var{c})
-## @deftypefnx {Function File} {} trnd (@var{n}, @var{sz})
-## Return an @var{r} by @var{c} matrix of random samples from the t
-## (Student) distribution with @var{n} degrees of freedom.  @var{n} must
-## be a scalar or of size @var{r} by @var{c}.  Or if @var{sz} is a
-## vector create a matrix of size @var{sz}.
+## @deftypefn  {Function File} {} trnd (@var{n})
+## @deftypefnx {Function File} {} trnd (@var{n}, @var{r})
+## @deftypefnx {Function File} {} trnd (@var{n}, @var{r}, @var{c}, @dots{})
+## @deftypefnx {Function File} {} trnd (@var{n}, [@var{sz}])
+## Return a matrix of random samples from the t (Student) distribution with
+## @var{n} degrees of freedom.
 ##
-## If @var{r} and @var{c} are omitted, the size of the result matrix is
-## the size of @var{n}.
+## When called with a single size argument, return a square matrix with
+## the dimension specified.  When called with more than one scalar argument the
+## first two arguments are taken as the number of rows and columns and any
+## further arguments specify additional matrix dimensions.  The size may also
+## be specified with a vector of dimensions @var{sz}.
+## 
+## If no size arguments are given then the result matrix is the size of
+## @var{n}.
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
 ## Description: Random deviates from the t distribution
 
-function rnd = trnd (n, r, c)
-
-  if (nargin == 3)
-    if (! (isscalar (r) && (r > 0) && (r == round (r))))
-      error ("trnd: R must be a positive integer");
-    endif
-    if (! (isscalar (c) && (c > 0) && (c == round (c))))
-      error ("trnd: C must be a positive integer");
-    endif
-    sz = [r, c];
+function rnd = trnd (n, varargin)
 
-    if (any (size (n) != 1)
-        && ((length (size (n)) != length (sz)) || any (size (n) != sz)))
-      error ("trnd: N must be scalar or of size SZ");
-    endif
-  elseif (nargin == 2)
-    if (isscalar (r) && (r > 0))
-      sz = [r, r];
-    elseif (isvector(r) && all (r > 0))
-      sz = r(:)';
-    else
-      error ("trnd: R must be a positive integer or vector");
-    endif
-
-    if (any (size (n) != 1)
-        && ((length (size (n)) != length (sz)) || any (size (n) != sz)))
-      error ("trnd: N must be scalar or of size SZ");
-    endif
-  elseif (nargin == 1)
-    sz = size (n);
-  else
+  if (nargin < 1)
     print_usage ();
   endif
 
+  if (nargin == 1)
+    sz = size (n);
+  elseif (nargin == 2)
+    if (isscalar (varargin{1}) && varargin{1} >= 0)
+      sz = [varargin{1}, varargin{1}];
+    elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
+      sz = varargin{1};
+    else
+      error ("trnd: dimension vector must be row vector of non-negative integers");
+    endif
+  elseif (nargin > 2)
+    if (any (cellfun (@(x) (!isscalar (x) || x < 0), varargin)))
+      error ("trnd: dimensions must be non-negative integers");
+    endif
+    sz = [varargin{:}];
+  endif
+
+  if (!isscalar (n) && !isequal (size (n), sz))
+    error ("trnd: N must be scalar or of size SZ");
+  endif
+
+  if (iscomplex (n))
+    error ("trnd: N must not be complex");
+  endif
+
+  if (isa (n, "single"))
+    cls = "single";
+  else
+    cls = "double";
+  endif
+
   if (isscalar (n))
-    if (!(n > 0) || !(n < Inf))
-      rnd = NaN (sz);
-    elseif ((n > 0) && (n < Inf))
-      rnd = randn(sz) ./ sqrt(2*randg(n/2,sz)./n);
+    if ((n > 0) && (n < Inf))
+      rnd = randn (sz) ./ sqrt (2*randg (n/2, sz) / n);
     else
-      rnd = zeros (size (n));
+      rnd = NaN (sz, cls);
     endif
   else
-    rnd = zeros (size (n));
+    rnd = NaN (sz, cls);
 
-    k = find (!(n > 0) | !(n < Inf));
-    if (any (k))
-      rnd(k) = NaN;
-    endif
-
-    k = find ((n > 0) & (n < Inf));
-    if (any (k))
-      rnd(k) = randn(size(k)) ./ sqrt(2*randg(n(k)/2,size(k))./n(k));
-    endif
+    k = (n > 0) & (n < Inf);
+    rnd(k) = randn (sum (k(:)), 1) ./ sqrt (2*randg (n(k)/2) ./ n(k))(:);
   endif
 
 endfunction
+
+
+%!assert(size (trnd (2)), [1, 1]);
+%!assert(size (trnd (ones(2,1))), [2, 1]);
+%!assert(size (trnd (ones(2,2))), [2, 2]);
+%!assert(size (trnd (1, 3)), [3, 3]);
+%!assert(size (trnd (1, [4 1])), [4, 1]);
+%!assert(size (trnd (1, 4, 1)), [4, 1]);
+
+%% Test class of input preserved
+%!assert(class (trnd (1)), "double");
+%!assert(class (trnd (single(1))), "single");
+%!assert(class (trnd (single([1 1]))), "single");
+
+%% Test input validation
+%!error trnd ()
+%!error trnd (1, -1)
+%!error trnd (1, ones(2))
+%!error trnd (i)
+%!error trnd (1, [2 -1 2])
+%!error trnd (1, 2, ones(2))
+%!error trnd (1, 2, -1)
+%!error trnd (ones(2,2), 3)
+%!error trnd (ones(2,2), [3, 2])
+%!error trnd (ones(2,2), 2, 3)
+