Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 1 | .. _tut-fp-issues: |
| 2 | |
| 3 | ************************************************** |
| 4 | Floating Point Arithmetic: Issues and Limitations |
| 5 | ************************************************** |
| 6 | |
| 7 | .. sectionauthor:: Tim Peters <tim_one@users.sourceforge.net> |
| 8 | |
| 9 | |
| 10 | Floating-point numbers are represented in computer hardware as base 2 (binary) |
| 11 | fractions. For example, the decimal fraction :: |
| 12 | |
| 13 | 0.125 |
| 14 | |
| 15 | has value 1/10 + 2/100 + 5/1000, and in the same way the binary fraction :: |
| 16 | |
| 17 | 0.001 |
| 18 | |
| 19 | has value 0/2 + 0/4 + 1/8. These two fractions have identical values, the only |
| 20 | real difference being that the first is written in base 10 fractional notation, |
| 21 | and the second in base 2. |
| 22 | |
| 23 | Unfortunately, most decimal fractions cannot be represented exactly as binary |
| 24 | fractions. A consequence is that, in general, the decimal floating-point |
| 25 | numbers you enter are only approximated by the binary floating-point numbers |
| 26 | actually stored in the machine. |
| 27 | |
| 28 | The problem is easier to understand at first in base 10. Consider the fraction |
| 29 | 1/3. You can approximate that as a base 10 fraction:: |
| 30 | |
| 31 | 0.3 |
| 32 | |
| 33 | or, better, :: |
| 34 | |
| 35 | 0.33 |
| 36 | |
| 37 | or, better, :: |
| 38 | |
| 39 | 0.333 |
| 40 | |
| 41 | and so on. No matter how many digits you're willing to write down, the result |
| 42 | will never be exactly 1/3, but will be an increasingly better approximation of |
| 43 | 1/3. |
| 44 | |
| 45 | In the same way, no matter how many base 2 digits you're willing to use, the |
| 46 | decimal value 0.1 cannot be represented exactly as a base 2 fraction. In base |
| 47 | 2, 1/10 is the infinitely repeating fraction :: |
| 48 | |
| 49 | 0.0001100110011001100110011001100110011001100110011... |
| 50 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 51 | Stop at any finite number of bits, and you get an approximation. On most |
| 52 | machines today, floats are approximated using a binary fraction with |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 53 | the numerator using the first 53 bits starting with the most significant bit and |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 54 | with the denominator as a power of two. In the case of 1/10, the binary fraction |
| 55 | is ``3602879701896397 / 2 ** 55`` which is close to but not exactly |
| 56 | equal to the true value of 1/10. |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 57 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 58 | Many users are not aware of the approximation because of the way values are |
| 59 | displayed. Python only prints a decimal approximation to the true decimal |
| 60 | value of the binary approximation stored by the machine. On most machines, if |
| 61 | Python were to print the true decimal value of the binary approximation stored |
| 62 | for 0.1, it would have to display :: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 63 | |
| 64 | >>> 0.1 |
| 65 | 0.1000000000000000055511151231257827021181583404541015625 |
| 66 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 67 | That is more digits than most people find useful, so Python keeps the number |
| 68 | of digits manageable by displaying a rounded value instead :: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 69 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 70 | >>> 1 / 10 |
| 71 | 0.1 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 72 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 73 | Just remember, even though the printed result looks like the exact value |
| 74 | of 1/10, the actual stored value is the nearest representable binary fraction. |
| 75 | |
| 76 | Interestingly, there are many different decimal numbers that share the same |
| 77 | nearest approximate binary fraction. For example, the numbers ``0.1`` and |
| 78 | ``0.10000000000000001`` and |
| 79 | ``0.1000000000000000055511151231257827021181583404541015625`` are all |
| 80 | approximated by ``3602879701896397 / 2 ** 55``. Since all of these decimal |
Georg Brandl | b58f46f | 2009-04-24 19:06:29 +0000 | [diff] [blame] | 81 | values share the same approximation, any one of them could be displayed |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 82 | while still preserving the invariant ``eval(repr(x)) == x``. |
| 83 | |
| 84 | Historically, the Python prompt and built-in :func:`repr` function would chose |
Georg Brandl | 7d1e880 | 2009-08-13 08:47:18 +0000 | [diff] [blame^] | 85 | the one with 17 significant digits, ``0.10000000000000001``. Starting with |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 86 | Python 3.1, Python (on most systems) is now able to choose the shortest of |
| 87 | these and simply display ``0.1``. |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 88 | |
| 89 | Note that this is in the very nature of binary floating-point: this is not a bug |
| 90 | in Python, and it is not a bug in your code either. You'll see the same kind of |
| 91 | thing in all languages that support your hardware's floating-point arithmetic |
| 92 | (although some languages may not *display* the difference by default, or in all |
| 93 | output modes). |
| 94 | |
Mark Dickinson | 934896d | 2009-02-21 20:59:32 +0000 | [diff] [blame] | 95 | Python's built-in :func:`str` function produces only 12 significant digits, and |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 96 | you may wish to use that instead. It's unusual for ``eval(str(x))`` to |
| 97 | reproduce *x*, but the output may be more pleasant to look at:: |
| 98 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 99 | >>> str(math.pi) |
| 100 | '3.14159265359' |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 101 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 102 | >>> repr(math.pi) |
| 103 | '3.141592653589793' |
| 104 | |
| 105 | >>> format(math.pi, '.2f') |
| 106 | '3.14' |
| 107 | |
| 108 | It's important to realize that this is, in a real sense, an illusion: you're |
| 109 | simply rounding the *display* of the true machine value. |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 110 | |
Raymond Hettinger | f0320c7 | 2009-04-26 20:10:50 +0000 | [diff] [blame] | 111 | One illusion may beget another. For example, since 0.1 is not exactly 1/10, |
Raymond Hettinger | 4af3629 | 2009-04-26 21:37:46 +0000 | [diff] [blame] | 112 | summing three values of 0.1 may not yield exactly 0.3, either:: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 113 | |
Raymond Hettinger | 4af3629 | 2009-04-26 21:37:46 +0000 | [diff] [blame] | 114 | >>> .1 + .1 + .1 == .3 |
| 115 | False |
| 116 | |
| 117 | Also, since the 0.1 cannot get any closer to the exact value of 1/10 and |
| 118 | 0.3 cannot get any closer to the exact value of 3/10, then pre-rounding with |
| 119 | :func:`round` function cannot help:: |
| 120 | |
| 121 | >>> round(.1, 1) + round(.1, 1) + round(.1, 1) == round(.3, 1) |
| 122 | False |
| 123 | |
| 124 | Though the numbers cannot be made closer to their intended exact values, |
| 125 | the :func:`round` function can be useful for post-rounding so that results |
Georg Brandl | 7d1e880 | 2009-08-13 08:47:18 +0000 | [diff] [blame^] | 126 | with inexact values become comparable to one another:: |
Raymond Hettinger | 4af3629 | 2009-04-26 21:37:46 +0000 | [diff] [blame] | 127 | |
Georg Brandl | 7d1e880 | 2009-08-13 08:47:18 +0000 | [diff] [blame^] | 128 | >>> round(.1 + .1 + .1, 10) == round(.3, 10) |
Raymond Hettinger | 4af3629 | 2009-04-26 21:37:46 +0000 | [diff] [blame] | 129 | True |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 130 | |
| 131 | Binary floating-point arithmetic holds many surprises like this. The problem |
| 132 | with "0.1" is explained in precise detail below, in the "Representation Error" |
| 133 | section. See `The Perils of Floating Point <http://www.lahey.com/float.htm>`_ |
| 134 | for a more complete account of other common surprises. |
| 135 | |
| 136 | As that says near the end, "there are no easy answers." Still, don't be unduly |
| 137 | wary of floating-point! The errors in Python float operations are inherited |
| 138 | from the floating-point hardware, and on most machines are on the order of no |
| 139 | more than 1 part in 2\*\*53 per operation. That's more than adequate for most |
Georg Brandl | 7d1e880 | 2009-08-13 08:47:18 +0000 | [diff] [blame^] | 140 | tasks, but you do need to keep in mind that it's not decimal arithmetic and |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 141 | that every float operation can suffer a new rounding error. |
| 142 | |
| 143 | While pathological cases do exist, for most casual use of floating-point |
| 144 | arithmetic you'll see the result you expect in the end if you simply round the |
| 145 | display of your final results to the number of decimal digits you expect. |
Benjamin Peterson | e6f0063 | 2008-05-26 01:03:56 +0000 | [diff] [blame] | 146 | :func:`str` usually suffices, and for finer control see the :meth:`str.format` |
| 147 | method's format specifiers in :ref:`formatstrings`. |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 148 | |
Raymond Hettinger | eba99df | 2008-10-05 17:57:52 +0000 | [diff] [blame] | 149 | For use cases which require exact decimal representation, try using the |
| 150 | :mod:`decimal` module which implements decimal arithmetic suitable for |
| 151 | accounting applications and high-precision applications. |
| 152 | |
| 153 | Another form of exact arithmetic is supported by the :mod:`fractions` module |
| 154 | which implements arithmetic based on rational numbers (so the numbers like |
| 155 | 1/3 can be represented exactly). |
| 156 | |
Guido van Rossum | 0616b79 | 2007-08-31 03:25:11 +0000 | [diff] [blame] | 157 | If you are a heavy user of floating point operations you should take a look |
| 158 | at the Numerical Python package and many other packages for mathematical and |
| 159 | statistical operations supplied by the SciPy project. See <http://scipy.org>. |
Raymond Hettinger | 9fce0ba | 2008-10-05 16:46:29 +0000 | [diff] [blame] | 160 | |
| 161 | Python provides tools that may help on those rare occasions when you really |
| 162 | *do* want to know the exact value of a float. The |
| 163 | :meth:`float.as_integer_ratio` method expresses the value of a float as a |
| 164 | fraction:: |
| 165 | |
| 166 | >>> x = 3.14159 |
| 167 | >>> x.as_integer_ratio() |
Georg Brandl | 7d1e880 | 2009-08-13 08:47:18 +0000 | [diff] [blame^] | 168 | (3537115888337719, 1125899906842624) |
Raymond Hettinger | 9fce0ba | 2008-10-05 16:46:29 +0000 | [diff] [blame] | 169 | |
| 170 | Since the ratio is exact, it can be used to losslessly recreate the |
| 171 | original value:: |
| 172 | |
| 173 | >>> x == 3537115888337719 / 1125899906842624 |
| 174 | True |
| 175 | |
| 176 | The :meth:`float.hex` method expresses a float in hexadecimal (base |
| 177 | 16), again giving the exact value stored by your computer:: |
| 178 | |
| 179 | >>> x.hex() |
| 180 | '0x1.921f9f01b866ep+1' |
| 181 | |
| 182 | This precise hexadecimal representation can be used to reconstruct |
| 183 | the float value exactly:: |
| 184 | |
| 185 | >>> x == float.fromhex('0x1.921f9f01b866ep+1') |
| 186 | True |
| 187 | |
| 188 | Since the representation is exact, it is useful for reliably porting values |
| 189 | across different versions of Python (platform independence) and exchanging |
| 190 | data with other languages that support the same format (such as Java and C99). |
| 191 | |
Raymond Hettinger | 5afb5c6 | 2009-04-26 22:01:46 +0000 | [diff] [blame] | 192 | Another helpful tool is the :func:`math.fsum` function which helps mitigate |
| 193 | loss-of-precision during summation. It tracks "lost digits" as values are |
| 194 | added onto a running total. That can make a difference in overall accuracy |
| 195 | so that the errors do not accumulate to the point where they affect the |
| 196 | final total: |
| 197 | |
| 198 | >>> sum([0.1] * 10) == 1.0 |
| 199 | False |
| 200 | >>> math.fsum([0.1] * 10) == 1.0 |
| 201 | True |
Raymond Hettinger | 9fce0ba | 2008-10-05 16:46:29 +0000 | [diff] [blame] | 202 | |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 203 | .. _tut-fp-error: |
| 204 | |
| 205 | Representation Error |
| 206 | ==================== |
| 207 | |
| 208 | This section explains the "0.1" example in detail, and shows how you can perform |
| 209 | an exact analysis of cases like this yourself. Basic familiarity with binary |
| 210 | floating-point representation is assumed. |
| 211 | |
| 212 | :dfn:`Representation error` refers to the fact that some (most, actually) |
| 213 | decimal fractions cannot be represented exactly as binary (base 2) fractions. |
| 214 | This is the chief reason why Python (or Perl, C, C++, Java, Fortran, and many |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 215 | others) often won't display the exact decimal number you expect. |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 216 | |
| 217 | Why is that? 1/10 is not exactly representable as a binary fraction. Almost all |
| 218 | machines today (November 2000) use IEEE-754 floating point arithmetic, and |
| 219 | almost all platforms map Python floats to IEEE-754 "double precision". 754 |
| 220 | doubles contain 53 bits of precision, so on input the computer strives to |
Benjamin Peterson | 5c6d787 | 2009-02-06 02:40:07 +0000 | [diff] [blame] | 221 | convert 0.1 to the closest fraction it can of the form *J*/2**\ *N* where *J* is |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 222 | an integer containing exactly 53 bits. Rewriting :: |
| 223 | |
| 224 | 1 / 10 ~= J / (2**N) |
| 225 | |
| 226 | as :: |
| 227 | |
| 228 | J ~= 2**N / 10 |
| 229 | |
| 230 | and recalling that *J* has exactly 53 bits (is ``>= 2**52`` but ``< 2**53``), |
| 231 | the best value for *N* is 56:: |
| 232 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 233 | >>> 2**52 <= 2**56 // 10 < 2**53 |
| 234 | True |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 235 | |
| 236 | That is, 56 is the only value for *N* that leaves *J* with exactly 53 bits. The |
| 237 | best possible value for *J* is then that quotient rounded:: |
| 238 | |
| 239 | >>> q, r = divmod(2**56, 10) |
| 240 | >>> r |
Georg Brandl | bae1b94 | 2008-08-10 12:16:45 +0000 | [diff] [blame] | 241 | 6 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 242 | |
| 243 | Since the remainder is more than half of 10, the best approximation is obtained |
| 244 | by rounding up:: |
| 245 | |
| 246 | >>> q+1 |
Georg Brandl | bae1b94 | 2008-08-10 12:16:45 +0000 | [diff] [blame] | 247 | 7205759403792794 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 248 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 249 | Therefore the best possible approximation to 1/10 in 754 double precision is:: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 250 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 251 | 7205759403792794 / 2 ** 56 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 252 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 253 | Dividing both the numerator and denominator by two reduces the fraction to:: |
| 254 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 255 | 3602879701896397 / 2 ** 55 |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 256 | |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 257 | Note that since we rounded up, this is actually a little bit larger than 1/10; |
| 258 | if we had not rounded up, the quotient would have been a little bit smaller than |
| 259 | 1/10. But in no case can it be *exactly* 1/10! |
| 260 | |
| 261 | So the computer never "sees" 1/10: what it sees is the exact fraction given |
| 262 | above, the best 754 double approximation it can get:: |
| 263 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 264 | >>> 0.1 * 2 ** 55 |
| 265 | 3602879701896397.0 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 266 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 267 | If we multiply that fraction by 10\*\*55, we can see the value out to |
| 268 | 55 decimal digits:: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 269 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 270 | >>> 3602879701896397 * 10 ** 55 // 2 ** 55 |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 271 | 1000000000000000055511151231257827021181583404541015625 |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 272 | |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 273 | meaning that the exact number stored in the computer is equal to |
| 274 | the decimal value 0.1000000000000000055511151231257827021181583404541015625. |
| 275 | Instead of displaying the full decimal value, many languages (including |
| 276 | older versions of Python), round the result to 17 significant digits:: |
| 277 | |
| 278 | >>> format(0.1, '.17f') |
| 279 | '0.10000000000000001' |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 280 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 281 | The :mod:`fractions` and :mod:`decimal` modules make these calculations |
| 282 | easy:: |
Georg Brandl | 116aa62 | 2007-08-15 14:28:22 +0000 | [diff] [blame] | 283 | |
Raymond Hettinger | 8bd1d4f | 2009-04-24 03:09:06 +0000 | [diff] [blame] | 284 | >>> from decimal import Decimal |
| 285 | >>> from fractions import Fraction |
Raymond Hettinger | b8aff6a | 2009-06-29 01:50:51 +0000 | [diff] [blame] | 286 | |
| 287 | >>> Fraction.from_float(0.1) |
| 288 | Fraction(3602879701896397, 36028797018963968) |
| 289 | |
| 290 | >>> (0.1).as_integer_ratio() |
| 291 | (3602879701896397, 36028797018963968) |
| 292 | |
| 293 | >>> Decimal.from_float(0.1) |
| 294 | Decimal('0.1000000000000000055511151231257827021181583404541015625') |
| 295 | |
| 296 | >>> format(Decimal.from_float(0.1), '.17') |
| 297 | '0.10000000000000001' |