# HG changeset patch
# User Arun Giridhar <arungiridhar@gmail.com>
# Date 1713712957 14400
# Node ID b6bb53ff12b41fa4bbce520b2c43cfc3e3a1ad79
# Parent  066a27297ba37172ea9466903b0fc75bb0c7c2fa# Parent  60772ed3e9203768f759734c3c60a49afb3cbd68
maint: Merge default to bytecode-interpreter

diff -r 066a27297ba3 -r b6bb53ff12b4 scripts/specfun/nchoosek.m
--- 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));