Mercurial > octave
changeset 33436:60772ed3e920
nchoosek: Avoid overflow sacrificing precision if necessary (bug #65495).
* scripts/specfun/nchoosek.m: Try to avoid integer saturation or floating point
overflow by reordering operations. Try to keep precision otherwise. Add
self-test involving high floating point numbers.
author | Hendrik Koerner <koerhen@web.de> |
---|---|
date | Wed, 20 Mar 2024 07:35:25 +0100 |
parents | 69eb4c27d8c8 |
children | b6bb53ff12b4 b2333101ca29 |
files | scripts/specfun/nchoosek.m |
diffstat | 1 files changed, 20 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/nchoosek.m Fri Apr 19 18:30:39 2024 +0200 +++ b/scripts/specfun/nchoosek.m Wed Mar 20 07:35:25 2024 +0100 @@ -140,11 +140,23 @@ g1 = gcd (C, i); g2 = gcd (v - k + i, i/g1); C /= g1; - C *= (v - k + i)/g2; - if (is_int && (C == imax)) + + ## In theory and (always for integers) i/(g1 * g2) is identical to 1 by + ## design. Or for floats and beyond flintmax, the gcd may not be + ## correctly derived by the gcd function and i/(g1 * g2) may not be 1. + C_next = C * ((v - k + i)/g2); + if (is_int || (i/(g1 * g2) == 1) || ! isinf (C_next)) + C = C_next; + C /= i/(g1 * g2); + else + C /= i/(g1 * g2); + ## We have potential precision loss by dividing (too) early, but + ## advantage is that we prevent possible interim overflows + C *= (v - k + i)/g2; + endif + if (is_int && (C == imax)) || (! is_int && isinf (C)) break; # Stop here; saturation reached. endif - C /= i/(g1 * g2); else C *= (v - k + i); C /= i; @@ -288,6 +300,11 @@ %! assert (x, uint8 (252)); %! assert (class (x), "uint8"); +## Floating point number above flintmax +%!test <*65495> +%! warning ("off", "Octave:nchoosek:large-output-float", "local"); +%! assert (! isinf (nchoosek (1024, 512))) + ## Test combining rules for integers and floating point %!test %! x = nchoosek (uint8 (10), single (5));