How Did a Rounding Error Stop a Missile Defense?

Use the widget above, with Two ways pressed (as the battery was).

• a.Move the hours from 0 to 8. How far has the gate moved, and is the Scud still inside it?
• b.Find the hour at which the radar stops tracking the Scud.
• c.Go to 100 hours. Now press Same way. What happens to the gate, and to the clock’s error?

Where the radar looks next

A Patriot battery has a radar that sends out pulses and listens for echoes. When it finds something with the look of a Scud, it doesn’t keep scanning the whole sky. From the object’s speed and the time of its last echo, the computer works out where it should be at the next look, and checks only a stretch of the beam around that place: the range gate. If the object is found inside the gate, it is confirmed as a Scud and tracked. If the gate is in the wrong place, the radar is looking at empty sky.

So the prediction needs two numbers: the speed, and the time. A Scud flies at about Mach 5, the report says, 3,750 miles an hour, which is about 1,676 meters a second. An error of a tenth of a second in the time puts the prediction about 168 m out.

Time kept as a count

The computer kept time as a whole number of tenths of a second since it was switched on: 32, 33, 34… To use it in a calculation, it turned the count into seconds by multiplying it by 0.1. The report gives the reason the answer wasn’t exact: the Patriot’s weapons control computer was a 1970s design whose registers, the small stores where it does its arithmetic, were only 24 bits long.

You met the trouble with 0.1 in Grade 11: in binary it never ends.

$$0.1 = 0.0\,0011\,0011\,0011\,0011\ldots_2$$

The computer had to cut it off somewhere. Robert Skeel, a computer scientist at the University of Illinois, worked out exactly where. He noticed that the report’s errors were always the same share of the time, very nearly 2⁻²⁰ of it (about 0.0001%), and that this is just what you get if 0.1 is kept as a fixed-point number with its first 23 binary digits after the point, and the rest simply dropped. Fixed point means the point stays in one place, here at the start: every one of the 23 digits is a fraction (a half, a quarter, an eighth…), and there is no exponent, unlike a float. Dropping the rest, without rounding, is called chopping (or truncating). Chopping always makes a positive number smaller, so every time the computer worked out was a little too small: the clock fell behind.

An error that should have canceled

Now Skeel’s puzzle, the one the widget showed you. To predict where the Scud goes next, the computer doesn’t need the time since it was switched on. It needs the time between two echoes: the time of the later one minus the time of the earlier one. If both times are a little too small by the same share, the difference is too small by that same tiny share, about a millionth. For two echoes a second apart, that is a millionth of a second, and the Scud moves about 2 millimeters in that time.

So the clock could have been a third of a second wrong and still done no harm. Why did it matter?

“Further inquiries cleared up the mystery,” Skeel wrote. The program was written in assembly language about twenty years before. When the batteries were sent to the Gulf, the software was changed several times so that the system could cope with ballistic missiles as fast as a Scud, which it had not been designed for. One change added a new, more accurate way of turning the clock’s count into seconds. That conversion was needed in about half a dozen places in the program, but the call to the new one was not put in every place it was needed. So at some point the program subtracted a time worked out the old way from a time worked out the new way. One was too small by the whole drift; the other was nearly right. The difference carried the whole drift, 0.34 seconds after 100 hours, and the gate moved with it.

That is the Two ways button. A better conversion, put in some places and not others, turned a harmless error into a deadly one.

You will rebuild the clock in five steps and get the GAO’s own numbers out of it.

Step 1: 0.1 in binary

In Grade 11 you found the digits by hand: double the fraction; if it reaches 1, the next digit is 1 and you take the 1 away. Python’s floats are themselves only close to 0.1, so to get the true digits the program keeps the fraction exactly, as a numerator and a denominator, both whole numbers.

num, den = 1, 10          # the fraction 1/10digits = ''for i in range(32):    num = num * 2         # double what is left    if num >= den:        # it reached a whole one: the digit is 1        digits = digits + '1'        num = num - den    else:        digits = digits + '0'print('0.1 in binary is 0.' + digits + '...')

0.1 in binary is 0.00011001100110011001100110011001...

After the first 0, the four digits 0011 come round again and again, for ever.

Step 2: what the register kept

Keep the first 23 digits and drop the rest. int(text, 2) reads a string of binary digits as a whole number; dividing by 2²³ puts the point back in front of them.

digits = '00011001100110011001100110011001'   # from Step 1kept23 = int(digits[:23], 2)     # the first 23 binary digits, as a whole numbertenth = kept23 / 2 ** 23         # what the register's 0.1 was worthprint('Kept bits:', digits[:23])print('As a whole number:', kept23, 'out of', 2 ** 23)print('Worth:', f'{tenth:.15f}')print('Short of 0.1 by:', f'{0.1 - tenth:.15f}', 'seconds each tick')print('As a share of 0.1:', (0.1 - tenth) / 0.1, '=', 2 ** -20, '?')

Kept bits: 00011001100110011001100As a whole number: 838860 out of 8388608Worth: 0.099999904632568Short of 0.1 by: 0.000000095367432 seconds each tickAs a share of 0.1: 9.536743164617612e-07 = 9.5367431640625e-07 ?

Each tick is short by about 0.000000095 seconds, 95 billionths. The last line is Skeel’s clue: the share lost is 2⁻²⁰ (the two numbers differ only in their last few digits, because Python’s own 0.1 isn’t exact either). The kept number is 838,860 ÷ 8,388,608, so Python holds it exactly: it is a whole number divided by a power of 2.

Step 3: the clock over 100 hours

Now run the clock. The count after a number of hours is hours × 60 × 60 × 10; the true time is the count ÷ 10; the computer’s time is the count × the kept tenth.

tenth = 838860 / 2 ** 23          # the register's 0.1 for hours in [1, 8, 20, 48, 72, 100]:    ticks = hours * 60 * 60 * 10  # the clock counts tenths of a second    true = ticks / 10    clock = ticks * tenth         # how the program turned the count into seconds    print(f'{hours:>3} hours: {true:>9.0f} s, computed {clock:.4f} s, behind by {true - clock:.4f} s')

  1 hours:      3600 s, computed 3599.9966 s, behind by 0.0034 s  8 hours:     28800 s, computed 28799.9725 s, behind by 0.0275 s 20 hours:     72000 s, computed 71999.9313 s, behind by 0.0687 s 48 hours:    172800 s, computed 172799.8352 s, behind by 0.1648 s 72 hours:    259200 s, computed 259199.7528 s, behind by 0.2472 s100 hours:    360000 s, computed 359999.6567 s, behind by 0.3433 s

These are the numbers in the table at the back of the GAO report, Effect of Extended Run Time on Patriot Operation, to the last digit: 0.0034, 0.0275, 0.0687, 0.1648, 0.2472, 0.3433 seconds. Your program has found the report’s arithmetic from nothing but the binary digits of 0.1. It also finds a misprint: the report’s table gives the computed time after 100 hours as 359,999.6667, but 360,000 − 0.3433 is 359,999.6567, as your program says.

The drift grows in a straight line, 0.0034 seconds an hour, because every tick adds the same error: 36,000 ticks an hour × 0.000000095 seconds.

Step 4: two echoes, two ways

Now the part that mattered: the time between two echoes. old_way is the clock’s first conversion; new_way stands for the more accurate one added for the Gulf War (we don’t know exactly how it worked, so here it is exact). The echoes are a second apart, a number we made up: the report doesn’t give the time between looks. The speed is the report’s Mach 5.

TENTH = 838860 / 2 ** 23          # the register's 0.1, from Step 2SPEED = 3750 * 1609.344 / 3600    # a Scud: about Mach 5, 3,750 miles an hour, in meters a second def old_way(count):               # tenths to seconds, as the 1970s program did    return count * TENTH def new_way(count):               # the more accurate conversion added for the Gulf War    return count / 10 print(f'A Scud covers about {SPEED:,.0f} meters a second')for hours in [1, 8, 20, 100]:    first = hours * 36000         # the clock's count at the first radar pulse    second = first + 10           # the second pulse, a second later (made up)    same = old_way(second) - old_way(first)    mixed = new_way(second) - old_way(first)    print(f'{hours:>3} h: same way {same:.7f} s ({(same - 1) * SPEED:+.3f} m), '          f'two ways {mixed:.4f} s ({(mixed - 1) * SPEED:+.0f} m)')

A Scud covers about 1,676 meters a second  1 h: same way 0.9999990 s (-0.002 m), two ways 1.0034 s (+6 m)  8 h: same way 0.9999990 s (-0.002 m), two ways 1.0275 s (+46 m) 20 h: same way 0.9999990 s (-0.002 m), two ways 1.0687 s (+115 m)100 h: same way 0.9999990 s (-0.002 m), two ways 1.3433 s (+576 m)

Converted the same way, the gap between the echoes is wrong by a millionth of a second, and the prediction by 2 millimeters, after 1 hour or after 100: the drift cancels. Converted two ways, the gap carries the whole drift, and the prediction is 576 m out after 100 hours. Skeel put it as “more than half a kilometer”.

The report’s table gives 687 m, not 576 m. Its meters are about 2,000 times its seconds of error, so it must have used a speed of about 2,000 meters a second; it doesn’t say which speed, or why (perhaps a different speed of the Scud relative to the radar). Either way, the size of the miss is the drift times a speed.

Step 5: the fixes

Three ways to keep the gate on the Scud, each from the report or from Skeel’s note: restart the computer every few hours (the report says a restart, about 60 to 90 seconds, set the clock back to zero), use the old conversion everywhere, or use the new one everywhere.

TENTH = 838860 / 2 ** 23SPEED = 3750 * 1609.344 / 3600 def old_way(count):    return count * TENTH def new_way(count):    return count / 10 def gate_off(hours, fix):    first = hours * 36000    if fix == 'restart every 8 hours':        first = first % (8 * 36000)       # the clock starts again from 0    second = first + 10    if fix == 'new way everywhere':        gap = new_way(second) - new_way(first)    elif fix == 'old way everywhere':        gap = old_way(second) - old_way(first)    else:                                 # as fielded, or restarted        gap = new_way(second) - old_way(first)    return abs(gap - 1) * SPEED           # meters from where the Scud really is for fix in ['as fielded', 'restart every 8 hours', 'old way everywhere', 'new way everywhere']:    print(f'{fix:>22}: after 100 hours the gate is {gate_off(100, fix):,.3f} m off')

            as fielded: after 100 hours the gate is 575.546 m off restart every 8 hours: after 100 hours the gate is 23.022 m off    old way everywhere: after 100 hours the gate is 0.002 m off    new way everywhere: after 100 hours the gate is 0.000 m off

A restart caps the drift (at 100 hours, the clock has run 4 hours since its last restart), but only for as long as someone remembers to do it. The two real fixes are both “the same way everywhere”. Even the old, chopped conversion is fine, as long as it is the only one. The Army’s corrected software, the report says, “compensated for the inaccurate time calculation”.

Your programs found, from the binary digits of 0.1:

• a.A tiny error, repeated, becomes a big one. The register’s 0.1 was short by 0.000000095. After 3,600,000 ticks, the clock was 0.3433 seconds behind: the GAO’s own numbers, to the last digit (Step 3).
• b.The size of an error isn’t the whole story; its direction is. All of the clock’s errors went the same way, so they added up instead of canceling. Chopping always errs the same way; rounding to the nearest would have made each tick’s error four times smaller and the other way (Going further still).
• c.A consistent error can be harmless, and a partial fix can make it harmful. Every time worked out the old way was wrong by the same share, so the time between two echoes was right (Step 4). Only when an improved conversion was put in some places and not others did the drift land in the difference, and move the gate by hundreds of meters.

Programmers meet the third lesson far from missiles. A shop’s software that keeps money in rupees as floats in one place and in paisa as whole numbers in another, a phone app that reads the time from two clocks, a spreadsheet that rounds in one column and not the next: each can give two answers for the same thing, and the trouble shows only where the two meet.

Chop, or round?

Chopping 0.1 at 23 binary digits keeps 838,860 ÷ 2²³. Rounding to the nearest would keep 838,861 ÷ 2²³ instead, because the next digits, 11 0011…, are worth more than half of the last kept place.

More bits don’t always help

Look at 0.1’s digits again: 0.0 0011 0011… The 22nd and 23rd digits after the point are both 0. So keeping 21 digits, 22 or 23 gives exactly the same number, and the same drift. The 24th digit is a 1, so a 24th bit would have cut the error to three-eighths of it. Challenge ★★★ asks for the fewest bits for a given drift, and this pattern is the catch.

A float would round, not chop

The Patriot kept 0.1 in fixed point. A floating-point number rounds to the nearest number it can hold, which may be above or below. In the 8-bit float from Grade 11, store each of these and look at the rounding error’s sign.

An open question: a large program may convert the same quantity in many places, written by different people over many years. How could a computer find every place a program works out “the same thing” in different ways, without running the program? Tools that read code for such patterns exist, but no tool can find every case for every program.

Three challenges, each checked on the page: write your function in the program above the checker (press Edit), then press Check my program. The examples show what your function should give; there are more tests you can’t see.

★ The digits of any fraction

Write binary_digits(num, den, n): the first n binary digits after the point of the fraction num/den (where num is smaller than den), as a string of '0's and '1's.

def binary_digits(num, den, n):    # The first n binary digits of num/den after the point.    return '' print(binary_digits(1, 10, 8))

★★ How long can it run?

A clock counts tenths of a second and turns its count into seconds with 0.1 chopped to bits binary digits after the point. Write safe_hours(bits, limit): the most whole hours it can run before it is more than limit seconds behind.

def safe_hours(bits, limit):    # The most whole hours before the clock is more than limit seconds behind.    return 0 print(safe_hours(23, 0.34))

★★★ The fewest bits

Write bits_needed(hours, limit): the fewest binary digits of 0.1 to keep (chopped) so that a clock counting tenths is at most limit seconds behind after hours hours. The tests ask for clocks that run for years, with limits as small as a billionth of a second.

def bits_needed(hours, limit):    return 0 print(bits_needed(100, 0.35))

• Grade 11, Numbers with fractions: 0.1 in binary by hand, floating point, and why 0.1 + 0.2 isn’t 0.3.
• Grade 9, Number systems: binary place values, the start of everything here.
• Grade 12, Tests that run themselves: choosing test data, including the longest run the program will ever have.
• Why Did a Rocket Destroy Itself 37 Seconds After Launch?: the companion story, a number too big for its box, and a program reused without checking its numbers.