lgamma1p {DPQ} | R Documentation |

## Accurate `log(gamma(a+1))`

### Description

Compute

`l\Gamma_1(a) := \log\Gamma(a+1) = \log(a\cdot \Gamma(a)) = \log a + \log \Gamma(a),`

which is “in principle” the same as
`log(gamma(a+1))`

or `lgamma(a+1)`

,
accurately also for (very) small `a`

`(0 < a < 0.5)`

.

### Usage

```
lgamma1p (a, tol_logcf = 1e-14, f.tol = 1, ...)
lgamma1p.(a, cutoff.a = 1e-6, k = 3)
lgamma1p_series(x, k)
lgamma1pC(x)
```

### Arguments

`a` , `x` |
a numeric vector. |

`tol_logcf` |
for |

`f.tol` |
numeric ( |

`...` |
further optional arguments passed on to |

`cutoff.a` |
for |

`k` |
an integer, the number of terms in the series expansion used internally. |

### Details

`lgamma1p()`

is an **R** translation of the function (in Fortran) in
Didonato and Morris (1992) which uses a 40-degree polynomial approximation.

`lgamma1p_series(x, k)`

is Taylor series approximation of order `k`

,
(derived via Maple), which is ```
-\gamma x + \pi^2 x^2/ 12 +
O(x^3)
```

, where `\gamma`

is Euler's constant 0.5772156649. ...

`lgamma1pC()`

is an interface to **R** C API (‘Mathlib’ / ‘Rmath.h’) function.

### Value

a numeric vector with the same attributes as `a`

.

### Author(s)

Morten Welinder (C code of Jan 2005, see R's bug issue
PR#7307) for `lgamma1p()`

.

Martin Maechler, notably for `lgamma1p_series()`

which works
with package Rmpfr but otherwise may be *much* less
accurate than Morten's 40 term series!

### References

Didonato, A. and Morris, A., Jr, (1992)
Algorithm 708: Significant digit computation of the incomplete beta function ratios.
*ACM Transactions on Mathematical Software*, **18**, 360–373;
see also `pbeta`

.

### See Also

### Examples

```
curve(-log(x*gamma(x)), 1e-30, .8, log="xy", col="gray50", lwd = 3,
axes = FALSE, ylim = c(1e-30,1))
sfsmisc::eaxis(1); sfsmisc::eaxis(2)
at <- 10^(1-4*(0:8))
abline(h = at, v = at, col = "lightgray", lty = "dotted")
curve(-lgamma( 1+x), add=TRUE, col="red2", lwd=1/2)# underflows even earlier
curve(-lgamma1p (x), add=TRUE, col="blue") -> lgxy
curve(-lgamma1p.(x), add=TRUE, col=adjustcolor("forest green",1/4),
lwd = 5, lty = 2)
for(k in 1:7)
curve(-lgamma1p_series(x, k=k), add=TRUE, col=paste0("gray",30+k*8), lty = 3)
stopifnot(with(lgxy, all.equal(y, -lgamma1pC(x))))
```

*DPQ*version 0.5-8 Index]