MATLAB: Do I see strange results when I use the FPRINTF command with the %g format on a SunOS

fprintflinuxMATLABunix

Why do I see strange results when I use the FPRINTF command with the %g format on a SunOS?
The following command displays strange results as shown below:
>> fprintf('%.30g',1.0003e-28)
1.00030000000000000769130446023e-28
The expected output is
1.0003e-28

Best Answer

The answer that you see is the result of the following 3 facts:
1. Floating point numbers use a binary representation.
2. 1.0003e-28 does not have an exact binary representation in finite precision.
If the IEEE format had a longer mantissa then we would have:
3a1fb365d1ac67e8ee6fa88b480dd8...
But MATLAB gives the following result for a standard IEEE double which is the correctly rounded result for the stored precision:
3a1fb365d1ac67e9
When you convert that rounded value back into decimal you get what is produced.
3. Never print more than 16 decimal digits. Any more digits are often times not very useful.
Consider the following MATLAB script and its output:
format hex
a = 1.0003e-28;
am = a - .5 * eps * 1e-28
a
ap = a + .5 * eps * 1e-28
%
















fprintf('%.30g\n',am);
fprintf('%.30g\n',a);
fprintf('%.30g\n',ap);
%
%
b = zeros(32,1);
b(28) = 1;
b(32) = 3;
%
d(b,200)
am =
3a1fb365d1ac67e8 (This is smaller by 1 significant bit)
a =
3a1fb365d1ac67e9 (MATLAB converted value)
ap =
3a1fb365d1ac67ea (This is larger by 1 significant bit)
1.00029999999999989558742731424e-28 (smaller by 1 binary bit)
1.00030000000000000769130446023e-28 (value)
1.00030000000000011979518160621e-28 (bigger by 1 binary bit)
1234567890123456 <- digits of precision (stop at 16)
ans =
3a1fb365d1ac67e8ee6fa88b480dd8 (multiple precision 'IEEE number')
The routine d is:
function v = d(u, digits)
%
% abstract: d takes a positive decimal fraction less than 1
% and converts it to a binary fraction of length
% digits. It then converts it to a IEEE double
% precision number but possibly with a longer
% mantissa.
%
% u is a vector of decimal digits in the range [0-9]
% with the decimal point to the 'left' of u(1).
%
% digits is the length in binary digits to save.
%
% v is a string of hexidecimal digits of the IEEE
% double number, but with possibly longer mantiassa
%
% Algorithm: Knuth, Vol. 2, p. 319, Method 2a
%
% By: Martin Knapp-Cordes
% MathWorks, Inc.
%
% Date: Apr 2001 - 09
%
%-------------------------------------------------------------
%
u2 = zeros(size(u,1)+1,1);
u2(2:end) = u;
b = zeros(digits,1);
for i=1:digits
u2 = mult2(u2);
b(i) = u2(1);
u2(1) = 0;
end
%
% Always produce an least an IEEE double.
%
o = find(b == 1);
base = o(1);
s = size(b,1);
v = dec2hex(1023-base,3);
n = max(ceil((s - base) / 4), 13);
v2 = zeros(n * 4,1);
v2(1:(s - base)) = b((base+1):end);
for i=1:4:n*4
v = [ v dec2hex(bin2dec(sprintf('%d',v2(i:i+3))))];
end
v = lower(v);
function u2 = mult2(u)
%
% Multiply by 2 with carry
%
u2 = 2 .* u;
c = (u2>=10);
u2(1:end-1) = u2(1:end-1) + c(2:end) - c(1:end-1) .* 10;
u2(end) = u2(end) - c(end) .* 10;