Numerical precision and logarithms

Question

I'm having a problem (at least I think I have a problem) with the following calculation:

ppm <- 20
mDa <- 2
x <- c( 100, 100.002 )

base  <- 1 + ((x * ppm * 1E-6) + (mDa * 1E-3))/x
base
# [1] 1.00004 1.00004
base - 1.00004
# [1]  0.00000e+00 -3.99992e-10

logb( x[2], base[2] ) - logb( x[1], base[1] )
# [1] 1.651291

However, I would have expected that the result is approximately 0.5, since I expected the base to be in both cases to be approximately 1.00004:

logb( x[2], 1.00004 ) - logb( x[1], 1.00004 )
# [1] 0.500005

Although I have no proof at hand, I doubt that the result of logb( x[2], 1.00004 ) - logb( x[1], 1.00004 ) is mathematically correct and I assume that I hit a numerical precision issue. Any ideas how to avoid this problem are highly appreciated.

Edit

What I'm actually trying to do

I need to rescale positive numbers (a, b) -> (a',b') with b > a, such that the difference of two numbers on the new scale d'( a', b' ) = b' - a' is larger 1 iff the difference on the original scale d(a, b) = b -[ a + ( a * ppm * 1E-6) + (mDa * 1E-3)] is larger zero. I know that there might be a problem, because d(a, b) ≠ d(b, a). Typical ranges for the values are a,b ∈ [50, 1500], mDa ∈ [0, 10] and ppm ∈ [1, 50].

Answer 1

When you're taking logarithms of large numbers with a base very close to 1, small differences in that base can lead to noticeable differences in the final value. Your bases differ by 0.0000000004, but that can make a difference with a base very close to 1:

logb(100, 1.0000399996)
# 115132.7
logb(100, 1.00004)
# 115131.6
logb(100, 1.0000400004)
# 115130.4

Answer 2

Try Rmpfr :

Rgames> rfoo<-mpfr(100,100)
Rgames> log100<-log(rfoo)
Rgames> log100
1 'mpfr' number of precision  100   bits 
[1] 4.6051701859880913680359829093676
Rgames> logbase<-log(mpfr(1.0004,100))
Rgames> log100/logbase
1 'mpfr' number of precision  100   bits 
[1] 11515.227896589510924644721707849
Rgames> logbase<-log(mpfr(1.00004,100))
Rgames> log100/logbase
1 'mpfr' number of precision  100   bits 
[1] 115131.55721932987847380223102368

Thus showing that josilber's answer is spot-on.

Answer 3

@josiber is right, log(1.00004) is approximately 0.00004 so you divide by a tiny number and get huge results... So the difference you observe is relatively small.

If you seek better accuracy without extended precision library, you could also try using log1p(x) which compute log(1+x)

baseM1  <- ((x * ppm * 1E-6) + (mDa * 1E-3))/x

then

log( x[2] )/log1p( baseM1[2] ) - log( x[1] )/log1p( baseM1[1] )