how to solve n!/(j!*(n-j)!)
Juan Rufilanchas
I have to solve this combinatorial for values of n up to 1000 (and any j),
I have been told that there is a simplification formula (named after a
certain Sterling I thing).
Does anybody has an m file for that or at least the formula?
Any help is really appreciated.
Juan.
Michael Glora
>
> I have to solve this combinatorial for values of n up to 1000 (and any j),
> I have been told that there is a simplification formula (named after a
> certain Sterling I thing).
>
> Does anybody has an m file for that or at least the formula?
I don't know anything about Sterling, but you can try this function
(it will not be extremely fast due to the loops, but it works)
HTH,
Michael Glora
function f = combinatoric(n,k)
% computes n!/(k! (n-k)!)
if 2*k >= n,
%f = prod(k+1:n)/prod(1:n-k);
f = 1;
for i = 1:n-k,
f = f*(k+i)/i;
end
else
%f = prod(n-k+1:n)/prod(1:k);
f = 1;
for i = 1:k,
f = f*(n-k+i)/i;
end
end
The Vaughans
Paul Skoczylas wrote:
> I don't know anything about Stirling, but if I wanted to evaluate
> n!/(n-j)! for large values of n and/or j (but still assuming n>j), I
> would use the gammaln function.
If you want to use some of Paul's erudite advice, please note that he
slightly misread your subject and left out a factor of j! that you had. It
is easy to follow through and correct this omission.
Paul Skoczylas
"The Vaughans" <vaug...@dellnet.com> wrote in message
news:3A02EA2A...@dellnet.com...
OOPS! Yes, I missed that extra factor. That's what happens when people
don't include all the relevant text in the body of the message. (I
really hate it when people do that!) My newsreader prints the subject
in a smaller font than the text--and I obviously didn't read it as
carefully as I could have. (I thought something was missing, as I
remember there being a third term in the log form, but instead of
rereading the question, I just assumed that he was asking about a
different problem than I had seen before.)
Anyway, I think that what I did post was accurate; I just didn't answer
the question that was asked... But since it was probably a homework
problem, the original poster will still have to do a little work, so I
don't feel bad about it at all!
-Paul
Juan Rufilanchas
thanks a lot for all these answers, I am trying them and will see which one
is better for my binomial trees.
Again, thanks a lot.
Juan.
"Juan Rufilanchas" <jru...@retemailxxx.es> wrote in message
news:hYrM5.165$z3.2304@uchinews...
Edric M Ellis
JR> Hi friends, I have to solve this combinatorial for values of n
JR> up to 1000 (and any j), I have been told that there is a
JR> simplification formula (named after a certain Sterling I
JR> thing).
I think that the other response may be more useful to you, but in any
case Stirling's approximation is this:
n! = ( n.^n ) .* exp( -n ) .* sqrt( 2 * pi * n );
This approximation becomes more accurate as n tends to infinity.
What you can also do is take the log of both sides (the resulting
expression is often used in statistical mechanics):
log( n! ) = n .* ( log( n ) - 1 ) + ( .5 * log( 2 * pi * n ) );
[ the last term is often omitted ]
This is useful when you need to deal with factorials of very large
numbers.
Edric.
--
Edric M Ellis
The MathWorks, Ltd.
Tel: +44 (0) 1223 423 200 Ext: 218
Fax: +44 (0) 1223 423 289
Thomas Pieper
using Sterling only works for large n's and j's. Why not doing some low-level
math and hack it into Matlab like this:
=====
function c=fact(n,j);
if j<n/2,
j=n-j;
end;
c=prod((j+1:n)./(1:n-j));
==========
Of course, c=prod(j+1:n)/prod(1:n-j) is more accurate, but the former formula
works for n=1028 and j=514 which is the worst case.
Have fun,
Thomas
Juan Rufilanchas schrieb:
Thomas Pieper
even simpler, why not use the Matlab function nchoosek(n,j)?
Juan Rufilanchas schrieb:
The Vaughans
> I have to solve this combinatorial for values of n up to 1000 (and any j),
> I have been told that there is a simplification formula (named after a
> certain Sterling I thing).
Some kind folks have given you some formulas and code that will probably be
quite helpful. Let me add just a few ideas to help round out your
understanding.
The quantity that you are calculating is a binomial coefficient, as you may
know. What makes it numerically difficult is that each of the three individual
factorials gets huge as "n" and "j" get big. In particular, they will be
vastly bigger than the final result, possibly overflowing your ability to deal
with big numbers, even though the final answer is not that huge.
Here's where Stirling's formula is useful. By expanding each of your
factorials using the formula, you will notice that the hugest pieces nicely
cancel each other out (which is why the final result is not as huge). By
applying the formula, you get rid of the part of the calculation that makes it
challenging from a numerical point of view.
Finally, as others have mentioned, this formula is an approximation that gets
better with larger "n". But what I have not noticed anyone mention is that you
have only been given the zeroth-order term in the expansion. In fact,
Stirling's formula can be used with more terms, allowing you to get accurate
results even for moderate values of "n". You will have to explore numerically
to see how many terms you want to keep.
A reference that I know has Stirling's formula with the first few terms is
"Handbook of Mathematical Function", by Abramowitz and Stegun.
Final note: Stirling's formula is often expressed in terms of the "gamma
function", rather than the factorial. You will need to know that z! =
gamma(z+1).
Hope this helps.
Paul Skoczylas
>
> I have to solve this combinatorial for values of n up to 1000 (and
any j),
> I have been told that there is a simplification formula (named after a
> certain Sterling I thing).
>
> Does anybody has an m file for that or at least the formula?
>
I don't know anything about Stirling, but if I wanted to evaluate
n!/(n-j)! for large values of n and/or j (but still assuming n>j), I
would use the gammaln function.
gamma(m+1)=m!
gammaln(m+1)=log(m!)
Rewrite and manipulate the equation:
A=n!/(n-j)!
log(A)=log(n!/(n-j)!)
log(A)=log(n!)-log((n-j)!)
A=exp(log(n!)-log((n-j)!))
so,
A=exp(gammaln(n+1)-gammaln(n-j+1))
Test this with the following Matlab code:
» n=10;j=6;
» A1=prod(1:n)/prod(1:n-j)
A1 =
151200
» A2=exp(gammaln(n+1)-gammaln(n-j+1))
A2 =
1.5120e+005
Note of course that the prod method won't work for larger n and/or j,
but this example just shows that the gammaln method gives the same
answer. Note also that the gammaln method may not give exact answers
due to the inherent problems with variable precision arithmetic.
However, if you know that your answers are integers (and they will be if
n and j are both integers), you can round your result to the nearest
integer.
-Paul