R: how to have gamma() return the actual number instead of Inf

The issue is that you cannot represent it with double precision:

gamma(200) # too large value
#R> [1] Inf
lgamma(200) # but log is not
#R> [1] 857.9337
exp(857) # the issue!
#R> [1] Inf
.Machine$double.xmax # maximum double value
#R> [1] 1.797693e+308
gamma(171) # almost there!
#R> [1] 7.257416e+306

You can work with the log of the gamma function instead using lgamma. Otherwise you will need to use a third party library which has higher precision than R's floating points.

A google search suggests that the Rmpfr::igamma function might be what you want if you cannot work with the log of the gamma function:

Rmpfr::igamma(171, 0)
#R> 1 'mpfr' number of precision  53   bits 
#R> [1] 7.257415615307999e+306
Rmpfr::igamma(200, 0)
#R> 1 'mpfr' number of precision  53   bits 
#R> [1] 3.9432893368239526e+372

Using lgamma as proposed by Benjamin Cristoffersen, you can calculate the significand and the exponent (base 10) as individual variables:

(res <- gamma(100))
9.332622e+155

# Natural logarithm of result
(ln_res <- lgamma(100))
359.1342

# log base 10 of result
(log10_res <- ln_res/log(10))
155.97

# decimal part of the number above, raised to the 10th power
(significand_res <- 10 ^ (log10_res %% 1))
9.332622

# non-decimal part
(exp_res <- log10_res %/% 1)
155

For gamma(200) this returns: 3.9432 * 10 ^ 372