Mercurial > octave-nkf
changeset 6902:d6d19fcc51b0
[project @ 2007-09-14 19:34:12 by jwe]
| author | jwe |
|---|---|
| date | Fri, 14 Sep 2007 19:34:12 +0000 |
| parents | 0baa196d93b5 |
| children | a56ab599ac4c |
| files | scripts/miscellaneous/bincoeff.m |
| diffstat | 1 files changed, 22 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/miscellaneous/bincoeff.m Fri Sep 14 17:16:42 2007 +0000 +++ b/scripts/miscellaneous/bincoeff.m Fri Sep 14 19:34:12 2007 +0000 @@ -75,16 +75,30 @@ b(ind) = 1; ind = ((k > 0) & ((n == real (round (n))) & (n < 0))); - b(ind) = (-1) .^ k(ind) .* exp (gammaln (abs (n(ind)) + k(ind)) ... - - gammaln (k(ind) + 1) - gammaln (abs (n(ind)))); + b(ind) = (-1) .^ k(ind) .* exp (gammaln (abs (n(ind)) + k(ind)) + - gammaln (k(ind) + 1) + - gammaln (abs (n(ind)))); - ind = ((k > 0) & ((n != real (round (n))) | (n >= k))); - b(ind) = exp (gammaln (n(ind) + 1) - gammaln (k(ind) + 1) ... - - gammaln (n(ind) - k(ind) + 1)); - - ## clean up rounding errors + ind = ((k > 0) & (n >= k)); + b(ind) = exp (gammaln (n(ind) + 1) + - gammaln (k(ind) + 1) + - gammaln (n(ind) - k(ind) + 1)); + + ind = ((k > 0) & ((n != real (round (n))) & (n < k))); + b(ind) = (1/pi) * exp (gammaln (n(ind) + 1) + - gammaln (k(ind) + 1) + + gammaln (k(ind) - n(ind)) + + log (sin (pi * (n(ind) - k(ind) + 1)))); + + ## Clean up rounding errors. ind = (n == round (n)); b(ind) = round (b(ind)); - + + ind = (n != round (n)); + b(ind) = real (b(ind)); + endfunction +%!assert(bincoeff(4,2), 6) +%!assert(bincoeff(2,4), 0) +%!assert(bincoeff(0.4,2), -.12, 8*eps)