master 0e00f19 2/2: Calc: fix binomial coefficients for negative arguments (bug#16999)

[Mattias Engdegård]

branch: master
commit 0e00f199cd3599977f75326bb7adc9d70390661e
Author: Mattias Engdegård <mattiase@acm.org>
Commit: Mattias Engdegård <mattiase@acm.org>

    Calc: fix binomial coefficients for negative arguments (bug#16999)
    
    For some values outside integers 0≤k≤n, (n choose k) gave wrong
    results, entered infinite recursion or used unreasonably amounts of
    stack space.  This change fixes that and extends the function to all
    integer arguments using the definitions from M. J. Kronenburg
    (https://arxiv.org/abs/1105.3689).
    
    * lisp/calc/calc-comb.el (calcFunc-choose):
    Fix sign error to prevent infinite recursion and extend function to
    handle all integer arguments.
    (math-choose-iter, math-choose-float-iter): Rewrite in iterative form;
    no TCO in elisp yet.
    * test/lisp/calc/calc-tests.el (calc-tests--fac, calc-tests--choose)
    (calc-tests--check-choose, calc-tests--explain-choose)
    (calc-tests--calc-to-number): New helper functions.
    (calc-choose): New test.
---
 lisp/calc/calc-comb.el       | 42 ++++++++++++++++++---------
 test/lisp/calc/calc-tests.el | 67 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 96 insertions(+), 13 deletions(-)

diff --git a/lisp/calc/calc-comb.el b/lisp/calc/calc-comb.el
index 2efeb7f..f7e29c6 100644
--- a/lisp/calc/calc-comb.el
+++ b/lisp/calc/calc-comb.el
@@ -439,12 +439,25 @@
           (math-div (calcFunc-fact (math-float n))
                     (math-mul (calcFunc-fact m)
                               (calcFunc-fact (math-sub n m))))))
-       ((math-negp m) 0)
-       ((math-negp n)
-        (let ((val (calcFunc-choose (math-add (math-add n m) -1) m)))
+        ;; For the extension to negative integer arguments we follow
+        ;; M. J. Kronenburg, The Binomial Coefficient for Negative Arguments,
+        ;; arXiv:1105.3689v2
+        ((and (math-negp n) (not (math-negp m)))
+         ;; n<0≤m: (n choose m) = (-1)^m (-n+m-1 choose m)
+        (let ((val (calcFunc-choose (math-add (math-sub m n) -1) m)))
           (if (math-evenp (math-trunc m))
               val
             (math-neg val))))
+        ((and (math-negp n) (math-num-integerp n))
+         (if (math-lessp n m)
+             0
+           ;; m≤n<0: (n choose m) = (-1)^(n-m) (-m-1 choose n-m)
+           (let ((val (calcFunc-choose (math-sub (math-neg m) 1)
+                                       (math-sub n m))))
+             (if (math-evenp (math-sub n m))
+                 val
+               (math-neg val)))))
+       ((math-negp m) 0)
        ((and (math-num-integerp n)
              (Math-lessp n m))
         0)
@@ -461,20 +474,23 @@
               (math-choose-float-iter tm n 1 1)))))))
 
 (defun math-choose-iter (m n i c)
-  (if (and (= (% i 5) 1) (> i 5))
+  (while (<= i m)
+    (when (and (= (% i 5) 1) (> i 5))
       (math-working (format "choose(%d)" (1- i)) c))
-  (if (<= i m)
-      (math-choose-iter m (1- n) (1+ i)
-                       (math-quotient (math-mul c n) i))
-    c))
+    (setq c (math-quotient (math-mul c n) i))
+    (setq n (1- n))
+    (setq i (1+ i)))
+  c)
 
 (defun math-choose-float-iter (count n i c)
-  (if (= (% i 5) 1)
+  (while (> count 0)
+    (when (= (% i 5) 1)
       (math-working (format "choose(%d)" (1- i)) c))
-  (if (> count 0)
-      (math-choose-float-iter (1- count) (math-sub n 1) (1+ i)
-                             (math-div (math-mul c n) i))
-    c))
+    (setq c (math-div (math-mul c n) i))
+    (setq n (math-sub n 1))
+    (setq i (1+ i))
+    (setq count (1- count)))
+  c)
 
 
 ;;; Stirling numbers.
diff --git a/test/lisp/calc/calc-tests.el b/test/lisp/calc/calc-tests.el
index 0909756..dce82b6 100644
--- a/test/lisp/calc/calc-tests.el
+++ b/test/lisp/calc/calc-tests.el
@@ -391,6 +391,73 @@ An existing calc stack is reused, otherwise a new one is 
created."
                     (var n var-n) -1 1))
                  8)))
 
+(defun calc-tests--fac (n)
+  (apply #'* (number-sequence 1 n)))
+
+(defun calc-tests--choose (n k)
+  "N choose K, reference implementation."
+  (cond
+   ((and (integerp n) (integerp k))
+    (if (<= 0 n)
+        (if (<= 0 k n)
+            (/ (calc-tests--fac n)
+               (* (calc-tests--fac k) (calc-tests--fac (- n k))))
+          0)    ; 0≤n<k
+      ;; n<0, n and k integers: use extension from M. J. Kronenburg
+      (cond
+       ((<= 0 k)
+        (* (expt -1 k)
+           (calc-tests--choose (+ (- n) k -1) k)))
+       ((<= k n)
+        (* (expt -1 (- n k))
+           (calc-tests--choose (+ (- k) -1) (- n k))))
+       (t  ; n<k<0
+        0))))
+   ((natnump k)
+    ;; Generalisation for any n, integral k≥0: use falling product
+    (/ (apply '* (number-sequence n (- n (1- k)) -1))
+       (calc-tests--fac k)))
+   (t (error "case not covered"))))
+
+(defun calc-tests--check-choose (n k)
+  (equal (calcFunc-choose n k)
+         (calc-tests--choose n k)))
+
+(defun calc-tests--explain-choose (n k)
+  (let ((got (calcFunc-choose n k))
+        (expected (calc-tests--choose n k)))
+    (format "(calcFunc-choose %d %d) => %S, expected %S" n k got expected)))
+
+(put 'calc-tests--check-choose 'ert-explainer 'calc-tests--explain-choose)
+
+(defun calc-tests--calc-to-number (x)
+  "Convert a Calc object to a Lisp number."
+  (pcase x
+    ((pred numberp) x)
+    (`(frac ,p ,q) (/ (float p) q))
+    (`(float ,m ,e) (* m (expt 10 e)))
+    (_ (error "calc object not converted: %S" x))))
+
+(ert-deftest calc-choose ()
+  "Test computation of binomial coefficients (bug#16999)."
+  ;; Integral arguments
+  (dolist (n (number-sequence -6 6))
+    (dolist (k (number-sequence -6 6))
+      (should (calc-tests--check-choose n k))))
+
+  ;; Fractional n, natural k
+  (should (equal (calc-tests--calc-to-number
+                  (calcFunc-choose '(frac 15 2) 3))
+                 (calc-tests--choose 7.5 3)))
+
+  (should (equal (calc-tests--calc-to-number
+                  (calcFunc-choose '(frac 1 2) 2))
+                 (calc-tests--choose 0.5 2)))
+
+  (should (equal (calc-tests--calc-to-number
+                  (calcFunc-choose '(frac -15 2) 3))
+                 (calc-tests--choose -7.5 3))))
+
 (provide 'calc-tests)
 ;;; calc-tests.el ends here