class 08 · hour eight of fifty
Why does the
order matter?
Nothing today overflows, nothing is subnormal, no NaN turns up. Every number stays in the well behaved middle of the range, and the answers are still wrong.
Eleven programs. By the end they will predict an error before measuring it.
Three numbers, added two ways.
/* 01-brackets.c */
#include <stdio.h>
int main(void)
{
float a = 1e8f;
float b = -1e8f;
float c = 1.0f;
printf("a = %g b = %g c = %g\n", a, b, c);
printf("(a + b) + c = %g\n", (a + b) + c);
printf("a + (b + c) = %g\n", a + (b + c));
printf("a + b == b + a %s\n",
a + b == b + a ? "true" : "false");
return 0;
}
a = 1e+08 b = -1e+08 c = 1
(a + b) + c = 1
a + (b + c) = 0
a + b == b + a true
remember
- Float addition is still commutative. Swapping two operands never changes anything.
- It is not associative. Moving the brackets changes the answer.
- That single fact is why a compiler is not allowed to rearrange your arithmetic.
Ask which one is right. Neither. Both are the correctly rounded result of the operations as written.
To explain it we need one number: how far apart floats are.
/* 02-gaps.c */
#include <stdio.h>
#include <math.h>
int main(void)
{
float x[] = { 1.0f, 100.0f, 10000.0f,
1000000.0f, 100000000.0f };
printf("at this size the next float is this far away\n");
for (int i = 0; i < 5; i++)
printf("%-15g %g\n", x[i],
nextafterf(x[i], 2 * x[i]) - x[i]);
return 0;
}
at this size the next float is this far away
1 1.19209e-07
100 7.62939e-06
10000 0.000976562
1e+06 0.0625
1e+08 8
remember
- The gap is always about the number divided by 8 million. It grows as the number grows.
- Near 100 million, consecutive floats are 8 apart. There is no float between 100000000 and 100000008.
- So
1e8 + 1 has nowhere to land, and comes straight back to 1e8. That is why the brackets mattered.
Go back to program 1 with this in hand. b + c is 1e8 plus one, which rounds back, so a + that is zero.
So how small does an addend have to be before it does nothing?
/* 03-swallowed.c */
#include <stdio.h>
int main(void)
{
float x[] = { 1.0f, 100.0f, 10000.0f,
1000000.0f, 100000000.0f };
printf("running total "
"anything smaller than this adds nothing\n");
for (int i = 0; i < 5; i++) {
float y = 1.0f;
while (x[i] + y != x[i]) y = y / 2.0f;
printf("%-15g %g\n", x[i], y * 2);
}
return 0;
}
running total anything smaller than this adds nothing
1 1.19209e-07
100 7.62939e-06
10000 0.000976562
1e+06 0.0625
1e+08 2
remember
- At 1e8 the gap is 8, but anything under 2 disappears. The threshold is a quarter of the gap here, not half, because 1e8 sits on a power of two boundary.
- The loop finds it by halving until the addition stops having an effect. No theory needed.
- Adding a small number to a large one is the operation that loses information. Keep running totals small.
Worth pausing on 1e8. Adding one to a hundred million, a million times, still leaves you at a hundred million.
Now something completely ordinary.
/* 04-million.c */
#include <stdio.h>
int main(void)
{
float sum = 0.0f;
for (int i = 0; i < 1000000; i++)
sum += 0.1f;
printf("added 0.1 a million times\n");
printf("expected 100000\n");
printf("got %.4f\n", sum);
printf("error %+.4f (%.3f%%)\n",
sum - 100000.0f,
(sum - 100000.0f) / 1000.0f);
return 0;
}
added 0.1 a million times
expected 100000
got 100958.3438
error +958.3438 (0.958%)
the question for this hour
- One rounding is at most half a gap, about one part in sixteen million.
- A million of them produced an error of one part in a hundred. Four orders of magnitude worse.
- Where did 958 come from? By the end of the hour you will work it out on paper before running the program.
Do not explain yet. Let them guess. Most say the errors piled up randomly, which is wrong and worth being wrong about.
Look at one addition, in slow motion.
/* 05-band.c */
#include <stdio.h>
#include <math.h>
int main(void)
{
float s = 40000.0f;
float gap = nextafterf(s, 2 * s) - s;
printf("running total %g\n", s);
printf("gap here %g\n", gap);
printf("0.1 is %.2f gaps\n", 0.1f / gap);
printf("so it must move %.0f gaps\n\n",
floorf(0.1f / gap + 0.5f));
for (int i = 0; i < 5; i++) {
float before = s;
s += 0.1f;
printf("%.4f + 0.1 -> %.4f added %.7f\n",
before, s, s - before);
}
return 0;
}
running total 40000
gap here 0.00390625
0.1 is 25.60 gaps
so it must move 26 gaps
40000.0000 + 0.1 -> 40000.1016 added 0.1015625
40000.1016 + 0.1 -> 40000.2031 added 0.1015625
40000.2031 + 0.1 -> 40000.3047 added 0.1015625
40000.3047 + 0.1 -> 40000.4062 added 0.1015625
40000.4062 + 0.1 -> 40000.5078 added 0.1015625
remember
- A float can only land on the grid. 0.1 is 25.6 gaps, and you cannot move 25.6 gaps, so it moves 26.
- Every addition adds 0.1015625 instead of 0.1. Too much by 0.0015625.
- The situation is identical each time, so the decision is identical each time. This is not a random walk. It is the same mistake, repeated.
This is the slide the whole class turns on. The five identical numbers in the right column are the evidence.
The direction changes as the total grows, in a pattern.
Each time the running total passes a power of two, the gap doubles, so 0.1 is worth half as many gaps as before.
When the fraction is .6 or .8 it rounds up and the total runs fast. When it is .2 or .4 it rounds down and the total runs slow.
remember
- The fractions cycle with period four: .2, .6, .8, .4, and round again.
- 0.1 in binary is 0.0001100110011... with a repeating block of four.
- The repeating block of one tenth turns up as the drift pattern of a loop.
Write 0.1 in binary on the board and circle the repeating 0011. Then point at the four boxes.
That is enough to work out the answer before running it.
/* 06-predict.c */
#include <stdio.h>
#include <math.h>
int main(void)
{
double predicted = 0;
for (int e = 10; e <= 16; e++) {
double hi = pow(2, e + 1) > 100000 ? 100000 : pow(2, e + 1);
double gap = pow(2, e - 23);
double got = floor(0.1 / gap + 0.5) * gap;
double err = got - 0.1;
double adds = (hi - pow(2, e)) / got;
predicted += err * adds;
printf("2^%-2d %-10.7f %-9.2f %-10.7f %+.7f %-8.0f %+.1f
",
e, gap, 0.1 / gap, got, err, adds, err * adds);
}
float sum = 0.0f;
for (int i = 0; i < 1000000; i++) sum += 0.1f;
printf("
predicted %+.1f
measured %+.1f
",
predicted, sum - 100000.0f);
return 0;
}
remember
- Each band takes twice as many additions as the one before, so the last two bands do almost all the damage.
- Predicted +943.7 against measured +958.3, which is 1.5% out.
- Floating point error is not luck. Given the magnitudes and the operation count, you can compute it in advance.
band gap 0.1 is rounds to error adds drift
2^10 0.0001221 819.20 0.0999756 -0.0000244 10243 -0.3
2^11 0.0002441 409.60 0.1000977 +0.0000977 20460 +2.0
2^12 0.0004883 204.80 0.1000977 +0.0000977 40920 +4.0
2^13 0.0009766 102.40 0.0996094 -0.0003906 82241 -32.1
2^14 0.0019531 51.20 0.0996094 -0.0003906 164483 -64.3
2^15 0.0039062 25.60 0.1015625 +0.0015625 322639 +504.1
2^16 0.0078125 12.80 0.1015625 +0.0015625 339338 +530.2
predicted +943.7
measured +958.3
The residual is real and explainable: the band boundaries are approximate because the sum accelerates. Say so rather than hiding it.
Same million numbers, added in groups.
/* 07-chunks.c */
#include <stdio.h>
int main(void)
{
printf("chunks sum error\n");
for (int chunks = 1; chunks <= 4096; chunks *= 8) {
int per = 1000000 / chunks;
float total = 0.0f;
for (int k = 0; k < chunks; k++) {
float part = 0.0f;
for (int i = 0; i < per; i++) part += 0.1f;
total += part;
}
printf("%6d %-12.4f %+.4f\n",
chunks, total, total - 100000.0f);
}
return 0;
}
chunks sum error
1 100958.3438 +958.3438
8 99910.3203 -89.6797
64 99985.0000 -15.0000
512 99995.0156 -4.9844
4096 99938.4922 -61.5078
remember
- Five ways of adding the same million numbers. Five different answers.
- Splitting the work keeps every running total small, so the grid stays fine and less is lost per step.
- This is what a parallel reduction does. Change the number of threads and you change the answer, without changing the program.
Note that more chunks is not simply better. 4096 is worse than 512, because combining that many partial sums has its own rounding.
Or keep track of what you dropped.
/* 08-kahan.c */
#include <stdio.h>
int main(void)
{
float naive = 0.0f;
for (int i = 0; i < 1000000; i++)
naive += 0.1f;
float sum = 0.0f, lost = 0.0f;
for (int i = 0; i < 1000000; i++) {
float y = 0.1f - lost; /* put back what we dropped */
float t = sum + y;
lost = (t - sum) - y; /* what got dropped now */
sum = t;
}
printf("naive %.4f error %+.4f\n",
naive, naive - 100000.0f);
printf("compensated %.4f error %+.4f\n",
sum, sum - 100000.0f);
return 0;
}
naive 100958.3438 error +958.3438
compensated 100000.0000 error +0.0000
remember
(t - sum) - y recovers exactly the part of y that did not fit. It is not an estimate.- That crumb is small, so it is added to the next small value rather than to the large total, where it would vanish again.
- Six extra lines and four times the additions, and the error goes to zero. This is Kahan summation.
Trace one iteration by hand at 40000. The crumb is exactly the 0.0015625 from program 5.
A different failure, with the same cause.
/* 09-cancel.c */
#include <stdio.h>
int main(void)
{
float x = 1e-8f;
float sum = 1.0f + x;
printf("x = %.10e\n", x);
printf("x == 0 %s\n",
x == 0.0f ? "true" : "false");
printf("1 + x = %.10e\n", sum);
printf("(1 + x) - 1 = %.10e\n", sum - 1.0f);
printf("gap at 1.0 = %.10e\n", 1.1920929e-07f);
return 0;
}
x = 9.9999999392e-09
x == 0 false
1 + x = 1.0000000000e+00
(1 + x) - 1 = 0.0000000000e+00
gap at 1.0 = 1.1920929e-07
remember
- x is nowhere near zero, but it is smaller than the gap at 1.0, so the addition threw it away.
- The subtraction is exact. It performed no rounding at all.
- The error was already there. Subtracting nearly equal numbers only removes the leading digits that were hiding it.
- When a result looks wrong, suspect the step before the subtraction, not the subtraction.
This is the one students misdiagnose. They blame the subtraction. The subtraction is the only honest operation on the slide.
The same arithmetic, on a system that was keeping time.
/* 10-dhahran.c */
#include <stdio.h>
int main(void)
{
/* the Patriot clock counted tenths of a second in a
fixed point register, keeping 23 bits after the point */
double stored = 0.0, bit = 0.5;
for (int i = 0; i < 23; i++) {
if (stored + bit <= 0.1) stored += bit;
bit /= 2.0;
}
double err = 0.1 - stored;
long ticks = 100L * 3600L * 10L;
double drift = err * ticks;
printf("0.1 kept in 23 bits = %.20f\n", stored);
printf("error per tick = %.4e seconds\n", err);
printf("ticks in 100 hours = %ld\n", ticks);
printf("clock is now off by = %.4f seconds\n", drift);
printf("a Scud at Mach 5 = %.0f metres\n",
drift * 1676.0);
return 0;
}
0.1 kept in 23 bits = 0.09999990463256835938
error per tick = 9.5367e-08 seconds
ticks in 100 hours = 3600000
clock is now off by = 0.3433 seconds
a Scud at Mach 5 = 575 metres
remember
- One tenth has no exact binary form, so a fixed point clock ticking in tenths drifts by a fixed amount every tick.
- Exactly the arithmetic from program 6: one error, multiplied by an operation count.
- Dhahran, 25 February 1991. The battery had been running about 100 hours. The interceptor missed. Twenty eight people were killed.
The figures match the US government report: 0.000000095 seconds per tick and 0.34 seconds after a hundred hours. This is the same program, not a retelling.
And then the part that has nothing to do with arithmetic.
The problem was already known. A modified version of the software, with the drift corrected, had been written and was on its way.
It reached Dhahran the following day.
remember
- The system was never designed to run for a hundred hours without a restart. Field use changed and the assumption did not.
- Rebooting reset the clock, so a machine that was restarted regularly never showed the fault.
- A numerical bug and a deployment failure are different problems, and this incident is both.
Slow down here. Do not editorialise, the facts are enough.
One more source of a different answer.
/* 11-reorder.c */
#include <stdio.h>
int main(void)
{
float sum = 0.0f;
for (int i = 0; i < 1000000; i++)
sum += 0.1f;
printf("%.4f\n", sum);
return 0;
}
$ gcc 11-reorder.c -o r -O0
$ ./r
100958.3438
$ gcc 11-reorder.c -o r -O2
$ ./r
100958.3438
$ gcc 11-reorder.c -o r -O2 -ffast-math
$ ./r
99759.8516
remember
- At -O2 the compiler optimises heavily and still does not touch the order, because it is not allowed to.
- -ffast-math is you telling it that it may pretend float addition is associative. It then splits the loop into parallel accumulators, which is faster and gives a different answer.
- Same source, same machine, same input. The flag is a promise you make, not an optimisation you get for free.
Run all three compiles live. The -O2 line matters: it shows restraint, not incapability.
Floats round. Rounding depends on what you are holding.
the three questions to ask
- What magnitudes is this running at, and what is the gap there?
- How many operations happen, and does the running total grow through them?
- Does anything subtract two numbers that are nearly equal?
and what to do about it
- Keep running totals small. Split the work, or sort smallest first.
- Keep the crumbs, if the accuracy is worth four times the additions.
- Do not compare floats for equality. Compare against a tolerance you chose deliberately.
None of this is a defect. It is a finite grid behaving exactly as specified, and it is predictable once you know where the grid lines are.
Send off: take program 4 and make it accurate without using Kahan and without changing the type. Sorting is one answer. There are others.