Number greater than the largest positive floating-point number (1.7977e+308)

조회 수: 11 (최근 30일)
Eric
Eric 2014년 9월 1일
답변: Steven Lord 2016년 5월 18일
I'm trying to calculate the following expression:
(l+m)! / (l!m!)
but when l=m > 85 I get infinity, because the result is greater than the largest positive floating-point number (1.7977e+308).
Is there a way of increasing this number limit, or even a trick to avoid this problem?

이 질문에 답변하려면 로그인하십시오.

채택된 답변

Roger Stafford
Roger Stafford 2014년 9월 1일
There is a much more satisfactory way to compute that expression:
round(prod((L+1:L+m)./(1:m)))
or
round(prod((m+1:m+L)./(1:L)))
whichever you prefer. For the case m = L = 85, the answer would be 9.1448e+49, which is far below the limit for double floating point numbers. You will not encounter difficulties until you reach m = L = 515, and that is unavoidable because the result really is in excess of 'realmax'.
  댓글 수: 4
이전 댓글 2개 표시
Roger Stafford
Roger Stafford 2014년 9월 2일
편집: Roger Stafford 2014년 9월 2일
@Tigo. I will repeat a comment here that I made in John's answer. You say "I figure out that (l+m)! / (l!m!) is equivalent to product_{j=1}^l (m+j)/j. Unfortunately I'm having problems to write this on Matlab." However, you appear not to recognize that matlab's equivalent form of this computation is what I already wrote for you:
prod((m+1:m+L)./(1:L))
For example, if m = 11 and L = 4 both expressions are computing the product
(11+1)/1 times (11+2)/2 times (11+3)/3 times (11+4)/4 = 1365
The two are identical computations, being translations of precisely the same algorithm expressed in two different languages. The 'round' is merely added to this to correct for round-off errors that have occurred during the computation.
Eric
Eric 2014년 9월 2일
편집: Eric 2014년 9월 2일
"Why do you need numbers that big? What is the real world application of this?"
These big numbers come in the numerator of ((l+m)!)/(l!m!), but the actual result is not that big because we still have the denominator. One thing that bugs me is that this same expression runs without problems in Mathematica.

댓글을 달려면 로그인하십시오.

추가 답변 (3개)

Roger Stafford
Roger Stafford 2014년 9월 2일
In case it is of interest, here is an alternate method of computation based on Pascal triangle summations.
m = 25;
L = 11;
x = ones(m+1,1);
for k = 2:L+1
x = cumsum(x);
end
y = x(m+1); % <-- This is the answer
This involves more operations than the other method presented earlier here, but, as opposed to this other method, there are no round-off errors until reaching numbers beyond 2^53 where not all integers can be represented. This is because only sums of integers are being computed in this method.
John D'Errico
John D'Errico 2014년 9월 2일
I'd suggest working in logs, which will be as insensitive as possible to problems with exponents. Only at the very end do you need to exponentiate.
If you truly need to compute a number too large to represent in a double, then you will have no choice in the end. Use either the symbolic toolbox, or a tool like my own HPF toolbox.
  댓글 수: 2
Eric
Eric 2014년 9월 2일
Thanks. But before I try to learn that (I'm a beginner on Matlab), do you have any ideas of how to write this on Matlab?
product_{j=1}^l (m+j)/j
(see http://www.wolframalpha.com/input/?i=product_%7Bj%3D1%7D%5El+%28m%2Bj%29%2Fj)
Roger Stafford
Roger Stafford 2014년 9월 2일
product_{j=1}^l (m+j)/j
I gave you the matlab version of that in the line:
round(prod((m+1:m+L)./(1:L)))
The added 'round' merely corrects for rounding errors occurring in the division process.

댓글을 달려면 로그인하십시오.

Steven Lord
Steven Lord 2016년 5월 18일
(L+M)! / (L! * M!) is just nchoosek(L+M, L) or nchoosek(L+M, M). For most cases this uses the approach that others in this discussion have recommended.
But if L+M is on the order of 5000, the answer will not be exact unless one of L or M is much, much smaller than 5000 and will be infinite (too large to represent as a finite double) if the second input is greater than about 150 or 160.

카테고리

Help Center 및 File Exchange에서 Loops and Conditional Statements에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by