Why does 10^61 not equal 1e61 in MATLAB?

I noticed that 10^61==1e61 evaluates to false (0) in MATLAB. This implies that a^b and 1eX are computed differently. I am wondering how each one is computed, and why are they computed differently? In other words, why is 1eX not just computed as 10^X using the existing exponentiation function?

So far, I have figured out that 10^61-1e61==eps(1e61), meaning that 10^61 and 1e61 are adjacent floating point numbers, and that the true value is in between these two, and that 1e61 is closer than 10^61 to the true value. But I am not sure why they are different -- in other programming languages, like Python and even Octave, they are the same. But in MATLAB, they seem to be different?

Answer 1

[Roger Stafford]

Matlab uses two different functions for these computations, the exponentiation function for 10^61 and the decimal-to-binary conversion function for 1e61. The two can easily arrive at a slightly different answer due to differences in rounding in the respective functions, since in general the two functions have to deal with different types of numbers. This should not be surprising to anyone using digital computers, since rounding errors are an ever-present element in the great majority of computations.

On my own computer the two quantities you mention actually produce identical results. However if I change to 1e64 and 10^64, the results differ by a single bit out of the 53 bits of their respective significands. Why should such a minuscule difference matter? In any case it is impossible to represent either of these quantities exactly as a binary floating point number.

[dpb]

"_... it is impossible to represent either of these quantities exactly as a binary floating point number."_

?? They're integers well within the range of FP.

I'm still surprised Matlab doesn't special case it as an optimization if nothing else.

[Roger Stafford]

I disagree dpb! The integer 5^61 expressed in binary has 142 bits of which the top and bottom bits are ones, and there are only 53 significand bits available in double precision, so 10^61, which would be just the same with shifted binary point, could not possibly be represented exactly.

[Michael]

Roger, I understand that 10^61 cannot be represented exactly as a floating point number. I would like to understand how these quantities (10^61 and 1e61) are computed, and why they are computed differently. Could you elaborate more on how the decimal-to-binary conversion works for 1e61? Is this method used because it is more efficient than just the regular exponentiation function for arbitrary bases?

Thank you.

[dpb]

@Roger -- Yeah, brain cramp, sorry! Wasn't thinking clearly.

@Michael -- I don't know how ML implements the exponentiation operator, but the conversion will revert to the RTL support for the underlying compiler I'd presume.

[Walter Roberson]

You can see the processing of exponents somewhere in https://github.com/zerovm/glibc/blob/master/stdio-common/vfscanf.c . The processing was not obvious to me though.

10^61 is probably handled by power(), not by interpretation-time parsing. In turn power() probably calls C / C++'s pow(), which probably calls machine instructions. In older material I found that a then-typical implementation involved taking a logarithm of the base, multiplying by the exponent, and then taking the exponentiation of the result, and the documentation in MATLAB of how negative numbers to a power are processed would tend to suggest that path is more generally in use.

Answer 2

[dpb]

>> for i=3:40, if 10^i~=eval(['1E' num2str(i)]),disp(i),break,end,end
  29
>>

I may be wrong but I think most languages special-case the case 1^N; apparently Matlab doesn't.