--- a/etc/NEWS.10.md	Tue Mar 12 15:05:37 2024 -0400
+++ b/etc/NEWS.10.md	Tue Mar 12 14:27:13 2024 -0700
@@ -17,6 +17,11 @@
 
 - `nchoosek` algorithm is now ~2x faster and provides greater precision.
 
+- `nextpow2` algorithm is now more accurate for inputs very close to a power
+of 2.  The output class now matches the input class for Matlab compatibility.
+The function no longer accepts complex inputs and will emit an error for these
+inputs.
+
 ### Graphical User Interface
 
 ### Graphics backend

--- a/scripts/general/nextpow2.m	Tue Mar 12 15:05:37 2024 -0400
+++ b/scripts/general/nextpow2.m	Tue Mar 12 14:27:13 2024 -0700
@@ -25,16 +25,17 @@
 
 ## -*- texinfo -*-
 ## @deftypefn {} {@var{n} =} nextpow2 (@var{x})
-## Compute the exponent for the smallest power of two larger than the input.
+## Compute the exponent of the next power of two not smaller than the input.
 ##
-## For each element in the input array @var{x}, return the first integer
+## For each element in the input array @var{x}, return the smallest integer
 ## @var{n} such that
 ## @tex
 ## $2^n \ge |x|$.
 ## @end tex
 ## @ifnottex
-## 2^n @geq{} abs (x).
+## @code{2^@var{n} @geq{} abs (@var{x})}.
 ## @end ifnottex
+## For input elements equal to zero, return zero.
 ##
 ## @seealso{pow2, log2}
 ## @end deftypefn
@@ -47,10 +48,23 @@
 
   if (! isnumeric (x))
     error ("nextpow2: X must be numeric");
+  elseif (! isreal (x))
+    error ("nextpow2: X must be real");
   endif
 
-  n = ceil (log2 (abs (x)));
-  n(x == 0) = 0;  # special case
+  [f, n] = log2 (abs (x));
+  n(f == 0.5)--;
+
+  if (isfloat (x))
+    idx_nan = isnan (x);
+    n(idx_nan) = x(idx_nan);
+    n(isinf (x)) = Inf;
+  else
+    n = cast (n, class (x));
+    if (intmax (x) > flintmax ())
+      n((2 .^ n) < abs (x))++;
+    endif
+  endif
 
 endfunction
 
@@ -70,8 +84,75 @@
 %!assert (nextpow2 (Inf), Inf)
 %!assert (nextpow2 (-Inf), Inf)
 %!assert (nextpow2 (NaN), NaN)
-%!assert (nextpow2 ([1, Inf, 3, -Inf, 9, NaN]), [0, Inf, 2, Inf, 4, NaN])
+%!assert (nextpow2 (NA), NA)
+%!assert (nextpow2 ([1, Inf, 3, -Inf, 9, NaN, NA]), ...
+%!                  [0, Inf, 2,  Inf, 4, NaN, NA])
+
+%!test
+%! p = (-1074:1023).';
+%! x = 2 .^ p;
+%! x = [x, x + eps(x)];
+%! x = [x, -x];
+%! n = nextpow2 (x);
+%! assert (n(:, 1), p);
+%! assert (n(:, 3), p);
+%! assert (n(:, 2), p + 1);
+%! assert (n(:, 4), p + 1);
+
+%!assert (nextpow2 (realmax ()), 1024)
+%!assert (nextpow2 (-realmax ()), 1024)
+
+%!test
+%! p = single (-149:127).';
+%! x = 2 .^ p;
+%! x = [x, x + eps(x)];
+%! x = [x, -x];
+%! n = nextpow2 (x);
+%! assert (n(:, 1), p);
+%! assert (n(:, 3), p);
+%! assert (n(:, 2), p + 1);
+%! assert (n(:, 4), p + 1);
+
+%!assert (nextpow2 (realmax ('single')), single (128))
+%!assert (nextpow2 (-realmax ('single')), single (128))
+
+%!test
+%! p = int32 (0:30).';
+%! x = 2 .^ p;
+%! x = [x, x + 1];
+%! x = [x, -x];
+%! n = nextpow2 (x);
+%! assert (n(:, 1), p);
+%! assert (n(:, 3), p);
+%! assert (n(:, 2), p + 1);
+%! assert (n(:, 4), p + 1);
+
+%!assert (nextpow2 (int32 (0)), int32 (0))
+%!assert (nextpow2 (intmin ('int32')), int32 (31))
+%!assert (nextpow2 (intmax ('int32')), int32 (31))
+
+%!assert (nextpow2 (uint32 (0)), uint32 (0))
+%!assert (nextpow2 (intmax ('uint32')), uint32 (32))
+
+%!test
+%! p = int64 (0:62).';
+%! x = 2 .^ p;
+%! x = [x, x + 1];
+%! x = [x, -x];
+%! n = nextpow2 (x);
+%! assert (n(:, 1), p);
+%! assert (n(:, 3), p);
+%! assert (n(:, 2), p + 1);
+%! assert (n(:, 4), p + 1);
+
+%!assert (nextpow2 (int64 (0)), int64 (0))
+%!assert (nextpow2 (intmin ('int64')), int64 (63))
+%!assert (nextpow2 (intmax ('int64')), int64 (63))
+
+%!assert (nextpow2 (uint64 (0)), uint64 (0))
+%!assert (nextpow2 (intmax ('uint64')), uint64 (64))
 
 ## Test input validation
 %!error <Invalid call> nextpow2 ()
 %!error <X must be numeric> nextpow2 ("t")
+%!error <X must be real> nextpow2 (1 + 2i)