R devel - May 2020 - pbirthday() for larger number of classes

Hi,

pbirthday(, coincident = 2) starts to issue warnings (see (*) below)
for larger number of classes (R 4.0.0, R-devel
./src/library/stats/R/birthday.R:47).

The default coincident = 2 is computed as 1 - prod((c:(c - n +
1))/rep(c, n)) where c = classes.
Using exp(log(...)), one can derive the return value if(n > 0)  1 -
exp(sum(log1p(-(0:(n-1))/c))) else 0.
Simplifying this a bit further one obtains if(n >= 2) 1 -
exp(sum(log1p(-(1:(n-1))/c))) else 0.
For large c, sum(log1p(-(1:(n-1))/c)) is close to 0, so a more robust
version would be
to return if(n >= 2) -expm1(sum(log1p(-(1:(n-1))/c))) else 0 in the
default case 'coincident = 2'
(internally: if (k == 2) ...).

## Auxiliary function *just* considering 'coincident = 2'
pbirthday2 <- function(n, classes = 365) {
    c <- classes # as pbirthday()
    if(n >= 2) -expm1(sum(log1p(-(1:(n-1))/c))) else 0 # return value
suggested for the case 'if (k == 2) ...'
}

pbirthday (3, classes = 2) == pbirthday2(3, classes = 2) # identical

pbirthday (3, classes = 2^53) == pbirthday2(3, classes = 2^53) # not
identical anymore...
stopifnot(all.equal(pbirthday (3, classes = 2^53), pbirthday2(3,
classes = 2^53))) # ... but numerically equal

pbirthday (3, classes = 2^54) # warnings start to appear (*)
pbirthday2(3, classes = 2^54) # fine

pbirthday (3, classes = 2^56) == 0 # numerically indistinguishable from 0
pbirthday2(3, classes = 2^56) # fine

Cheers,
M

... and one should include the pigeonhole principle:

pbirthday2 <- function(n, classes = 365) {
    c <- classes # as pbirthday()
    if(n >= 2) {
        if(n > classes) 1 else -expm1(sum(log1p(-(1:(n-1))/classes)))
    } else 0
}