# Bank Format and Metric Socket Wrenches

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.

### Contents

#### Tax Time

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

```
format bank
```

into MATLAB. Now

interest

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.

#### Breaking Ties

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

```
format short
```

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?

#### None of the above.

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.

#### Symbolic Toolbox

To understand what is going on, the Symbolic Toolbox function

sym(x,'e')

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.

#### Socket Wrenches

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.

make_chart

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.

#### 10 millimeters

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.

#### One-tenth

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

In decimal

vpa(c)

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.

**カテゴリ:**- Fun,
- Numerical Analysis,
- Precision

## コメント

コメントを残すには、ここ をクリックして MathWorks アカウントにサインインするか新しい MathWorks アカウントを作成します。