The problem lies in calculating the norm. This is numerically zero for too high photon numbers, leading to the infinite value when renormalizing.
Instead, the norm and the sterling formula should be combined to give a better approximation. András (Vukics), haven't we implemented something like that when I last visited Budapest? Have you pushed these changes somewhere?
Last edit: Raimar Sandner 2013-12-13
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
I don’t see how the normalization factor could be combined with coherentElement (which is based on the Sterling formula for high n values), since the latter depends on n while the former doesn’t.
I did some checks a few days ago with coherentElement alone, and it transpired that the problem is not with the Sterling formula, but the exact formula used in the case of n < max_factorial<double>::value. Here, simply an overflow occurs in pow(alpha,n). We could combine the finiteness of pow(alpha,n) into the condition of switching to the Sterling formula, but the problem is that if this happens for too small an n value, then the Sterling formula will not be a good approximation at all.
Ultimately, I think this is not really a bug at all, simply an overflow which can very well be expected.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Ah, I see. Probably the only "solution" then is to use a Gaussian distribution as an approximation to the Poission distribution if this overflow occurs. This might be implemented in C++QED for convenience.
This bug can be closed then?
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The attached figure shows the difference between the exact and Sterling-formula based calculation of 1/\sqrt(n!)
Actually, the Sterling formula does give a good approximation already for small n …
The problem lies in calculating the norm. This is numerically zero for too high photon numbers, leading to the infinite value when renormalizing.
Instead, the norm and the sterling formula should be combined to give a better approximation. András (Vukics), haven't we implemented something like that when I last visited Budapest? Have you pushed these changes somewhere?
Last edit: Raimar Sandner 2013-12-13
I don’t see how the normalization factor could be combined with coherentElement (which is based on the Sterling formula for high n values), since the latter depends on n while the former doesn’t.
I did some checks a few days ago with coherentElement alone, and it transpired that the problem is not with the Sterling formula, but the exact formula used in the case of n < max_factorial<double>::value. Here, simply an overflow occurs in pow(alpha,n). We could combine the finiteness of pow(alpha,n) into the condition of switching to the Sterling formula, but the problem is that if this happens for too small an n value, then the Sterling formula will not be a good approximation at all.
Ultimately, I think this is not really a bug at all, simply an overflow which can very well be expected.
Ah, I see. Probably the only "solution" then is to use a Gaussian distribution as an approximation to the Poission distribution if this overflow occurs. This might be implemented in C++QED for convenience.
This bug can be closed then?
Ticket moved from /p/cppqed/bugs/12/
The attached figure shows the difference between the exact and Sterling-formula based calculation of 1/\sqrt(n!)
Actually, the Sterling formula does give a good approximation already for small n …