--- a/scripts/specfun/bessel.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/bessel.m	Sat Sep 20 06:31:27 2014 -0700
@@ -31,24 +31,21 @@
 ##
 ## @item bessely
 ## Bessel functions of the second kind.  If the argument @var{opt} is supplied,
-## the result is multiplied by @code{exp (-abs (imag (x)))}.
+## the result is multiplied by @w{@code{exp (-abs (imag (x)))}}.
 ##
 ## @item besseli
 ## Modified Bessel functions of the first kind.  If the argument @var{opt} is
-## supplied,
-## the result is multiplied by @code{exp (-abs (real (x)))}.
+## supplied, the result is multiplied by @w{@code{exp (-abs (real (x)))}}.
 ##
 ## @item besselk
 ## Modified Bessel functions of the second kind.  If the argument @var{opt} is
-## supplied,
-## the result is multiplied by @code{exp (x)}.
+## supplied, the result is multiplied by @w{@code{exp (x)}}.
 ##
 ## @item besselh
-## Compute Hankel functions of the first (@var{k} = 1) or second (@var{k}
-## = 2) kind.  If the argument @var{opt} is supplied, the result is multiplied
-## by
-## @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
-## @var{k} = 2.
+## Compute Hankel functions of the first (@var{k} = 1) or second (@var{k} = 2)
+## kind.  If the argument @var{opt} is supplied, the result is multiplied by
+## @w{@code{exp (-I*@var{x})}} for @var{k} = 1 or @w{@code{exp (I*@var{x})}}
+## for @var{k} = 2.
 ## @end table
 ##
 ## If @var{alpha} is a scalar, the result is the same size as @var{x}.
@@ -58,11 +55,10 @@
 ## @code{length (@var{alpha})} columns.  Otherwise, @var{alpha} and
 ## @var{x} must conform and the result will be the same size.
 ##
-## The value of @var{alpha} must be real.  The value of @var{x} may be
-## complex.
+## The value of @var{alpha} must be real.  The value of @var{x} may be complex.
 ##
-## If requested, @var{ierr} contains the following status information
-## and is the same size as the result.
+## If requested, @var{ierr} contains the following status information and is
+## the same size as the result.
 ##
 ## @enumerate 0
 ## @item
@@ -75,8 +71,8 @@
 ## Overflow, return @code{Inf}.
 ##
 ## @item
-## Loss of significance by argument reduction results in less than
-## half of machine accuracy.
+## Loss of significance by argument reduction results in less than half of
+## machine accuracy.
 ##
 ## @item
 ## Complete loss of significance by argument reduction, return @code{NaN}.
@@ -88,7 +84,7 @@
 ## @end deftypefn
 
 function bessel ()
-  error ("bessel: you must use besselj, bessely, besseli, or besselk");
+  error ("bessel: you must use besselj, bessely, besseli, besselk, or besselh");
 endfunction
 
 

--- a/scripts/specfun/beta.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/beta.m	Sat Sep 20 06:31:27 2014 -0700
@@ -18,7 +18,9 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Mapping Function} {} beta (@var{a}, @var{b})
-## For real inputs, return the Beta function,
+## Compute the Beta function for real inputs @var{a} and @var{b}. 
+##
+## The Beta function definition is
 ## @tex
 ## $$
 ##  B (a, b) = {\Gamma (a) \Gamma (b) \over \Gamma (a + b)}.
@@ -31,7 +33,12 @@
 ## @end example
 ##
 ## @end ifnottex
-## @seealso{betaln, betainc}
+##
+## The Beta function can grow quite large and it is often more useful to work
+## with the logarithm of the output rather than the function directly.
+## @xref{XREFbetaln,,betaln}, for computing the logarithm of the Beta function
+## in an efficient manner.
+## @seealso{betaln, betainc, betaincinv}
 ## @end deftypefn
 
 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
@@ -44,12 +51,10 @@
     print_usage ();
   endif
 
-  if (any (size (a) != size (b)) && numel (a) != 1 && numel (b) != 1)
-    error ("beta: inputs A and B have inconsistent sizes");
-  endif
-
   if (! isreal (a) || ! isreal (b))
-    error ("beta: inputs A and B must be real");
+    error ("beta: A and B must be real");
+  elseif (! size_equal (a, b) && numel (a) != 1 && numel (b) != 1)
+    error ("beta: A and B must have consistent sizes");
   endif
 
   retval = real (exp (gammaln (a) + gammaln (b) - gammaln (a+b)));
@@ -60,8 +65,8 @@
 %!test
 %! a = [1, 1.5, 2, 3];
 %! b = [4, 3, 2, 1];
-%! v1 = beta (a,b);
-%! v2 = beta (b,a);
+%! v1 = beta (a, b);
+%! v2 = beta (b, a);
 %! v3 = gamma (a).*gamma (b) ./ gamma (a+b);
 %! assert (v1, v2, sqrt (eps));
 %! assert (v2, v3, sqrt (eps));
@@ -82,4 +87,10 @@
 
 %!error beta ()
 %!error beta (1)
+%!error beta (1,2,3)
+%!error <A and B must be real> beta (1i, 2)
+%!error <A and B must be real> beta (2, 1i)
+%!error <A and B must have consistent sizes> beta ([1 2], [1 2 3])
+%!error <A and B must have consistent sizes> beta ([1 2 3], [1 2])
+%!error <A and B must have consistent sizes> beta ([1 2 3], [1 2 3]')
 

--- a/scripts/specfun/nthroot.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/nthroot.m	Sat Sep 20 06:31:27 2014 -0700
@@ -23,8 +23,13 @@
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} nthroot (@var{x}, @var{n})
 ##
-## Compute the n-th root of @var{x}, returning real results for real
-## components of @var{x}.  For example:
+## Compute the real (non-complex) @var{n}-th root of @var{x}.
+##
+## @var{x} must have all real entries and @var{n} must be a scalar.
+## If @var{n} is an even integer and @var{x} has negative entries then
+## @code{nthroot} aborts and issues an error.
+##
+## Example:
 ##
 ## @example
 ## @group
@@ -34,10 +39,6 @@
 ## @result{} 0.50000 - 0.86603i
 ## @end group
 ## @end example
-##
-## @var{x} must have all real entries.  @var{n} must be a scalar.
-## If @var{n} is an even integer and @var{X} has negative entries, an
-## error is produced.
 ## @seealso{realsqrt, sqrt, cbrt}
 ## @end deftypefn
 
@@ -47,12 +48,12 @@
     print_usage ();
   endif
 
-  if (any (iscomplex (x(:))))
+  if (iscomplex (x))
     error ("nthroot: X must not contain complex values");
   endif
 
-  if (! isscalar (n) || n == 0)
-    error ("nthroot: N must be a nonzero scalar");
+  if (! isreal (n) || ! isscalar (n) || n == 0)
+    error ("nthroot: N must be a real nonzero scalar");
   endif
 
   if (n == 3)
@@ -63,16 +64,17 @@
     y = 1 ./ nthroot (x, -n);
   else
     ## Compute using power.
-    if (n == fix (n) && mod (n, 2) == 1)
+    integer_n = n == fix (n);
+    if (integer_n && mod (n, 2) == 1)
       y = abs (x) .^ (1/n) .* sign (x);
     elseif (any (x(:) < 0))
-      error ("nthroot: if X contains negative values, N must be an odd integer");
+      error ("nthroot: N must be an odd integer if X contains negative values");
     else
       y = x .^ (1/n);
     endif
 
-    if (finite (n) && n > 0 && n == fix (n))
-      ## Correction.
+    if (integer_n && n > 0 && finite (n))
+      ## FIXME: What is this correction for?
       y = ((n-1)*y + x ./ (y.^(n-1))) / n;
       y = merge (finite (y), y, x);
     endif
@@ -82,18 +84,19 @@
 endfunction
 
 
-%!assert (nthroot (-32,5), -2);
-%!assert (nthroot (81,4), 3);
-%!assert (nthroot (Inf,4), Inf);
-%!assert (nthroot (-Inf,7), -Inf);
-%!assert (nthroot (-Inf,-7), 0);
+%!assert (nthroot (-32, 5), -2)
+%!assert (nthroot (81, 4), 3)
+%!assert (nthroot (Inf, 4), Inf)
+%!assert (nthroot (-Inf, 7), -Inf)
+%!assert (nthroot (-Inf, -7), 0)
 
 %% Test input validation
-%!error (nthroot ())
-%!error (nthroot (1))
-%!error (nthroot (1,2,3))
-%!error (nthroot (1+j,2))
-%!error (nthroot (1,[1 2]))
-%!error (nthroot (1,0))
-%!error (nthroot (-1,2))
+%!error nthroot ()
+%!error nthroot (1)
+%!error nthroot (1,2,3)
+%!error <X must not contain complex values> nthroot (1+j, 2)
+%!error <N must be a real nonzero scalar> nthroot (1, i)
+%!error <N must be a real nonzero scalar> nthroot (1, [1 2])
+%!error <N must be a real nonzero scalar> nthroot (1, 0)
+%!error <N must be an odd integer> nthroot (-1, 2)
 

--- a/scripts/specfun/perms.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/perms.m	Sat Sep 20 06:31:27 2014 -0700
@@ -19,15 +19,17 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} perms (@var{v})
+## Generate all permutations of @var{v} with one row per permutation.
 ##
-## Generate all permutations of @var{v}, one row per permutation.  The
-## result has size @code{factorial (@var{n}) * @var{n}}, where @var{n}
+## The result has size @code{factorial (@var{n}) * @var{n}}, where @var{n}
 ## is the length of @var{v}.
 ##
-## As an example, @code{perms ([1, 2, 3])} returns the matrix
+## Example
 ##
 ## @example
 ## @group
+## perms ([1, 2, 3])
+## @result{}
 ##   1   2   3
 ##   2   1   3
 ##   1   3   2
@@ -36,12 +38,18 @@
 ##   3   2   1
 ## @end group
 ## @end example
+##
+## Programming Note: The maximum length of @var{v} should be less than or
+## equal to 10 to limit memory consumption.
+## @seealso{permute, randperm, nchoosek}
 ## @end deftypefn
 
 function A = perms (v)
+
   if (nargin != 1)
     print_usage ();
   endif
+
   vidx = uint8 ([1:length(v)]');
   n = length (vidx);
 
@@ -62,7 +70,9 @@
       endfor
     endfor
   endif
+
   A = v(p);
+
 endfunction
 
 

--- a/scripts/specfun/pow2.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/pow2.m	Sat Sep 20 06:31:27 2014 -0700
@@ -17,9 +17,9 @@
 ## <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Mapping Function} {} pow2 (@var{x})
-## @deftypefnx {Mapping Function} {} pow2 (@var{f}, @var{e})
-## With one argument, computes
+## @deftypefn  {Function File} {} pow2 (@var{x})
+## @deftypefnx {Function File} {} pow2 (@var{f}, @var{e})
+## With one input argument, compute
 ## @tex
 ## $2^x$
 ## @end tex
@@ -28,14 +28,14 @@
 ## @end ifnottex
 ## for each element of @var{x}.
 ##
-## With two arguments, returns
+## With two input arguments, return
 ## @tex
 ## $f \cdot 2^e$.
 ## @end tex
 ## @ifnottex
 ## f .* (2 .^ e).
 ## @end ifnottex
-## @seealso{log2, nextpow2}
+## @seealso{log2, nextpow2, power}
 ## @end deftypefn
 
 ## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
@@ -67,4 +67,5 @@
 %! assert (pow2 (x,y), z, sqrt (eps));
 
 %!error pow2 ()
+%!error pow2 (1,2,3)
 

--- a/scripts/specfun/reallog.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/reallog.m	Sat Sep 20 06:31:27 2014 -0700
@@ -18,12 +18,15 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} reallog (@var{x})
-## Return the real-valued natural logarithm of each element of @var{x}.  Report
-## an error if any element results in a complex return value.
+## Return the real-valued natural logarithm of each element of @var{x}.
+##
+## If any element results in a complex return value @code{reallog} aborts
+## and issues an error.
 ## @seealso{log, realpow, realsqrt}
 ## @end deftypefn
 
 function y = reallog (x)
+
   if (nargin != 1)
     print_usage ();
   elseif (iscomplex (x) || any (x(:) < 0))
@@ -31,13 +34,17 @@
   else
     y = log (x);
   endif
+
 endfunction
 
 
-%!assert (log (1:5), reallog (1:5))
+%!assert (reallog (1:5), log (1:5))
 %!test
-%! x = rand (10,10);
-%! assert (log (x),reallog (x));
+%! x = rand (10, 10);
+%! assert (reallog (x), log (x));
 
+%!error reallog ()
+%!error reallog (1,2)
+%!error <produced complex result> reallog (2i)
 %!error <produced complex result> reallog (-1)
 

--- a/scripts/specfun/realpow.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/realpow.m	Sat Sep 20 06:31:27 2014 -0700
@@ -18,31 +18,37 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} realpow (@var{x}, @var{y})
-## Compute the real-valued, element-by-element power operator.  This is
-## equivalent to @w{@code{@var{x} .^ @var{y}}}, except that @code{realpow}
-## reports an error if any return value is complex.
-## @seealso{reallog, realsqrt}
+## Compute the real-valued, element-by-element power operator.
+##
+## This is equivalent to @w{@code{@var{x} .^ @var{y}}}, except that
+## @code{realpow} reports an error if any return value is complex.
+## @seealso{power, reallog, realsqrt}
 ## @end deftypefn
 
 function z = realpow (x, y)
+
   if (nargin != 2)
     print_usage ();
-  else
-    z = x .^ y;
-    if (iscomplex (z))
-      error ("realpow: produced complex result");
-    endif
   endif
+
+  z = x .^ y;
+  if (iscomplex (z))
+    error ("realpow: produced complex result");
+  endif
+
 endfunction
 
 
-%!assert (power (1:10, 0.5:0.5:5), realpow (1:10, 0.5:0.5:5))
-%!assert ([1:10] .^ [0.5:0.5:5], realpow (1:10, 0.5:0.5:5))
+%!assert (realpow (1:10, 0.5:0.5:5), power (1:10, 0.5:0.5:5))
+%!assert (realpow (1:10, 0.5:0.5:5), [1:10] .^ [0.5:0.5:5])
 %!test
-%! x = rand (10,10);
-%! y = randn (10,10);
+%! x = rand (10, 10);
+%! y = randn (10, 10);
 %! assert (x.^y, realpow (x,y));
 %!assert (realpow (1i,2), -1)
 
+%!error realpow ()
+%!error realpow (1)
+%!error realpow (1,2,3)
 %!error <produced complex result> realpow (-1, 1/2)
 

--- a/scripts/specfun/realsqrt.m	Sat Sep 20 06:21:29 2014 -0700
+++ b/scripts/specfun/realsqrt.m	Sat Sep 20 06:31:27 2014 -0700
@@ -18,12 +18,15 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {Function File} {} realsqrt (@var{x})
-## Return the real-valued square root of each element of @var{x}.  Report an
-## error if any element results in a complex return value.
+## Return the real-valued square root of each element of @var{x}.
+##
+## If any element results in a complex return value @code{realsqrt} aborts
+## and issues an error.
 ## @seealso{sqrt, realpow, reallog}
 ## @end deftypefn
 
 function y = realsqrt (x)
+
   if (nargin != 1)
     print_usage ();
   elseif (iscomplex (x) || any (x(:) < 0))
@@ -31,13 +34,16 @@
   else
     y = sqrt (x);
   endif
+
 endfunction
 
 
-%!assert (sqrt (1:5), realsqrt (1:5))
+%!assert (realsqrt (1:5), sqrt (1:5))
 %!test
-%! x = rand (10,10);
-%! assert (sqrt (x), realsqrt (x));
+%! x = rand (10, 10);
+%! assert (realsqrt (x), sqrt (x));
 
+%!error realsqrt ()
+%!error realsqrt (1,2)
 %!error <produced complex result> realsqrt (-1)