{"id":508,"date":"2015-06-14T00:12:51","date_gmt":"2015-06-14T00:12:51","guid":{"rendered":"http:\/\/www.vidarholen.net\/contents\/blog\/?p=508"},"modified":"2015-06-14T00:12:51","modified_gmt":"2015-06-14T00:12:51","slug":"technically-correct-floating-point-calculations-in-bc","status":"publish","type":"post","link":"https:\/\/www.vidarholen.net\/contents\/blog\/?p=508","title":{"rendered":"Technically correct: floating point calculations in bc"},"content":{"rendered":"

Whenever someone asks how to do floating point math in a shell script, the answer is typically bc<\/code>:<\/p>\n

\r\n$  echo \"scale=9; 22\/7\" | bc\r\n3.142857142\r\n<\/pre>\n
However, this is technically wrong: bc<\/code> does not support floating point at all! What you see above is arbitrary precision FIXED point arithmetic. <\/p>\n
The user’s intention is obviously to do math with fractional numbers, regardless of the low level implementation, so the above is a good and pragmatic answer. However, technically correct is the best kind of correct, so let’s stop being helpful and start pedantically splitting hairs instead!<\/p>\n
Fixed vs floating point<\/b><\/p>\n
There are many important things that every programmer should know about floating point<\/a>, but in one sentence, the larger they get, the less precise they are.<\/p>\n
In fixed point you have a certain number of digits, and a decimal point fixed in place like on a tax form: 001234.56<\/code>. No matter how small or large the number, you can always write down increments of 0.01, whether it’s 000000.01 or 999999.99.<\/p>\n
Floating point, meanwhile, is basically scientific notation. If you have 1.23e-4 (0.000123), you can increment by a millionth to get 1.24e-4. However, if you have 1.23e4 (12300), you can’t add less than 100 unless you reserve more space for more digits.<\/p>\n
We can see this effect in practice in any language that supports floating point, such as Haskell:<\/p>\n
\r\n> truncate (16777216 - 1 :: Float)\r\n16777215\r\n> truncate (16777216 + 1 :: Float)\r\n16777216\r\n<\/pre>\n
Subtracting 1 gives us the decremented number, but adding 1 had no effect with floating point math! bc, with its arbitrary precision fixed points, would instead correctly give us 16777217! This is clearly unacceptable!<\/p>\n
Floating point in bc<\/b><\/p>\n
The problem with the bc<\/code> solution is, in other words, that the math is too correct. Floating point math always introduces and accumulates rounding errors in ways that are hard to predict. Fixed point doesn’t, and therefore we need to find a way to artificially introduce the same type of inaccuracies! We can do this by rounding a number to a N significant bits, where N = 24 for float<\/code> and 52 for double<\/code>. Here is some bc code for that:<\/p>\n
\r\nscale=30\r\n\r\ndefine trunc(x) {\r\n  auto old, tmp\r\n  old=scale; scale=0; tmp=x\/1; scale=old\r\n  return tmp\r\n}\r\ndefine fp(bits, x) {\r\n  auto i\r\n  if (x < 0) return -fp(bits, -x);\r\n  if (x == 0) return 0;\r\n  i=bits\r\n  while (x < 1) { x*=2; i+=1; }\r\n  while (x >= 2) { x\/=2; i-=1; }\r\n  return trunc(x * 2^bits + 0.5) \/ 2^(i)\r\n}\r\n\r\ndefine float(x) { return fp(24, x); }\r\ndefine double(x) { return fp(52, x); }\r\ndefine test(x) {\r\n  print \"Float:  \", float(x), \"\\n\"\r\n  print \"Double: \", double(x), \"\\n\"\r\n}\r\n<\/pre>\n
With this file named fp<\/code>, we can try it out:<\/p>\n
\r\n$ bc -ql fp <<< \"22\/7\"\r\n3.142857142857142857142857142857\r\n\r\n$ bc -ql fp <<< \"float(22\/7)\"\r\n3.142857193946838378906250000000\r\n<\/pre>\n
The first number is correct to 30 decimals. Yuck! However, with our floating point simulator applied, we get the desired floating point style errors after ~7 decimals!<\/p>\n
Let's write a similar program for doing the same thing but with actual floating point, printing them out up to 30 decimals as well:<\/p>\n
\r\n{-# LANGUAGE RankNTypes #-}\r\nimport Control.Monad\r\nimport Data.Number.CReal\r\nimport System.Environment\r\n\r\nmain = do\r\n    input <- liftM head getArgs\r\n    putStrLn . (\"Float:  \" ++) $ showNumber (read input :: Float)\r\n    putStrLn . (\"Double: \" ++) $ showNumber (read input :: Double)\r\n  where\r\n    showNumber :: forall a. Real a => a -> String\r\n    showNumber = showCReal 30 . realToFrac\r\n<\/pre>\n
Here's a comparison of the two:<\/p>\n
\r\n$ bc -ql fp <<< \"x=test(1000000001.3)\"\r\nFloat:  1000000000.000000000000000000000000000000\r\nDouble: 1000000001.299999952316284179687500000000\r\n\r\n$ .\/fptest 1000000001.3\r\nFloat:  1000000000.0\r\nDouble: 1000000001.2999999523162841796875\r\n<\/pre>\n
Due to differences in rounding and\/or off-by-one bugs, they're not always identical like here, but the error bars are similar. <\/p>\n
Now we can finally start doing floating point math in bc!<\/p>\n","protected":false},"excerpt":{"rendered":"
Whenever someone asks how to do floating point math in a shell script, the answer is typically bc: $ echo “scale=9; 22\/7” | bc 3.142857142 However, this is technically wrong: bc does not support floating point at all! What you see above is arbitrary precision FIXED point arithmetic. The user’s intention is obviously to do … Continue reading “Technically correct: floating point calculations in bc”<\/span><\/a><\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":"","jetpack_publicize_message":"","jetpack_is_tweetstorm":false,"jetpack_publicize_feature_enabled":true},"categories":[5,4,23],"tags":[11,51,40],"class_list":["post-508","post","type-post","status-publish","format-standard","hentry","category-advanced-linux","category-linux","category-programming","tag-bash","tag-bc","tag-technically-correct"],"jetpack_featured_media_url":"","_links":{"self":[{"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/posts\/508","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=508"}],"version-history":[{"count":29,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/posts\/508\/revisions"}],"predecessor-version":[{"id":537,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=\/wp\/v2\/posts\/508\/revisions\/537"}],"wp:attachment":[{"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=508"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=508"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.vidarholen.net\/contents\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=508"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}