Up: Testing
Previous: Interface and Implementation
Next: Fixtures

Floating Point

July 7th, 2011 Leave a comment Go to comments
slide 01 Hello, and welcome to the fifth episode of the Software Carpentry lecture on testing. In this episode, we’ll take a look at how computers represent real numbers, and some of the traps waiting for anyone who uses them.
slide 02 Let’s go back to the fields of Saskatchewan. Instead of working directly with pixels in photographs…
slide 03 …we would probably actually work with some kind of geographic information system that would store each field’s coordinates and ground cover.
slide 04 Those coordinates would probably be recorded as latitude and longitude.
slide 05 Which the computer would store as floating point numbers rather than integers.
slide 06 And that’s where our problems start. To understand, we have to step back and look at how computers store numbers.
slide 07 Finding a good representation for floating point numbers is hard.
slide 07 The root of the problem is that we cannot represent an infinite number of real values with a finite set of bit patterns. And unlike integers, no matter what values we do represent, there will be an infinite number of values between each of them that we can’t.
slide 09 The explanation that follows is simplified—possibly over-simplified—to keep it manageable.
slide 10 Please, if you’re doing any calculation on a computer at all, take half an hour to read Goldberg’s excellent paper, “What Every Computer Scientist Should Know About Floating-Point Arithmetic”.
slide 11 So: floating point numbers are usually represented using sign, magnitude, and an exponent.
slide 12 In a 32-bit word, the IEEE 754 standard calls for 1 bit of sign, 23 bits for the magnitude (or mantissa), and 8 bits for the exponent.
slide 13 To illustrate the problems with floating point, we’ll use a much dumber representation.
slide 14 We’ll only worry about positive values without fractional parts.
slide 15 And to make it all fit on our slides, we’ll only use 5 bits: 3 for the magnitude, and 2 for the exponent.
slide 16 Here are the possible values (in decimal) that we can represent this way. Each one is the mantissa times two to the exponent.
slide 17 For example, the decimal values 48 is binary 110 times 2 to the binary 11 power.
slide 18 Which is 6 times 2 to the third…
slide 19 …or 6 times 8.
slide 20 Of course, real floating point representations like the IEEE 754 standard don’t have all the redundancy that you see in this table. Again, if you’d like to know how they get rid of it, and what implications this has, please read Goldberg’s paper.
slide 21 Here’s a clearer view of the values our scheme can represent.
slide 22 The first thing you should notice is that there are a lot of values we can’t store. We can do 8 and 10, for example, but not 9.
slide 23 This is exactly like the problems hand calculators have with fractions like 1/3: in decimal, we have to round that to 0.3333 or 0.3334.
slide 24 But if this scheme has no representation for 9…
slide 25 …then 8+1 must be stored as either 8 or 10.
slide 26 Which raises an interesting question: if 8+1 is 8, what is 8+1+1?
slide 27 If we add from the left, 8+1 is 8, plus another 1 is 8 again.
slide 28 If we add from the right, though, 1+1 is 2, and 2+8 is 10. Changing the order of operations can make the difference between right and wrong.
slide 29 This is the sort of thing that numerical analysts spend their time on. In this case, if we sort the values we’re adding, then add from smallest to largest, it gives us the best chance of getting the best possible answer. In other situations, like inverting a matrix, the rules are much, much more complicated—so complicated that we’ll skip past them completely.
slide 30 Instead, here’s another observation about our uneven number line:
slide 31 The spacing between the values we can represent is uneven.
slide 32 However, the relative spacing between each set of values stays the same, i.e., the first group is separated by 1, then the separation becomes 2, then 4, then 8, so that the ratio of the spacing to the values stays roughly constant.
slide 33 This happens because we’re multiplying the same fixed set of mantissas by ever-larger exponents, and it points us at a couple of useful definitions.
slide 34 The absolute error in some approximation to a value is simply the absolute value of the difference between the two.
slide 35 The relative error, on the other hand, is the ratio of the absolute error to the value we’re approximating.
slide 37 This example might help. If we’re off by 1 in approximating 8+1 and 56+1, that’s the same absolute error, but the relative error is much larger in the first case than in the second.
slide 38 When we’re thinking about floating point numbers, relative error is almost always more useful than absolute—after all, it makes little sense to say that we’re off by a hundredth when the value in question is a billionth.
slide 39 What does this have to do with testing?
slide 40 Well, let’s have a look at a little program
slide 41 The loop runs over the integers from 1 to 9 inclusive.
slide 42 Using those values, we create the numbers 0.9, 0.09, 0.009, and so on, and put them in the list vals. (Remember, the double star means “raised to the power” in Python.)
slide 43 We then calculate the sum of those numbers. Clearly, this should be 0.9, 0.99, 0.999, and so on.
slide 44 But is it?
slide 45 Let’s calculate the same value a different way, by subtracting .1 from 1, then subtracting .01 from 1, and so on.
slide 46 This should create exactly the same sequence of numbers.
slide 47 Let’s check by finding the difference and printing it out.
slide 48 Here’s our answer. The first column is the loop index; the second, what we got by totaling up 0.9, 0.99, and so on, and the third is the difference between that number and what we got by subtracting from 1.0.
slide 49 The first thing you should notice is that the very first value contributing to our sum is already slightly off. Even with 23 bits for a mantissa, we cannot exactly represent 0.9 in base 2, any more than we can exactly represent 1/3 in base 10. Doubling the size of the mantissa would reduce the error, but we can’t ever eliminate it.
slide 50 The good news is, 9 times 10-1 and 1 minus 0.1 are exactly the same—it might not be precisely right, but at least it’s consistent.
slide 51 The same cannot be said for some of the later values, though: at some points, there actually is a small difference between what we get from our two different formulas.
slide 52 And if you weren’t confused enough, sometimes the accumulation of addition errors actually makes the result more accurate once again.
slide 53 So what does this have to do with testing?
slide 54 Well, if the function you’re testing uses floating point numbers, what do you compare its result to?
slide 55 If we compared the sum of the first few numbers in vals to what it’s supposed to be, the answer could be False.
slide 56 Even if we’re initializing the list with the right values…
slide 57 …and calculating the sum correctly.
slide 58 This is a hard problem—a very hard problem.
slide 59 No one has a good generic answer, because the root of our problem is that we’re using approximations, and each approximation has to be judged on its own merits.
slide 60 So what can you do to test your programs?
slide 61 The first rule is, compare what you get to analytic solutions whenever you can. For example, if you’re looking at the behavior of drops of liquid helium, start by checking your program’s output on a stationary spherical drop in zero gravity. You should be able to calculate the right answer in that case, and if your program doesn’t work for it, it probably won’t work for anything else.
slide 62 The second rule is to compare more complex versions of your code to simpler ones. If you’re about to replace a simple algorithm for calculating heat transfer with one that’s more complex, but hopefully faster, don’t throw the old code away: use its output as a check on the correctness of the new code. And if you bump into someone at a conference who has a program that can calculate some of the same results as yours, swap data sets: it’ll help you both.
slide 63 The third rule is, never use == (or !=) on floating point numbers, because two numbers calculated in different ways will probably not have exactly the same bits.
slide 64 (It’s OK to use less than, greater than or equals, and other orderings, though—we’ll explore why in the exercises.)
slide 65 If you want to compare floating point numbers, use a function like nose.assert_almost_equals. As its documentation says:
slide 66 This checks whether two numbers are within some tolerance of each other. If they are, it declares them equal.
slide 67 And if you thought to wonder whether that tolerance is absolute or relative, you’re ready to start testing floating-point code.
slide 68

  1. Jonathan Dursi
    August 25th, 2010 at 15:31 | #1
    Reply | Quote

    This is great and important stuff for students to be learning, and testing is a good place for it. This covers the important points well.

    Some thoughts –

    If you’re looking for something to drop, I’ve not found that knowing the twos-complement representation of negative integers has improved my life in any significant way.

    The points about truncating the number line for integers and discretizing the number line for floats (and thus introducing the possibility of mathematically equal expressions not being numerically equal) can be made before walking through the entire floating point number representation, and indeed can motivate it.

    The bit about “No one has a good answer yet for the problem of testing scientific programs that use floating point arithmetic” is a little overwrought. There’s no automatic way to generate tests for programs that use only integer arithmetic, either; you need to understand your problem to generate good tests. Period.

    This gets complicated when you introduce any sort of approximation, but floating point is only one way that can happen. If you use (say) heuristics to find near-optimal cases in large search spaces (I’m looking at you, bioinformatics!) than changing parameters in your search strategy — or getting the same data in a different order — can also change your answers, even with purely integer arithmetic! And once again, you need to know whether the answer is “correct enough”.

    If you’re looking for something to add, noting the special floating-point cases of +/-Inf and NaN is possibly useful.

  2. Rachel Slaybaugh
    December 5th, 2010 at 17:51 | #2
    Reply | Quote

    Great concepts to introduce. There is a type in the small number line example table, the last entry in the left-most column should be 111.

  3. Mark Betnel
    February 1st, 2011 at 22:35 | #3
    Reply | Quote

    Eye-opening!

  4. Mark Betnel
    February 3rd, 2011 at 16:45 | #4
    Reply | Quote

    I think that should be nose.test.assert_almost_equals() in the last couple of slides

  5. Mark Betnel
    February 3rd, 2011 at 16:48 | #5
    Reply | Quote

    oops, nose.tools.assert_almost_equals()

  6. Rivo Randrianarivony
    June 11th, 2011 at 17:51 | #6
    Reply | Quote

    Great episode. But isn’t there a built-in (or from SciPy/NumPy) function to check floating-point “equality” like R’s all.equal() ??