Getting NaN while using the ratio of power and factorial

1

With the following piece of code, I get NaN

l=0;
lp=2;
d=200;
a2= ((-sqrt(d))^lp);
a3 =((-sqrt(d))^l);
sum1=0

for m=0:500
    
    a1= ((2^m)*(d^m))/(factorial(m))
    
    mul=a1*a2*a3;
    sum1=sum1+mul;
end

Is there any way to modify the code to get rid of the issue?

[James Tursa]

Cross posted on stackoverflow and answered there over 4 hours ago. If you had simply checked the answer there you could have saved yourself the trouble of posting here and getting essentially the same answer.

https://stackoverflow.com/questions/72961454/nan-in-matlab-while-using-ratio-of-power-and-factorial

Answer 1

[Torsten]

exp(400)*200 is quite large ...

format long
l=0;
lp=2;
d=200;
a2= ((-sqrt(d))^lp);
a3 =((-sqrt(d))^l);
a1 = 1.0;
sum1 = a1;
for m=0:500
    a1 = a1* 2*d/(m+1);
    sum1=sum1+a1;
end
sum1 = sum1*a2*a3

sum1 = 
    1.044293406962028e+176

exp(400)*200

ans = 
    1.044293937952829e+176

Answer 2

[Steven Lord]

m = 500;
d = 200;

The numerator and denominator of your expression (for sufficiently large m) both overflow to inf. Dividing infinity by infinity results in a NaN.

numerator = 2^m*d^m

numerator = Inf

denominator = factorial(m)

denominator = Inf

x = numerator / denominator

x = NaN

One potential approach to avoid this is to avoid explicitly computing 2^m, d^m, or factorial(m).

numeratorVector = repmat(2*d, 1, m); % prod(numeratorVector) would effectively give numerator
denominatorVector = 1:m; % prod(denominatorVector) would effectively give denominator
xVector = numeratorVector./denominatorVector;
format longg
x2 = prod(xVector)

x2 = 
     8.78187252741829e+166

Let's check symbolically.

numeratorSymbolic = sym(2*d)^m;
vpa(numeratorSymbolic) % Pretty big

ans = 

1.07150860718626732094842504906e+1301

denominatorSymbolic = factorial(sym(m));
vpa(denominatorSymbolic) % Also pretty big

ans = 

1.220136825991110068701238785423e+1134

x3 = vpa(numeratorSymbolic/denominatorSymbolic) % Big but not quite as big as above

x3 = 

8.7818725274182843596325298923868e+166

You could be a little more sophisticated / clever if you wanted (preemptively cancelling out factors of 2 in numeratorVector by dividing even values in denominatorVector by 2.) Or you could keep track of x for each value of m then figure out what you need to multiply it by to get x for the next value of m.

Answer 3

[Benjamin Thompson]

The factorial function output increases very fast as input increases. See "doc factorial" for details. The output is "inf" for input of 171 or larger.

[Steven Lord]

But 400^m overflows for m > 118 so this whole expression becomes infinite or undefined beyond that point.

d = 200;
(2*d)^118

ans = 1.1043e+307

(2*d)^119

ans = Inf

Answer 4

[Vijeta Singh Yadav]

Data overflow and underflow problem exists in the code

visit this link for avoid data overflow and underflow problems.