Why is 10^(log10(x)) not x?

Hello,

When I enter the equation x = 10^(log10(x)), I should get exactly the same value for x again. However, it differs slightly in the decimal places. It is important for my calculation, that it ends up with exactly the same value. Do I have a wrong input in matlab?

a = 14.25;
a2 = 10^(log10(a));
disp(cell2mat(compose('%10.20f',a2)))
    14.24999999999999822364

Has anybody an explanation? Thanks a lot!

Answer 1

[Stephen23]

"Has anybody an explanation?"

The simple explanation is: log10(a) is not exactly representable using a binary floating point number.

Any operations on binary floating point numbers can accumulate floating point error. It is the author's responsibility to understand this and take into account the finite precision of all numeric calculations and the results.

a = 14.25;
b = log10(a)

b = 1.1538

Question: can the value b be exactly represented by a finite binary number? (hint: no)

Question: what happens to this floating point error when you perform operations on that value? (hint: it accumulates)

"It is important for my calculation, that it ends up with exactly the same value."

That is not how numeric computing works. Computers do not have infinite memory, so you cannot expect them to store infinite precision numeric data. It is exactly the same as if I ask you to write 1/3 as a decimal fraction on a piece of paper. Then perform some arithmetic operations on the number that you have written down (not on the fraction 1/3): does the result of those operations have infinite precision? (hint: no, it is limited to what you could write down).

If you really want to work with exactly the same value then you could use symbolic variables/operations. But this is more complex and much slower than numeric computations.

s = sym(14.25)

s = 

t = simplify(10^(log10(s)))

t = 

[Juditha Schmidt]

Thanks a lot! I don't have a mathematical background and this explanation really helped me!

[Steven Lord]

One explanation I think explains roundoff well is to ask the person to divide 1 by 3 using pencil and paper. They may write as many decimal places of the result as they want, but only the decimal places they explicitly write in this step can be used in the next.

Now multiply that result by 3. In exact arithmetic, (one divided by three) times three is exactly one. But due to the roundoff of the result caused by not writing the (infinite number) of decimal places their result of (1 divided by 3) times 3 is just slightly smaller than one.

Answer 2

[Walter Roberson]

a = 14.25;
a2 = 10^(log10(a));
for K = 1 : 100
    a3 = 10^(log10(a+K*eps(a)));
    if a3 ~= a2; break ;end
end
fprintf('a2 = %10.20f\n', a2);

a2 = 14.24999999999999822364

fprintf('a3 = %10.20f after %d*eps\n', a3, K);

a3 = 14.25000000000000532907 after 2*eps

fprintf('a2 difference = %10.20f\n', abs(a2 - a));

a2 difference = 0.00000000000000177636

fprintf('a3 difference = %10.20f\n', abs(a3 - a));

a3 difference = 0.00000000000000532907

You can see from this that the value reconstructed into a2 is the closest that the 10^log10 process can get.

The reason that I do the loop over K is that if you try 10^log10(a+eps(a)) you will get back exactly the same as a2. This means that log10(a) is the same as log10(a+eps(a)) in this case, but log10(a+2*eps(a)) is different.

a+eps(a) is the next representable number after a; a+2*eps(a) is the representable number after that.

It is important for my calculation, that it ends up with exactly the same value.

You will need to switch to using the symbolic toolbox.

digits(16)
EPS = 10^(-digits());
a = sym(14.25);
a2 = 10^(log10(a));
for K = 1 : 100
    a3 = 10^(log10(a+K*EPS));
    if a3 ~= a2; break ;end
end
fprintf('a2 = %10.20f\n', double(a2));

a2 = 14.25000000000000000000

fprintf('a3 = %10.20f after %d*eps\n', double(a3), K);

a3 = 14.25000000000000000000 after 1*eps

disp(abs(a2-a))

fprintf('a2 difference = %10.20f\n', double(abs(a2 - a)));

a2 difference = -0.00000000000000000000

fprintf('a3 difference = %10.20f\n', double(abs(a3 - a)));

a3 difference = 0.00000000000000010000

[Juditha Schmidt]

Great! Thanks for all the explanations and examples. I will implement it this way! Thanks a lot!

[Paul]

@Walter Roberson

a+eps(a) is the next representable number after a; a+2*eps(a) is the representable number after that.

For ever postive value a, Isn't the 2nd representable number after a: a + eps(a) + eps(a + eps(a)) ?

[Walter Roberson]

No. eps(a) generally expresses the magnitude of the last bit position in the double precision representation. (I say generally because when you get very close to zero, there are denormalized numbers.)

format long g
a = rand() * 10^randi([-300 300])

a = 
     1.08196220699052e+146

num2hex(a)

ans = '5e4154554262a237'

nexta = typecast(typecast(a, 'uint64') + uint64(1),'double')

nexta = 
     1.08196220699052e+146

num2hex(nexta)

ans = '5e4154554262a238'

nexta - a

ans = 
     2.21813575529665e+130

eps(a)

ans = 
     2.21813575529665e+130

(nexta - a) - eps(a)

ans = 
     0