The formula has a 0^0, I’m sure this is what’s causing the NaN.
But the limiting distribution puts all the mass on one point, and this limit is what I was hoping the function would return.
This works for the binomial case, where again the simple formula has an indefinite 0^0, but I guess someone has caught the special case.