evaluating a numerically unstable quotient

Hello!

I am trying to evaluate the following quotient:

   exp((a/c)*(1-c)^k)/[(1-c;1-c)_k * exp(a/c)],

where 'a' is about 1.3, k is a natural number, say 100, c is between 0 and 1 (not inclusive), and (1-c;1-c)_k is the q-Pocchammer symbol. The problem is to evaluate the fraction for c->0, say when c is about 0.01 or even 0.001. In that case the quotient is numerically unstable and MATLAB returns a super large number. What would be a smart way to evaluate the quotient?

Thanks in advance, --Alex

Answer 1

[Matt J]

Looks to me like the quotient should in fact go to infinity as c--->0. I rewrite the quotient expression as follows

   exp(a*f(c))/(1-c;1-c)_k

where

   f(c) = [(1-c)^k  -  1]/c

By L'hopital's rule, f(c)-->-1 as c--->0, so that the numerator goes to exp(-a) asymptotically. The q-Pocchammer quantity in the denominator goes to zero, however.

[Alex Poe]

Right, but say for c=0.001 MATLAB's answer is NaN, while it should be a large but finite number, unlikely outside of the range that MATLAB can handle. The problem is the way the problem is fed to MATLAB, which is I successively compute the exponent in the numerator, one in the denominator, and the QP symbol, and then find the quotient. I think that's not a very intelligent way.

[Matt J]

Well, you could try this. It does give a finite non-NaN result for c=0.001,

 a=1.3;
 k=100;
 c=.001;

 f=@(c) ((1-c)^k  -  1)/c;

 logqpoc=@(aa,qq) sum(log((1-aa*qq.^(0:k-1))));

 >> quotient = exp(a*f(c) - logqpoc(1-c,1-c)),

 quotient =

   2.2211e+89

I can't see the use for such a huge number, however. Are you sure you wouldn't be more interested in the log() of the quotient?