{"id":3839,"date":"2018-02-28T15:50:46","date_gmt":"2018-02-28T15:50:46","guid":{"rendered":"http:\/\/www.blopig.com\/blog\/?p=3839"},"modified":"2018-03-07T15:53:23","modified_gmt":"2018-03-07T15:53:23","slug":"a-short-intro-to-machine-precision-and-how-to-beat-it","status":"publish","type":"post","link":"https:\/\/www.blopig.com\/blog\/2018\/02\/a-short-intro-to-machine-precision-and-how-to-beat-it\/","title":{"rendered":"A short intro to machine precision and how to beat it"},"content":{"rendered":"

Most people who\u2019ve ever sat through a Rounding Error is Bad lecture will be familiar with the following example:<\/p>\n

> (0.1+0.1+0.1) == 0.3\r\nFALSE<\/pre>\n
The reason this is so unsettling is because most of the time we think about numbers in base-10. This means we use ten digits \"\{0,, and we perform arithmetic based on this ten digit notation. This doesn\u2019t always matter much for pen and paper maths but it\u2019s an integral part of how we think about more complex operations and in particular how we think about accuracy. We see \"0.1\" as a finite decimal fraction, and so it\u2019s only natural that we should be able to do accurate sums with it. And if we can do simple arithmetic, then surely computers can too? In this blog post I\u2019m going to try and briefly explain what causes rounding errors such as the one above, and how we might get away with going beyond machine precision.<\/p>\n
Take a number \"x, say \"x=1\/3\". The decimal representation of \"x\" is of the form \"x=\sum_{i=1}^{\infty}. The \"a_i here are the digits that go after the radix point. In the case of \"x=1\/3\" these are all equal \"a_i=3\", or \"x=0.333\dots. Some numbers, such as our favourite \"x\", don\u2019t have a finite decimal expansion. Others, such as \"0.3\", do, meaning that after some \"i, all \"a_{i+j}=0\". When we talk about rounding errors and accuracy, what we actually mean is that we only care about the first few digits, say \"i\leq, and we\u2019re happy to approximate to \"x\approx, potentially rounding up at the last digit.<\/p>\n
Computers, on the other hand, store numbers in base-2 rather than base-10, which means that they use a different series expansion \"x=\sum_{i=1}^{\infty}, \"b_i to represent the same number. Our favourite number \"x\" is actually stored as \"0.1010101\dots rather than \"0.3333333\dots, despite the fact it appears as the latter on a computer screen. Crucially, arithmetic is done in base-2 and, since only a finite number of binary digits are stored (\"i\leq for most purposes these days), rounding errors also occur in base-2.<\/p>\n
All numbers with a finite binary expansion, such as \"0.25_{10}=0\times also have a finite decimal expansion, meaning we can do accurate arithmetic with them in both systems. However, the reverse isn\u2019t true, which is what causes the issue with \"0.1+0.1+0.1\neq. In binary, the nice and tidy \"0.1_{10}. We observe the rounding error because unlike us, the computer is trying to sum over infinite series.<\/p>\n
While it\u2019s not possible to do infinite sums with finite resources, there is a way to go beyond machine precision if you wanted to, at least for rational \"x=p\/q\", where \"p,. In the example above, the issue comes from dividing by 10 on each side of the (in)equality. Luckily for us, we can avoid doing so. Integer arithmetic is easy in any base, and so<\/p>\n
> (1+1+1) == 3 \r\nTRUE<\/pre>\n
Shocking, I know. On a more serious note, it is possible to write an algorithm which calculates the binary expansion of \"x=p\/q\" using only integer arithmetic. Usually binary expansion happens in this way:<\/p>\n
set x, maxIter\r\ninitialise b, i=1\r\nwhile x>0 AND i<=maxIter {\r\n if 2*x>=1\r\n    b[i]=1\r\n else\r\n    b[i]=0\r\n x = 2*x-b[i]\r\n i = i+1 \r\n}\r\nreturn b<\/pre>\n
Problems arise whenever we try to compute something non-integer (the highlighted lines 3, 4, and 8). However, we can rewrite these using \"x and\u00a0 shifting the division by \"q\" to the right-hand side of each inequality \/ assignment operator:<\/p>\n
set p, q, maxIter\r\ninitialise b, i=1 \r\nwhile p>0 AND i<=maxIter { \r\n if 2*p>=q \r\n    b[i]=1 \r\n else \r\n    b[i]=0 \r\n p = 2*p-b[i]*q \r\n i = i+1 \r\n} \r\nreturn b<\/pre>\n
Provided we’re not dealing with monstrously large integers (i.e. as long as we can safely double \"p\"), implementing the above lets us compute \"p\/q\" with arbitrary precision given by maxIter<\/code>. So we can beat machine precision for rationals! And the combination of arbitrarily accurate rationals and arbitrarily accurate series approximations (think Riemann zeta function, for example) means we can also get the occasional arbitrarily accurate irrational.<\/p>\n
To sum up, rounding errors are annoying, partly because it’s not always intuitive when and how they happen. As a general rule the best way to avoid them is to make your computer do as little work as possible, and to avoid non-integer calculations whenever you can. But you already knew that, didn’t you?<\/p>\n
This post was partially inspired by the undergraduate course on Simulation and Statistical Programming lectured by Prof Julien Berestycki and Prof Robin Evans. It was also inspired by my former maths teacher who used to mark us down for doing more work than necessary even when our solutions were correct. He had a point.
\n<\/em><\/p>\n","protected":false},"excerpt":{"rendered":"
Most people who\u2019ve ever sat through a Rounding Error is Bad lecture will be familiar with the following example: > (0.1+0.1+0.1) == 0.3 FALSE The reason this is so unsettling is because most of the time we think about numbers in base-10. This means we use ten digits , and we perform arithmetic based on […]<\/p>\n","protected":false},"author":40,"featured_media":0,"comment_status":"closed","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","wikipediapreview_detectlinks":true,"_monsterinsights_skip_tracking":false,"_monsterinsights_sitenote_active":false,"_monsterinsights_sitenote_note":"","_monsterinsights_sitenote_category":0,"ngg_post_thumbnail":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[15],"tags":[],"ppma_author":[529],"class_list":["post-3839","post","type-post","status-publish","format-standard","hentry","category-technical"],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"authors":[{"term_id":529,"user_id":40,"is_guest":0,"slug":"lyuba","display_name":"Lyuba","avatar_url":"https:\/\/secure.gravatar.com\/avatar\/e10fff614041616d2611ef048154e92d6d6a208f36f8634b41a9abcec271c8c1?s=96&d=mm&r=g","0":null,"1":"","2":"","3":"","4":"","5":"","6":"","7":"","8":""}],"_links":{"self":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/3839","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/users\/40"}],"replies":[{"embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/comments?post=3839"}],"version-history":[{"count":15,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/3839\/revisions"}],"predecessor-version":[{"id":3862,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/posts\/3839\/revisions\/3862"}],"wp:attachment":[{"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/media?parent=3839"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/categories?post=3839"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/tags?post=3839"},{"taxonomy":"author","embeddable":true,"href":"https:\/\/www.blopig.com\/blog\/wp-json\/wp\/v2\/ppma_author?post=3839"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}