Thread Subject: How to calculate large number

the largest integer in matlab is around 170!. I wonder how to calculate the question below which has several parts larger than 170!:

b(k) = the summary from k=0 to n: (((-1)^k)/(k!))*((n-k)^k)* (e^((n-k)*2))*(2^k)

assume n=5000.

You may want to give John D'Errico's FEX submission a try on this one.

http://www.mathworks.com/matlabcentral/fileexchange/22725




byiGy[[ibfnyc]hyfiny[\hb?o[sy_[[p(4Soe:y!ni\a]ijg[_i_!jyg_c

netbsd8 wrote:
> the largest integer in matlab is around 170!. I wonder how to calculate the question below which has several parts larger than 170!:

> b(k) = the summary from k=0 to n: (((-1)^k)/(k!))*((n-k)^k)* (e^((n-k)*2))*(2^k)

> assume n=5000.

For n = 5000, you will need to evaluate the expression to approximately 1368 decimal places
in order to get the first decimal place correct in the answer; my trials appear to indicate
that below 1365 decimal places for the calculation, even the sign of the result may
be wrong.

--
.signature note: I am now avoiding replying to unclear or ambiguous postings.
Please review questions before posting them. Be specific. Use examples of what you mean,
of what you don't mean. Specify boundary conditions, and data classes and value
relationships -- what if we scrambled your data or used -Inf, NaN, or complex(rand,rand)?

Thanks. It seems good for integers, but how to deal with float number? for example: exp(5000)? Thanks.

netbsd8 <netbsd8@gmail.com> wrote in message <26997467.1233209248647.JavaMail.jakarta@nitrogen.mathforum.org>...
> Thanks. It seems good for integers, but how to deal with float number? for example: exp(5000)? Thanks.

Hi,
why dont you do the exp(log()) trick!
If you rewrite every term in the sum to
exp(log((((-1)^k)/(k!))*((n-k)^k)* (e^((n-k)*2))*(2^k)))

it can be written as

exp( log(-1^k) + log((n-k)^k) + log(exp((2*(n-k))) + log(2^k) - log(k!))

if you then apply the rules for taking the logarithm of powers and if you now that gamma(n+1) = n! and gammaln(n+1) = log(n!) then you should avoid dealing with large numbers. At least it should reduce some problems.
BTW the same holds for integers and real numbers.

Hope that helps!

Thanks a lot!

if I do as below:

if mod(k,2)=0
     bn=exp(k*log(n-k)+(n-k)*2*log(e)+k*log2-gammaln(k-1))
else
     bn=-exp(k*log(n-k)+(n-k)*2*log(e)+k*log2-gammaln(k-1))
end

then, when k is small, for example, k=2, I still need to calculate exp(>5000).

Any suggestion?

Yes you are right, it is still tricky!
 My first idea was simplify the whole thing

bn = exp( 2*n + cn) where cn is rest in the exponent. For huge n this will dominate the exponent. Because you have an alternating series, you have something like

... + exp(2*n + c(n)) - exp(2*(n+1) + c(n+1)) + ... . For huge n c(n) will be similar to c(n+1), hence you have something like
... + exp(2*n + c(n)) - exp(2*(n+1) + c(n)) + ...

the next step was to do a Taylor expansion around x = 2*n+c(n)
by this I ended up with something like
(1-2*x+2*n+2*x^2-4*x*n+2*n^2)*(exp(2*x+2+c)-exp(2*x+c))
and therefore the same mess like before.

Sorry, I had my last lectures in statistical physics some year ago. I have no clue how to solve this problem (Luckily it isnt my problem ;-) ).

You might search google for "exponential sum". Maybe this can help you!

BTW what kind of sum is it, related to number theory?

Sorry mate

netbsd8 <netbsd8@gmail.com> wrote in message <201661.1233201430355.JavaMail.jakarta@nitrogen.mathforum.org>...
> the largest integer in matlab is around 170!. I wonder how to calculate the question below which has several parts larger than 170!:
>
> b(k) = the summary from k=0 to n: (((-1)^k)/(k!))*((n-k)^k)* (e^((n-k)*2))*(2^k)
>
> assume n=5000.

I'll assume that you don't really need the exact
integer solution. My vpi tools will not do exponentials
as you found out anyway, although an interesting
question that I might take up is how to solve for
the correct integer part of exp(x), for large integer
values of x. Why anybody would like to see the
integer part of a huge exponential, I don't know.
But then I'm often surprised at what people use
in my submissions to the FEX, and how they use
them. Regardless, assuming that you want a double
precision result...

I'd definitely recommend logs to compute the
terms within the sum. Ignoring the sign term,
we have

log(1/(k!)) + k*log(2*(n-k)) + 2*(n-k)

Next, for small k, log(1/k!) is trivial to compute. But
for large k, either use a Stirling series, or gammaln.
gammaln makes sense here. Next, we need to know
how big these terms get.

n = 5000;
term = zeros(1,n+1);
term(n+1) = -inf;
for k = 0:(n-1)
  term(k+1) = k*log(2*(n-k)) + 2*(n-k) - gammaln(k+1);
end

Note that when k == n, we get a log(zero) in the
expression. How large do these terms get?

format long g
maxterm = max(term)
maxterm =
          11084.2586420162

Do the sum, subtract out the maximum from each
term. Do it by summing the positive members of
the series, then the negative members. Expect to
see some nasty subtractive cancellation anyway.

posterms = sum(exp(term(1:2:n) - maxterm))
posterms =
          30.7768674230547

negterms = sum(exp(term(2:2:n) - maxterm))
negterms =
           30.776867423051

So the sum is exp(maxterm)*(posterms-negterms).
Since maxterm is on the order of 11000, this is
still a problem. But a question is if those positive
and negative terms cancel each other out. They
are too close for my comfort.

Having convinced myself that this sum is
probably a homework assignment and actually
reduces to something I might find in the
literature somewhere, I'll stop here. Google
may be your best friend now.

John

Thanks again.

The equation comes from a paper related with M/D/1/N queue modeling

I tried to plot the graph when n<=200 which could be calculated by Matlab. at the beginning, the bn's quantity keep positive and increased a little, then it fluctuated between very large negative number and positive number, then it will be convergent to 1.924e+047 after n>=88. surprising!

netbsd8 wrote:
> I tried to plot the graph when n<=200 which could be calculated by Matlab. at the beginning,
> the bn's quantity keep positive and increased a little, then it fluctuated between very
> large negative number and positive number, then it will be convergent to 1.924e+047
> after n>=88. surprising!

You should not have been even slightly surprised if you had read my earlier reply
in which I wrote:

>> For n = 5000, you will need to evaluate the expression to approximately 1368 decimal places
>> in order to get the first decimal place correct in the answer; my trials appear to indicate
>> that below 1365 decimal places for the calculation, even the sign of the result may
>> be wrong.

The exact number of decimal places you will need for the calculation will depend upon
exactly how you break up the calculation. If you break it up as

sum( evalf((-1)^k/k!) * evalf((n-k)^k) *
  evalf( evalf(exp((n-k)*2)) * evalf(2^k)) ), k=0..n)

then for n from 1 to 200, each of those evalf's will need to be carried
out at a minimum of 54 decimal places in order to get the signs all
positive and to get the first decimal place correct for n=200. Evaluating
at 53 decimal places would get all positive values but the first
decimal place would be quite wrong for n=200. Evaluating at 50 decimal
places would have some of the results come out negative.


By graphing the log of the expression calculated at high precision
and eyeballing the result, it appears that the log of the expression
is relatively linear for increasing n (but not necessarily exactly
linear, I didn't plot fine enough grained to tell). *Possibly* for your
purposes you could replace the expression by an approximation of the form
A + B * exp(C * k). This still has the problem that C * k comes out
on the order of 300 (that is, that you would have exp(300)) for n = 200,
but with the substitution it might become practical to transform the
rest of the problem to one involving logs.

--
.signature note: I am now avoiding replying to unclear or ambiguous postings.
Please review questions before posting them. Be specific. Use examples of what you mean,
of what you don't mean. Specify boundary conditions, and data classes and value
relationships -- what if we scrambled your data or used -Inf, NaN, or complex(rand,rand)?