A report about a possible bug in format bank and a visit to a local hardware store made me realize that doing decimal arithmetic with binary floating point numbers is like tightening a European bolt with an American socket wrench.
It's mid-April and so those of us who file United States income taxes have that chore to do. Years ago, I used MATLAB to help with my taxes. I had a program named form1040.m that had one statement for each line on the tax form. I just had to enter my income and deductions. Then MATLAB would do all the arithmetic.
If we're really meticulous about our financial records, we keep track of things to the nearest penny. So line 28 of form1040 might have been something like
interest = 48.35
interest = 48.3500
I didn't like those trailing zeros in the output. So I introduced
into MATLAB. Now
interest = 48.35
is printed with just two digits after the decimal point.
format bank has turned out to be useful more broadly and is still in today's MATLAB.
We recently had a user ask about the rounding rule employed by format bank. What if a value fails halfway between two possible outputs? Which is chosen and why?
Here's the example that prompted the user's query. Start with
so we can see four decimal places.
x = (5:10:45)'/1000 y = 1+x z = 2+x
x = 0.0050 0.0150 0.0250 0.0350 0.0450 y = 1.0050 1.0150 1.0250 1.0350 1.0450 z = 2.0050 2.0150 2.0250 2.0350 2.0450
These values appear to fall halfway between pairs of two-digit decimal fractions. Let's see how the ties are broken.
format bank x y z
x = 0.01 0.01 0.03 0.04 0.04 y = 1.00 1.01 1.02 1.03 1.04 z = 2.00 2.02 2.02 2.04 2.04
Look at the last digits. What mysterious hand is at work here? Three of the x's, x(1), x(3), and x(4), have been rounded up. None of the y's have been rounded up. Two of the z's, z(2) and z(4), have been rounded up.
Email circulated internally at MathWorks for a few days after this was reported suggesting various explanations. Is it a bug in format bank? In the I/O library? Does it depend upon which compiler was used to build MATLAB? Do we see the same behavior on the PC and the Mac? Has it always been this way? These are the usual questions that we ask ourselves when we see curious behavior.
Do you know what's happening?
Well, none of the suspects I just mentioned is the culprit. The fact is none of the numbers fall exactly on a midpoint. Think binary, not decimal. A value like 0.005 expressed as a decimal fraction cannot be represented exactly as a binary floating point number. Decimal fractions fall in between the binary fractions, and rarely fall precisely half way.
To understand what is going on, the Symbolic Toolbox function
is your friend. The 'e' flag is provided for this purpose. The help entry says
'e' stands for 'estimate error'. The 'r' form is supplemented by a term involving the variable 'eps' which estimates the difference between the theoretical rational expression and its actual floating point value. For example, sym(3*pi/4,'e') is 3*pi/4-103*eps/249.
To see how this works for the situation encountered here.
symx = sym(x,'e') symy = sym(y,'e') symz = sym(z,'e')
symx = eps/2133 + 1/200 3/200 - eps/400 eps/160 + 1/40 (3*eps)/200 + 7/200 9/200 - (3*eps)/400 symy = 201/200 - (12*eps)/25 203/200 - (11*eps)/25 41/40 - (2*eps)/5 207/200 - (9*eps)/25 209/200 - (8*eps)/25 symz = 401/200 - (12*eps)/25 (14*eps)/25 + 403/200 81/40 - (2*eps)/5 (16*eps)/25 + 407/200 409/200 - (8*eps)/25
The output is not as clear as I would like to see it, but I can pick off the sign of the error terms and find
x + - + + - y - - - - - z - + - + -
Or, I can compute the signs with
format short sigx = sign(double(symx - x)) sigy = sign(double(symy - y)) sigz = sign(double(symz - z))
sigx = 1 -1 1 1 -1 sigy = -1 -1 -1 -1 -1 sigz = -1 1 -1 1 -1
The sign of the error term tells us whether the floating point numbers stored in x, y and z are larger or smaller than the anticipated decimal fractions. After this initial input, there is essentially no more roundoff error. format bank will round up or down from the decimal value accordingly. Again, it is just doing its job on the input it is given.
Photo credit: http://toolguyd.com/
A high-end socket wrench set includes both metric (left) and fractional inch (right) sizes. Again, think decimal and binary. Metric sizes of nuts and bolts and wrenches are specified in decimal fractions, while the denominators in fractional inch sizes are powers of two.
Conversion charts between the metric and fractional inch standards abound on the internet. Here is a link to one of them: Wrench Conversion Chart.
But we can easily compute our own conversion chart. And in the process compute delta, the relative error made when the closest binary wrench is used on a decimal bolt.
Inch Metric delta 1/64 0.016 1/32 0.031 3/64 0.047 1mm 0.039 -0.191 1/16 0.063 5/64 0.078 2mm 0.079 0.008 3/32 0.094 7/64 0.109 1/8 0.125 3mm 0.118 -0.058 9/64 0.141 5/32 0.156 4mm 0.157 0.008 11/64 0.172 3/16 0.188 13/64 0.203 5mm 0.197 -0.032 7/32 0.219 15/64 0.234 6mm 0.236 0.008 1/4 0.250 9/32 0.281 7mm 0.276 -0.021 5/16 0.313 8mm 0.315 0.008 11/32 0.344 9mm 0.354 0.030 3/8 0.375 13/32 0.406 10mm 0.394 -0.032 7/16 0.438 15/32 0.469 12mm 0.472 0.008 1/2 0.500 9/16 0.563 14mm 0.551 -0.021 5/8 0.625 16mm 0.630 0.008 11/16 0.688 18mm 0.709 0.030 3/4 0.750 13/16 0.813 20mm 0.787 -0.032 7/8 0.875 22mm 0.866 -0.010 15/16 0.938 24mm 0.945 0.008 1 1.000
Let's plot those relative errors. Except for the small sizes, where this set doesn't have enough wrenches, the relative error is only a few percent. But that's still enough to produce a damaging fit on a tight nut.
bar(k,d) axis([0 25 -.05 .05]) xlabel('millimeters') ylabel('delta')
You might notice that my conversion chart, like all such charts, and like the wrenches themselves, has a little bit of floating point character. The spacing of the entries is not uniform. The spacing between the binary values is 1/64, then 1/32, then 1/16. The spacing of the metric values is 1mm at the top of the chart and 2mm later on.
Suppose we want to tighten a 10mm nut and all we have are binary wrenches. The diameter of the nut in inches is
meter = 39.370079; d = 10*meter/1000
d = 0.3937
Consulting our chart, we see that a 13/32 wrench is the best fit, but it's a little too large. 10mm lies between these two binary values.
[floor(32*d) ceil(32*d)] b1 = 12/32; b2 = 13/32; [b1 d b2]'
ans = 12 13 ans = 0.3750 0.3937 0.4063
The fraction of the interval is
f = (d - b1)/(b2 - b1)
f = 0.5984
10mm is about 60% of the way from 12/32 inches to 13/32 inches.
Now let's turn to floating point numbers. What happens when execute this MATLAB statement?
x = 1/10
x = 0.1000
One is divided by ten and the closest floating point number is stored in x. The same value is produced by the statement
x = 0.1
x = 0.1000
The resulting x lies in the interval between 1/16 and 1/8. The floating point numbers in this interval are uniformly spaced with a separation of
e = eps(1/16)
e = 1.3878e-17
This is 2^-56. Let
e = sym(e)
e = 1/72057594037927936
This value e plays the role that 1/64 plays for my wrenches.
The result of x = 1/10 lies between these two binary fractions.
b1 = floor(1/(10*e))*e b2 = ceil(1/(10*e))*e
b1 = 7205759403792793/72057594037927936 b2 = 3602879701896397/36028797018963968
The tiny interval of length e is
c = [b1 1/10 b2]'
c = 7205759403792793/72057594037927936 1/10 3602879701896397/36028797018963968
ans = 0.099999999999999991673327315311326 0.1 0.10000000000000000555111512312578
Where in this interval does 1/10 lie?
f = double((1/10 - b1)/e)
f = 0.6000
So 1/10 is about 60% of the way from b1 to b2 and so is closer to b2 than to b1. The two statements
x = 1/10 x = double(b2)
x = 0.1000 x = 0.1000
store exactly the same value in x.
The 13/32 wrench is the closest tool in the binary collection to the 10mm nut.
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.