Search code examples
rfloating-pointlogarithm

Numerical precision and logarithms


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].


Solution

  • 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