0

I want to represent the following number using the log function:

2.5e-600/1.7e-500

here is what I did on paper, and what I would like to automate on R for any given number:

log(2.5e-600/1.7e-500) = log(2.5e-600)-log(1.7e-500)
                      = log(2.5)-600*log(10) - log(1.7) + 500*log(10)
                      = -229.8728

However, I am thinking that on R it won't be as straight forward to go from log(10^-600) to -600*log(10). Because R evaluates the inside expression first then applies the log function which gives -Inf instead of -1381.511

My question is how can I remedy to that problem? I am thinking that maybe there is a function that would allow me to retrieve the exponent part of a number ? same question for going from log(2.5e-600) to log(2.5)-600*log(10)

Imlerith
  • 479
  • 1
  • 7
  • 15

2 Answers2

2

Rs usual numeric format (double) will leave you SOL here, because the numbers you write down are already equal to 0 for R. (Type it out in console: ´2.5e-600´ to get [1] 0 as an answer)

You could take a look at arbitrary precision arithmetic packages for R. A quick google search brought up https://cran.r-project.org/web/packages/Rmpfr/vignettes/Rmpfr-pkg.pdf

> log(mpfr("2.5e-300")/mpfr("1.7e-500"))
1 'mpfr' number of precision  20   bits 
[1] 460.90283

Edit: now that the representation is clear, still use mpfr like this:

N <- mpfr(representable_part) * mpfr("1e-308")^reduction_number

For example

> mpfr(2.5) * mpfr("1e-308")^13
1 'mpfr' number of precision  17   bits 
[1] 2.500098e-4004
AlexR
  • 2,412
  • 16
  • 26
  • that's the point: because R gives 0 for 2.5e-600, I want to represent it using the log function – Imlerith Aug 11 '16 at 15:27
  • Please clarfiy what you mean by "represent" here. You can have R take the number as a string, for example. – AlexR Aug 11 '16 at 15:28
  • I meant rescale the number to avoid underflow – Imlerith Aug 11 '16 at 15:29
  • 1
    My question is **in what form** do you already **have** the number in R? Do you have it as a string (see my edit for a demo of said Package with string representations of the numbers) – AlexR Aug 11 '16 at 15:31
  • Allright, take a look at my edit for how to convert that to the proper mpfr number with wich you can calculate as usual. – AlexR Aug 11 '16 at 15:41
2

For this special case, you can do it fairly straightforwardly without using mpfr:

x1 <- "2.5e-600"
x2 <- "1.7e-500"

Split value into mantissa and exponent (log10):

get_num <- function(x) {
    as.numeric(strsplit(x,"e")[[1]])
}
val1 <- get_num(x1)
val2 <- get_num(x2)

Now take the ratio of the mantissas (log(v1)-log(v2)), but subtract the exponents:

log(val1[1])-log(val2[1])+log(10^(val1[2]-val2[2]))
## [1] -229.8728
Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • 2
    If you go down this route, be sure to simplify the last log-Term to `(val1[2]-val2[2])*log(10)` to be even more robust to under/overflow. Also, using regex for the splitting would increase robustness: `([+-]?[\d]+(?:\.\d+))[eE]([+-]?\d+)` will capture only well-formed numbers with groups 1 and 2 being the respective numbers. – AlexR Aug 11 '16 at 16:59