Mercurial > octave-libtiff
changeset 30217:f5aba4d7e651
nchoosek.m: New algorithm with more precision when N, K are scalars (bug #61199)
* nchoosek.m: Prepare numerator and denominator of binomial coefficient. Use
gcd() to find and filter out common factors.
author | Rik <rik@octave.org> |
---|---|
date | Tue, 28 Sep 2021 17:03:25 -0700 |
parents | c2132e0147ff |
children | 750de16fd35c |
files | scripts/specfun/nchoosek.m |
diffstat | 1 files changed, 30 insertions(+), 3 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/specfun/nchoosek.m Mon Sep 27 20:30:57 2021 -0700 +++ b/scripts/specfun/nchoosek.m Tue Sep 28 17:03:25 2021 -0700 @@ -111,10 +111,37 @@ n = numel (v); if (n == 1 && isnumeric (v)) - ## Improve precision at next step. + ## Improve precision over direct call to prod(). + ## Steps: 1) Make a list of integers for numerator and denominator, + ## 2) filter out common factors, 3) multiply what remains. k = min (k, v-k); - C = round (prod ((v-k+1:v)./(1:k))); - if (C*2*k*eps >= 0.5) + + ## For a ~25% performance boost, multiply values pairwise so there + ## are fewer elements in do/until loop which is the slow part. + ## Since Odd*Even is guaranteed to be Even, also take out a factor + ## of 2 from numerator and denominator. + if (rem (k, 2)) # k is odd + numer = [(v-k+1:v-(k+1)/2) .* (v-1:-1:v-(k-1)/2) / 2, v]; + denom = [(1:k/2) .* (k-1:-1:(k+1)/2) / 2, k]; + else # k is even + numer = (v-k+1:v-k/2) .* (v-k/2+1:v) / 2; + denom = (1:k/2) .* (k/2+1:k) / 2; + end + + # Remove common factors from numerator and denominator + do + for i = numel (denom):-1:1 + factors = gcd (denom(i), numer); + [f, j] = max (factors); + denom(i) /= f; + numer(j) /= f; + endfor + denom = denom(denom > 1); + numer = numer(numer > 1); + until (isempty (denom)) + + C = prod (numer); + if (C > flintmax) warning ("nchoosek: possible loss of precision"); endif elseif (k == 0)