In MATLAB, are variables REALLY double-precision by default?

gnovice picture gnovice · Nov 19, 2010 · Viewed 27.2k times · Source

This question arose out of something strange that I noticed after investigating this question further...

I always understood MATLAB variables to be double-precision by default. So, if I were to do something like declare a variable with 20 digits after the decimal point:

>> num = 2.71828182845904553488;
>> class(num)  % Display the variable type
ans =
double

I would expect the last 4 digits to be ignored, since the floating-point relative accuracy is on the order of 10-16:

>> eps(num)
ans =
    4.440892098500626e-016

If I try to display the number with more than 16 digits after the decimal point (using either fprintf or sprintf), I get what I expect to see:

>> fprintf('%0.20f\n', num)
2.71828182845904550000
>> sprintf('%0.20f', num)
ans =
2.71828182845904550000

In other words, digits 17 through 20 are all 0.

But things get weird when I pass num to the variable precision arithmetic function in the Symbolic Toolbox, telling it to represent the number using 21 digits of precision:

>> vpa(num, 21)
ans =
2.71828182845904553488

WHAT?! Those last 4 digits have reappeared! Shouldn't they have been lost when the original number I entered was stored as a double-precision variable num? Since num is a double-precision variable when it is passed to vpa, how did vpa know what they were?

My best guess as to what is happening is that MATLAB internally represents num with more precision than a double since I initialized it to a number with more digits past the decimal point than a double-precision variable could handle. Is this really what is happening, or is something else going on?



BONUS: And here's an additional source of confusion if you don't already have a migraine from the above...

>> num = 2.71828182845904553488;  % Declare with 20 digits past the decimal
>> num = 2.718281828459045531;    % Re-declare with 18 digits past the decimal
>> vpa(num, 21)
ans =
2.71828182845904553488  % It's the original 20-digit number!!!

Answer

Andrew Janke picture Andrew Janke · Nov 19, 2010

They're doubles. vpa() is simply choosing to display non-significant digits beyond the floating point relative accuracy, where printf() and disp() truncate or zero them out.

You're only getting your original four digits back out because the literal you chose to initialize num with just happens to be the exact decimal expansion of a binary double value, because it was copy and pasted from the output of the expansion of an actual double value from the other question. It won't work for other nearby values, as you show in your "BONUS" addendum.

More precisely, all numeric literals in Matlab produce values of type double. They get converted to the binary double value that is nearest to the decimal value they represent. In effect, digits in a literal beyond the limit of precision of the double type are silently dropped. When you copy and paste the output of vpa to create a new variable, as the other question's poster did with the e = ... statement, you're initializing a value from a literal, instead of dealing directly with the result of a previous expression.

The differences here are just in output formatting. I think what's going on is that vpa() is taking that double precision binary double and treating it as an exact value. For a given binary mantissa-exponent value, you can calculate the decimal equivalent to arbitrarily many decimal places. If you have a limited precision ("width") in the binary value, as you do with any fixed-size data type, only so many of those decimal digits are significant. printf() and Matlab's default display handle this by truncating the output or displaying non-significant digits as 0. vpa() is ignoring the limits of precision and continuing to calculate as many decimal places as you request.

Those additional digits are bogus, in the sense that if they were replaced by other values to produce a nearby decimal value, they would all get "rounded" to the same binary double value.

Here's a way to show it. These values of x are all the same when stored in doubles, and will all be represented the same by vpa().

x = [
    2.7182818284590455348848081484902650117874145507812500
    2.7182818284590455348848081484902650117874145507819999
    2.7182818284590455348848
    2.71828182845904553488485555555555555555555555555555
    exp(1)
    ]
unique(x)

Here's another way of demonstrating it. Here are two doubles that are very close to each other.

x0 = exp(1)
x1 = x0 + eps(x0)

vpa(x0) and vpa(x1) should produce outputs that differ a lot past the 16th digit. However, you shouldn't be able to create a double value x such that vpa(x) produces a decimal representation that falls between vpa(x0) and vpa(x1).

(UPDATE: Amro points out that you can use fprintf('%bx\n', x) to display an exact representation of the underlying binary value in hex format. You can use this to confirm the literals map to the same double.)

I suspect vpa() behaves this way because it treats its inputs as exact values, and polymorphically supports other Matlab types from the Symbolic Toolbox that have more precision than doubles. Those values will need to be initialized by means other than numeric literals, which is why sym() takes a string as an input and vpa(exp(1)) differs from vpa(sym('exp(1)')).

Make sense? Sorry for the long-windedness.

(Note I don't have the Symbolic Toolbox so I can't test vpa() myself.)