I have the following expression $$ \frac{a^{(k)}}{k!} $$ where $a^{(k)}$ is the rising factorial. Is it better to evaluate it using floating-point arithmetic separately, that is, call a function that returns the numerator and another function that returns the denominator, and then divide? Or should I evaluate it by expanding $$ \prod_{i=1}^k \frac{a+i-1}{i} $$ I tend to think that the second method is more stable because factorials can quickly grow really big, which may lead to a loss of accuracy. But I have never really checked.
$\begingroup$
$\endgroup$
0
Add a comment
|
2 Answers
$\begingroup$
$\endgroup$
2
If the programming environment provides an implementation of the beta function $\mathrm{B}(x,y)$ this computation is straightforward and usually accurate. We have
$$x^{(k)} = \frac{\Gamma(x+k)}{\Gamma(x)}, \ \ \ \ \ k!=\Gamma(k+1), \ \ \ \ \ \mathrm{B}(x,k)=\frac{\Gamma(x)\Gamma(k)}{\Gamma(x+k)}$$
Thus
$$f(x,k) = \frac{x^{(k)}}{k!} = \frac{\Gamma(x+k)}{\Gamma(x)\Gamma(k+1)}=\frac{\Gamma(x+k)}{k\Gamma(x)\Gamma(k)}=\frac{1}{k\mathrm{B}(x,k)}$$
For example $f(9,7) = 6435$. Examples of programming environments that provide the beta function are ISO C++17 or later, Julia, and the SciPy and Boost libraries. Best I can tell from a quick perusal, the beta function is not part of ISO Fortran 2018, but it is presumably available via the NAG and IMSL libraries.
Note the use of the weasel words "usually accurate" above. According to some quick ad-hoc tests some beta function implementations exhibit a maximum error on the order of 100 ulps. The likely cause of this is that they compute $\ln (\mathrm{B}(x,y))$ fully accurate to double precision and then simply exponentiate this intermediate result, thereby suffering from the error magnification property of $\exp$. One way around this is to compute $\ln (\mathrm{B}(x,y))$ with extended precision; for double-precision computation twelve additional bits would suffice to compute $\mathrm{B}(x,y)$ with a maximum error of a few ulps. This is not an unusual issue: similar scenarios occur in implementations of tgamma() and pow() from the C++ standard math library, for example.
-
$\begingroup$ When $k=0$, $f(x,k) = 1$. So, $B(x,k)$ goes to infinity when $k=0$? $\endgroup$– nougakoJan 14 at 20:45
-
$\begingroup$ @nougako $\mathrm{B}(x,k)$ goes to $\infty$ when $k=0$. If you allow $k=0$, you will need to special case that to avoid producing a NaN when using IEEEE-754 floating-point arithmetic. $\endgroup$– njuffaJan 14 at 20:50
$\begingroup$
$\endgroup$
Njuffa already answered satisfactorily, but let me comment that dealing with large numbers does not cause loss of precision in floating-point arithmetic: the error is a relative one, corresponding to a multiplicative term: the result of $ab$ is approximated with $ab(1+\varepsilon)$, with $|\varepsilon|$ smaller than machine precision, independently of the size of $a$ and $b$.
Computing numerator and denominator separately is a bad idea because of a different issue, that of overflow: after $a$ or $k$ are larger than something around 30, numerator and/or denominator become larger than the largest floating point number, which is about $10^{308}$ for float64 / double, and the computation fails. But if overflow does not happen, then the relative error is always bounded by $1+3k\mathsf{u} + O(\mathsf{u}^2)$ (that is, $\mathsf{u}$ for each operation).
answered Jan 13 at 8:37