Mercurial > octave
changeset 31878:942ca5f2c71c stable
mean: Process single as double to avoid precision limits on stable (bug #63848).
* mean.m: Change all sum calls to process as "double" to avoid single
precision limits from producing erroneous output. Update output
conversion block to preserve single output for 'default' after forcing
double processing and replace expensive call to cast with feval in
'native'. Add BISTs to test calls to sum that could lose precision. Add
xfail BIST showing the same problem with large doubles.
* NEWS.8.md: Update description of changes to mean.m.
| author | Nicholas R. Jankowski <jankowski.nicholas@gmail.com> |
|---|---|
| date | Thu, 02 Mar 2023 21:42:42 -0500 |
| parents | 715c212605b7 |
| children | fc40c78a8421 b8205aa00a05 |
| files | etc/NEWS.8.md scripts/statistics/mean.m |
| diffstat | 2 files changed, 24 insertions(+), 12 deletions(-) [+] |
line wrap: on
line diff
--- a/etc/NEWS.8.md Sun Feb 26 19:54:58 2023 -0800 +++ b/etc/NEWS.8.md Thu Mar 02 21:42:42 2023 -0500 @@ -39,6 +39,9 @@ - Octave is now compatible with PCRE2 (UTF-8). PCRE2 is preferred over PCRE if both are installed. Configure with `--without-pcre2` if you prefer Octave to use PCRE in this case. + +- `mean` now internally processes data as type double to reduce liklihood of +hitting overflow or precision limits with other types (bug #63848). ### Graphical User Interface
--- a/scripts/statistics/mean.m Sun Feb 26 19:54:58 2023 -0800 +++ b/scripts/statistics/mean.m Thu Mar 02 21:42:42 2023 -0500 @@ -139,7 +139,7 @@ n -= sum (idx(:)); x(idx) = 0; endif - y = sum (x(:), 1) ./ n; + y = sum (x(:), 1, "double") ./ n; else sz = size (x); ## Find the first non-singleton dimension. @@ -150,7 +150,7 @@ n = sum (! idx, dim); x(idx) = 0; endif - y = sum (x, dim) ./ n; + y = sum (x, dim, "double") ./ n; endif else @@ -169,7 +169,7 @@ n = sum (! idx, dim); x(idx) = 0; endif - y = sum (x, dim) ./ n; + y = sum (x, dim, "double") ./ n; else @@ -189,7 +189,7 @@ n -= sum (idx(:)); x(idx) = 0; endif - y = sum (x(:), 1) ./ n; + y = sum (x(:), 1, "double") ./ n; ## for 1 dimension left, return column vector case 1 @@ -200,7 +200,7 @@ if (omitnan) x_vec = x_vec(! isnan (x_vec)); endif - y(i) = sum (x_vec, 1) ./ numel (x_vec); + y(i) = sum (x_vec, 1, "double") ./ numel (x_vec); endfor y = ipermute (y, [misdim, dim]); @@ -214,7 +214,7 @@ if (omitnan) x_vec = x_vec(! isnan (x_vec)); endif - y(i,j) = sum (x_vec, 1) ./ numel (x_vec); + y(i,j) = sum (x_vec, 1, "double") ./ numel (x_vec); endfor endfor y = ipermute (y, [misdim, dim]); @@ -230,17 +230,15 @@ ## Convert output if requested switch (outtype) case "default" - ## do nothing, the operators already do the right thing. + if isa(x, "single") + y = single (y); + endif case "double" y = double (y); case "native" if (! islogical (x)) - y = cast (y, class (x)); + y = feval (class (x), y); endif - otherwise - ## FIXME: This is unreachable code. Valid values already - ## checked in input validation. - error ("mean: OUTTYPE '%s' not recognized", outtype); endswitch endfunction @@ -338,6 +336,17 @@ %! m(2,3) = 15.52301255230125; %! assert (mean (x, [3 2], "omitnan"), m, 4e-14); +## Test limits of single precision summation limits on each code path +%!assert <*63848> (mean (ones (80e6, 1, "single")), 1, eps) +%!assert <*63848> (mean (ones (80e6, 1, "single"), "all"), 1, eps) +%!assert <*63848> (mean (ones (80e6, 1, "single"), 1), 1, eps) +%!assert <*63848> (mean (ones (80e6, 1, "single"), [1 2]), 1, eps) +%!assert <*63848> (mean (ones (80e6, 1, "single"), [1 3]), 1, eps) + +## Test limits of double precision summation +%!assert <63848> (mean ([flintmax("double"), ones(1, 2^8-1, "double")]), ... +%! 35184372088833-1/(2^8), eps(35184372088833)) + ## Test input validation %!error <Invalid call> mean () %!error <Invalid call> mean (1, 2, 3)