[Rd] RFC: lchoose() vs lfactorial() etc

[Martin Maechler]

>>>>> "MM" == Martin Maechler <maechler at stat.math.ethz.ch>
>>>>>     on Tue, 15 Dec 2009 14:54:18 +0100 writes:

>>>>> "PS" == Petr Savicky <savicky at cs.cas.cz>
>>>>>     on Tue, 15 Dec 2009 10:52:43 +0100 writes:

    PS> On Tue, Dec 15, 2009 at 09:49:28AM +0100, Martin Maechler wrote:
    >>> lgamma(x) and lfactorial(x) are defined to return
    >>> 
    >>> ln|Gamma(x)| {= log(abs(gamma(x)))} or  ln|Gamma(x+1)| respectively.
    >>> 
    >>> Unfortunately,  we haven't chosen the analogous definition for 
    >>> lchoose().
    >>> 
    >>> So, currently 
    >>> 
    >>> > lchoose(1/2, 1:10)
    >>> [1] -0.6931472 -2.0794415        NaN -3.2425924        NaN -3.8869494
    >>> [7]        NaN -4.3357508        NaN -4.6805913
    >>> Warning message:
    >>> In lchoose(n, k) : NaNs produced
    >>> > 
    >>> 
    >>> which (the NaN's) is not particularly useful.
    >>> (I have use case where I really have to workaround those NaNs.)
    >>> 
    >>> I herebey propose to *amend* the definition of lchoose() such
    >>> that it behaves analogously to  lgamma() and lfactorial(),
    >>> i.e., to return
    >>> 
    >>> log(abs(choose(.,.))
    >>> 
    >>> Your comments are welcome.
    >>> Martin Maechler, ETH Zurich

    PS> I do not have a strong opinion, whether lchoose() should be log(choose(.,.))
    PS> or log(abs(choose(.,.))), although i would slightly prefere the latter (with abs()).
    PS> One of the reasons is that this would simplify the code of the function
    PS> lchoose() in src/nmath/choose.c. It produces the NA by explicit testing
    PS> the sign, so this test could be simply removed.

    PS> Let me also point out that the current implementation seems to contain
    PS> a bug, which causes that it is neither log(choose(.,.)) nor log(abs(choose(.,.))).

    PS> x <- choose(1/2, 1:10)
    PS> y <- lchoose(1/2, 1:10)
    PS> cbind(choose=x, lchoose=y, log.choose=log(abs(x)))
    PS> choose    lchoose log.choose
    PS> #  [1,]  0.500000000 -0.6931472 -0.6931472
    PS> #  [2,] -0.125000000 -2.0794415 -2.0794415
    PS> #  [3,]  0.062500000        NaN -2.7725887
    PS> #  [4,] -0.039062500 -3.2425924 -3.2425924
    PS> #  [5,]  0.027343750        NaN -3.5992673
    PS> #  [6,] -0.020507812 -3.8869494 -3.8869494
    PS> #  [7,]  0.016113281        NaN -4.1281114
    PS> #  [8,] -0.013092041 -4.3357508 -4.3357508
    PS> #  [9,]  0.010910034        NaN -4.5180723
    PS> # [10,] -0.009273529 -4.6805913 -4.6805913

    PS> So, besides x[1], we get NA for those cases, when choose() is positive.
    PS> The reason seems to be the test at line 89 of src/nmath/choose.c, which is

    PS> if (fmod(floor(n-k+1), 2.) == 0) /* choose() < 0 */
    PS> return ML_NAN;

    PS> but it should be negated. I tested it using

    PS> if (fmod(floor(n-k+1), 2.) != 0) /* choose() < 0 */
    PS> return ML_NAN;

    MM> Indeed !!! 
    MM> Pretty  "interesting" that this thinko has not been detected for
    MM> quite a few years....

    MM> This, of course, is an even more compelling reason to implement
    MM> the change of return  log(abs(choose(.,.)),
    MM> and at the moment, I'd even plan to "backport" that to R "2.10.1
    MM> patched", as the current behavior is clearly bugous.

This has now happened;
svn revisions  50775 and 50776  {R-devel and R-patched}.

Martin