Mercurial > octave

changeset 33437:b6bb53ff12b4 bytecode-interpreter

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Merge default to bytecode-interpreter
author Arun Giridhar <arungiridhar@gmail.com>
date Sun, 21 Apr 2024 11:22:37 -0400
parents 066a27297ba3 (current diff) 60772ed3e920 (diff)
children 2383f7553930
files
diffstat 1 files changed, 20 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/specfun/nchoosek.m	Fri Apr 19 12:57:20 2024 -0400
+++ b/scripts/specfun/nchoosek.m	Sun Apr 21 11:22:37 2024 -0400
@@ -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));