bug#16999: calc crashes when computation limit is increased

["Mattias Engdegård" <mattiase@acm.org>]

[-- Attachment #1: Type: text/plain, Size: 999 bytes --]

10 sep. 2020 kl. 18.02 skrev Mattias Engdegård <mattiase@acm.org>:

> Let me attempt a patch, and please forgive my reopening the bug.

Here is that patch. Sadly mail to Florian Beck's address bounced but perhaps someone else will see the proposal and subject it to scrutiny. In short:

1. A sign mistake in the original code caused the infinite recursion. Fixed.
2. Some parts of the computation used tail recursion, which limited their applicability in elisp. Changed to use an imperative loop (alas).
3. The binomial coefficients have now been extended for all integral arguments, using definitions from [1]; see also [2] and [3].

[1] M. J. Kronenburg; The Binomial Coefficient for Negative Arguments; https://arxiv.org/abs/1105.3689
[2] D. Loeb, E. Damiani and O. D’Antona; Getting Results with Negative Thinking; https://arxiv.org/abs/math/9502214
[3] David Fowler; The Binomial Coefficient Function; The American Mathematical Monthly, Vol. 103, No. 1 (Jan., 1996), pp. 1-17


[-- Attachment #2: 0001-Calc-fix-binomial-coefficients-for-negative-argument.patch --]
[-- Type: application/octet-stream, Size: 6320 bytes --]

From c451146726adfacdec69247c3da6dc4ebacd1b97 Mon Sep 17 00:00:00 2001
From: =?UTF-8?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=UTF-8
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 c5d4d0837e..ce676d17a9 100644
--- a/lisp/calc/calc-comb.el
+++ b/lisp/calc/calc-comb.el
@@ -445,12 +445,25 @@ calcFunc-choose
 	   (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)
@@ -467,20 +480,23 @@ calcFunc-choose
 	       (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 c8cb97a8bc..5030a55472 100644
--- a/test/lisp/calc/calc-tests.el
+++ b/test/lisp/calc/calc-tests.el
@@ -397,6 +397,73 @@ calc-sum-gcd
                     (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.21.1 (Apple Git-122.3)