Search code examples
gnuplotlogarithmrounding-error

gnuplot: how to get correct order of magnitude?


This question/issue is maybe related somehow to this topic.

If you type:

print log10(1e7) you will get 7.0.

print int(log10(1e7)) you will get 7.

However, if you type

print log10(1e6) you will get 6.0.

print int(log10(1e6)) you will get 5.

These are probably rounding errors related to log10 and cannot(?) be avoided.

Because if you type

print sprintf("%.20e",log10(1e6)) gives 5.99999999999999911182e+00

print sprintf("%.20e",log10(1e7)) gives 7.00000000000000000000e+00

You can extend and summarize this as a plot: Code:

### power problem in gnuplot
reset session
set colorsequence classic
set key left
set samples 41

set xrange[-20:20]
plot int(log10(10**x)) w lp pt 7,\
    x w lp pt 7
### end of code

Result:

enter image description here

You will see that in irregular distances there are differences between the expected and the obtained result.

So, I am still missing a function which gives me always the correct order of magnitude. Maybe rounding all numbers first to 15 decimal places? Any other ideas?


Solution

  • assuming that you are not dealing with more than 12-15 digits (or as @Ethan says more than 15-16 digits in a 64 bit system is nonsense anyway), the following function should give the correct order of magnitude as integer. I just tested a few examples and compared it with other "straightforward" methods. Please prove the function right or wrong.

    ### get the correct power of a number with gnuplot
    
    CorrectPower(n) = floor(log10(n*(1+1e-15)))
    IncorrectPower1(n) = floor(log10(n))
    IncorrectPower2(n) = floor(gprintf("%T",n))
    
    Numbers = "1e-6 1e-4 0.001 0.01 1e-2 1000 1000000 -1e-6 -1e-9 0.99 95 990"
    
    print " Number    cP  icP1 icP2"
    do for [i=1:words(Numbers)] {
        n = word(Numbers,i)
        print \
            sprintf("%7s:%5d%5d%5d", n, CorrectPower(n), IncorrectPower1(n), IncorrectPower2(n))
    }
    ### end of code
    

    Result:

     Number    cP  icP1 icP2
       1e-6:   -6   -5   -6
       1e-4:   -4   -3   -4
      0.001:   -3   -2   -3
       0.01:   -2   -1   -2
       1e-2:   -2   -1   -2
       1000:    3    2    3
    1000000:    6    5    6
      -1e-6:   -6   -5   -6
      -1e-9:   -9   -8   -9
       0.99:   -1    0    0
         95:    1    1    2
        990:    2    2    3
    

    Addition: For what it's worth, another function to get correct power as integer:

    CorrectPower2(n) = int(sprintf("%.15e",n)[19:])