From 0e00f199cd3599977f75326bb7adc9d70390661e Mon Sep 17 00:00:00 2001 From: =?utf8?q?Mattias=20Engdeg=C3=A5rd?= 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 %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