Improve numeric stability of `dtruncnorm()` when either bound is infinite

[eheinzen]

I'm having troubles where dtruncnorm() returns zero. I'm wondering two things:

1. If there can be a log=TRUE option
2. If so, if you could improve numeric stability by computing the normalization on the log-scale when one bound is infinite.

For instance, dtruncnorm() works well in this case, and matches the base R implementation I have:

> log(truncnorm::dtruncnorm(4.1, a = 4, b = Inf, mean = 3, sd = 2))
[1] -0.587424
> dnorm(4.1, mean = 3, sd = 2, log = TRUE) - pnorm(4, mean = 3, sd = 2, log.p = TRUE, lower.tail = FALSE)
[1] -0.587424
> # or as you have it in your C code
> dnorm((4.1 - 3)/2, mean = 0, sd = 1, log = TRUE) - log(2) - pnorm((4 - 3)/2, mean = 0, sd = 1, log.p = TRUE, lower.tail = FALSE)
[1] -0.587424

But when the standard deviation gets small, my base R implementation works while dtruncnorm() doesn't:

> log(truncnorm::dtruncnorm(4.1, a = 4, b = Inf, mean = 3, sd = 0.01))
[1] -Inf
> dnorm(4.1, mean = 3, sd = 0.01, log = TRUE) - pnorm(4, mean = 3, sd = 0.01, log.p = TRUE, lower.tail = TRUE)
[1] -6046.314

Happy to put through a pull request if you'd prefer.