Mercurial > octave

changeset 33436:60772ed3e920

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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));