Mercurial > octave

comparison scripts/general/nextpow2.m @ 33199:297fd823a953

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
nextpow2.m: Overhaul function for improved accuracy and Matlab compatibility (bug #65441) * nextpow2.m: Clarify documentation. Document return value of special input 0. Add input validation to exclude complex inputs. New algorithm for computing nextpow2 based on 2-output form of log2(). Use cast() to have output class match input class for Matlab compatibility. Add BIST tests. * NEWS.10.md: Announce changes to nextpow2 function.
author Rik <rik@octave.org>
date Tue, 12 Mar 2024 14:27:13 -0700
parents 2e484f9f1f18
children
comparison
equal deleted inserted replaced
33197:bee755660415 33199:297fd823a953
23 ## 23 ##
24 ######################################################################## 24 ########################################################################
25 25
26 ## -*- texinfo -*- 26 ## -*- texinfo -*-
27 ## @deftypefn {} {@var{n} =} nextpow2 (@var{x}) 27 ## @deftypefn {} {@var{n} =} nextpow2 (@var{x})
28 ## Compute the exponent for the smallest power of two larger than the input. 28 ## Compute the exponent of the next power of two not smaller than the input.
29 ## 29 ##
30 ## For each element in the input array @var{x}, return the first integer 30 ## For each element in the input array @var{x}, return the smallest integer
31 ## @var{n} such that 31 ## @var{n} such that
32 ## @tex 32 ## @tex
33 ## $2^n \ge |x|$. 33 ## $2^n \ge |x|$.
34 ## @end tex 34 ## @end tex
35 ## @ifnottex 35 ## @ifnottex
36 ## 2^n @geq{} abs (x). 36 ## @code{2^@var{n} @geq{} abs (@var{x})}.
37 ## @end ifnottex 37 ## @end ifnottex
38 ## For input elements equal to zero, return zero.
38 ## 39 ##
39 ## @seealso{pow2, log2} 40 ## @seealso{pow2, log2}
40 ## @end deftypefn 41 ## @end deftypefn
41 42
42 function n = nextpow2 (x) 43 function n = nextpow2 (x)
45 print_usage (); 46 print_usage ();
46 endif 47 endif
47 48
48 if (! isnumeric (x)) 49 if (! isnumeric (x))
49 error ("nextpow2: X must be numeric"); 50 error ("nextpow2: X must be numeric");
51 elseif (! isreal (x))
52 error ("nextpow2: X must be real");
50 endif 53 endif
51 54
52 n = ceil (log2 (abs (x))); 55 [f, n] = log2 (abs (x));
53 n(x == 0) = 0; # special case 56 n(f == 0.5)--;
57
58 if (isfloat (x))
59 idx_nan = isnan (x);
60 n(idx_nan) = x(idx_nan);
61 n(isinf (x)) = Inf;
62 else
63 n = cast (n, class (x));
64 if (intmax (x) > flintmax ())
65 n((2 .^ n) < abs (x))++;
66 endif
67 endif
54 68
55 endfunction 69 endfunction
56 70
57 71
58 %!assert (nextpow2 (16), 4) 72 %!assert (nextpow2 (16), 4)
68 %!assert (nextpow2 (0), 0) 82 %!assert (nextpow2 (0), 0)
69 %!assert (nextpow2 (1), 0) 83 %!assert (nextpow2 (1), 0)
70 %!assert (nextpow2 (Inf), Inf) 84 %!assert (nextpow2 (Inf), Inf)
71 %!assert (nextpow2 (-Inf), Inf) 85 %!assert (nextpow2 (-Inf), Inf)
72 %!assert (nextpow2 (NaN), NaN) 86 %!assert (nextpow2 (NaN), NaN)
73 %!assert (nextpow2 ([1, Inf, 3, -Inf, 9, NaN]), [0, Inf, 2, Inf, 4, NaN]) 87 %!assert (nextpow2 (NA), NA)
88 %!assert (nextpow2 ([1, Inf, 3, -Inf, 9, NaN, NA]), ...
89 %! [0, Inf, 2, Inf, 4, NaN, NA])
90
91 %!test
92 %! p = (-1074:1023).';
93 %! x = 2 .^ p;
94 %! x = [x, x + eps(x)];
95 %! x = [x, -x];
96 %! n = nextpow2 (x);
97 %! assert (n(:, 1), p);
98 %! assert (n(:, 3), p);
99 %! assert (n(:, 2), p + 1);
100 %! assert (n(:, 4), p + 1);
101
102 %!assert (nextpow2 (realmax ()), 1024)
103 %!assert (nextpow2 (-realmax ()), 1024)
104
105 %!test
106 %! p = single (-149:127).';
107 %! x = 2 .^ p;
108 %! x = [x, x + eps(x)];
109 %! x = [x, -x];
110 %! n = nextpow2 (x);
111 %! assert (n(:, 1), p);
112 %! assert (n(:, 3), p);
113 %! assert (n(:, 2), p + 1);
114 %! assert (n(:, 4), p + 1);
115
116 %!assert (nextpow2 (realmax ('single')), single (128))
117 %!assert (nextpow2 (-realmax ('single')), single (128))
118
119 %!test
120 %! p = int32 (0:30).';
121 %! x = 2 .^ p;
122 %! x = [x, x + 1];
123 %! x = [x, -x];
124 %! n = nextpow2 (x);
125 %! assert (n(:, 1), p);
126 %! assert (n(:, 3), p);
127 %! assert (n(:, 2), p + 1);
128 %! assert (n(:, 4), p + 1);
129
130 %!assert (nextpow2 (int32 (0)), int32 (0))
131 %!assert (nextpow2 (intmin ('int32')), int32 (31))
132 %!assert (nextpow2 (intmax ('int32')), int32 (31))
133
134 %!assert (nextpow2 (uint32 (0)), uint32 (0))
135 %!assert (nextpow2 (intmax ('uint32')), uint32 (32))
136
137 %!test
138 %! p = int64 (0:62).';
139 %! x = 2 .^ p;
140 %! x = [x, x + 1];
141 %! x = [x, -x];
142 %! n = nextpow2 (x);
143 %! assert (n(:, 1), p);
144 %! assert (n(:, 3), p);
145 %! assert (n(:, 2), p + 1);
146 %! assert (n(:, 4), p + 1);
147
148 %!assert (nextpow2 (int64 (0)), int64 (0))
149 %!assert (nextpow2 (intmin ('int64')), int64 (63))
150 %!assert (nextpow2 (intmax ('int64')), int64 (63))
151
152 %!assert (nextpow2 (uint64 (0)), uint64 (0))
153 %!assert (nextpow2 (intmax ('uint64')), uint64 (64))
74 154
75 ## Test input validation 155 ## Test input validation
76 %!error <Invalid call> nextpow2 () 156 %!error <Invalid call> nextpow2 ()
77 %!error <X must be numeric> nextpow2 ("t") 157 %!error <X must be numeric> nextpow2 ("t")
158 %!error <X must be real> nextpow2 (1 + 2i)