From 0e00f199cd3599977f75326bb7adc9d70390661e Mon Sep 17 00:00:00 2001
From: =?utf8?q?Mattias=20Engdeg=C3=A5rd?= <mattiase@acm.org>
Date: Fri, 11 Sep 2020 11:43:15 +0200
Subject: [PATCH] Calc: fix binomial coefficients for negative arguments
 (bug#16999)
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

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 2efeb7f0f22..f7e29c6e52c 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 09097564688..dce82b6f536 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
 
-- 
2.39.5