josdejong / mathjs

An extensive math library for JavaScript and Node.js
https://mathjs.org
Apache License 2.0
14.37k stars 1.24k forks source link

factorial(0^(-1)) #444

Closed sonnyk22 closed 9 years ago

sonnyk22 commented 9 years ago

It should not be able to evaluate this factorial(0^(-1)) because it is out of range or error, but MathJS evaluates it to 2.5066283.

FSMaxB commented 9 years ago

This is quite strange indeed. This might be related to how the factorial of real numbers get's calculated using the Gamma function.

This looks like a real bug.

FSMaxB commented 9 years ago

Looking at the source code it seems like this value was intentionally added, see https://github.com/josdejong/mathjs/blob/7b33c1e25473540266745db51c7011911910bf15/lib/function/probability/factorial.js#L34-L35

It's:

if (n === Number.POSITIVE_INFINITY) {
       return Math.sqrt(2 * Math.PI);
}

I'm not aware of the reason, but I'm not a mathematician.

FSMaxB commented 9 years ago

This change has apparently been introduced in commit 05c6dcaa. @BigFav, can you explain why the factorial of positive infinity should be sqrt(2*pi)?

ahrvoje commented 9 years ago

For POS_INF result should also be POS_INF, there is no reason to return sqrt(2pi), it's just plain wrong!

josdejong commented 9 years ago

One thing I noticed here is that the configuration option predictable is not applied by factorial and gamma, these should return NaN for invalid values. We need to fix that.

As for factorial(infinity) returning sqrt(2 pi), this is mathematically correct. It was introduced in PR https://github.com/josdejong/mathjs/pull/267. factorial(infinity) is the product of an infinite series, which is quite a topic in mathematics. Infinity is a very special number, and our "normal" rules for working with number do not always apply. You can google for it. For example read this question: http://math.stackexchange.com/questions/390435/factorial-of-infinity

I agree though that this may be very confusing and I'm not sure if it makes sense in practice. Just for your idea: wolfram gives undefined as result for factorial(0^(-1)), and returns Infinity as result for factorial(Infinity). Matlab also returns Inf for factorial(Inf). @BigFav what do you think?

ahrvoje commented 9 years ago

Wolfram returns POS_INF for POS_INF input, and that is a correct way to do it, mathjs should do the same. Classification of infinities is a separate problem, and that is where different outcomes for fact(INF), fact(1/0), fact(0^(-1)) are possible, maybe, but I believe that is a whole another discussion. Even in that case, result could be undefined, at best, never a finite value.

Also, one should be very very careful while evaluating inifinity expressions using zeta function. stackexchange post you refering to is a trick, one of many, to get a finite value using this technique. But there is an error in reasoning, the correct result is Infinity! There are many funny posts describing the process of getting this incorrect result. http://physicsbuzz.physicscentral.com/2014/01/redux-does-1234-112-absolutely-not.html https://plus.maths.org/content/infinity-or-just-112 https://www.youtube.com/watch?v=w-I6XTVZXww https://en.wikipedia.org/wiki/Zeta_function_regularization

There is no point in returning zeta-regularized results. If you do it here, the entire mathjs should be corrected to return zeta-reg result for every function it implements, and you could end up with no function returning infinite result, all would be regularized, with most users not being aware and not understanding what is going on.

In these extreme cases one should avoid using some 'super smart' definitions of factorial, and go back to basics. You get:

factorial(inf) = limit ( product(1, n), n -> inf) = inf

It's that simple!

BigFav commented 9 years ago

I am truly sorry for not commenting on this earlier. I chatted with someone else about it, and thought I got around to it here.

I think I got ahead of myself when I added this, as this was pushed around the time I learned that there was work done involving the Casimir force which, in my opinion, heavily imply that is the correct physical interpretation of the 1 + 2 + 3 + 4 ... = -1/12, but even this doesn't generalize to all of the other transformations into regularized results (like the one performed here). As an aside, if you are actually interested in this stuff, and have a high level of mathematical knowledge, you can check these out:

http://physics.stackexchange.com/questions/26877/regularization-of-the-casimir-effect https://terrytao.wordpress.com/2010/04/10/the-euler-maclaurin-formula-bernoulli-numbers-the-zeta-function-and-real-variable-analytic-continuation/ (This is the sited article in the above page)

I think this function should behave in an expected manner, so I think it should be changed as well. When most people use factorial with infinity, they almost definitely are not looking for a zeta-regularized result, and that it is unreasonable to expect the user to not only have this knowledge, but also want this usage by default. I will change it.

ahrvoje commented 9 years ago

Yes, that's the correct way to do it. I'm theoretical physicist, rather familliar with regularization techniques, the very reason I commented on the issue. Regularization is used (Casimir force calculation included) to avoid inifinite values at the intermediate calculation steps while obtaining the correct final result, if and when these inifinities cancel each other out correctly. Regularized infinity is not the correct limit of the expression, and it should not be used in the library.

As you (BigFav) point out, the function should behave in an expected manner returning POS_INF, and must be changed. Otherwise you will get incorrect limits.

BigFav commented 9 years ago

I think we can all agree limits is the way to go for this. It has been fixed in the develop branch, so it will be updated during the next release, and when it is I believe this issue will auto-close.

josdejong commented 9 years ago

Thanks @BigFav! I will do a release tonight.

josdejong commented 9 years ago

Done :)