{"id":1211,"date":"2015-10-13T21:17:41","date_gmt":"2015-10-14T04:17:41","guid":{"rendered":"http:\/\/www.elbeno.com\/blog\/?p=1211"},"modified":"2015-10-15T08:55:10","modified_gmt":"2015-10-15T15:55:10","slug":"floating-point-maths-constexpr-style","status":"publish","type":"post","link":"https:\/\/www.elbeno.com\/blog\/?p=1211","title":{"rendered":"Floating-point maths, constexpr style"},"content":{"rendered":"

(Start at the beginning of the series<\/a> – and all the source can be found in my github repo<\/a>)<\/p>\n

To ease into constexpr programming I decided to tackle some floating-point maths functions. Disclaimer: I’m not a mathematician and this code has not been rigorously tested for numeric stability or convergence in general. I wrote it more for my own learning than for serious mathematical purposes, and anyway, it’s perhaps likely that constexpr<\/code> maths functions will be in the standard library sooner or later (and possibly implemented by intrinsics). However, I did subject it to many ad-hoc tests using static_assert<\/code>. One of the nice things about compile-time computation: the tests are compile-time too, and if it compiles, one can be reasonably sure it works (at least for the values tested)!<\/p>\n

Anyway, I went to cppreference.com’s list of common mathematical functions<\/a> to see what to tackle. Starting with abs<\/code> is trivial:<\/p>\n

Starting with abs<\/code> is trivial:<\/p>\n

\r\ntemplate \r\nconstexpr T abs(T x)\r\n{\r\n  return x >= 0 ? x :\r\n    x < 0 ? -x :\r\n    throw err::abs_runtime_error;\r\n}\r\n<\/pre>\n
For the sake of clarity, I'm omitting some enable_if<\/code> boilerplate here that ensures that the argument has the right type.<\/p>\n
Iterative methods FTW<\/strong><\/p>\n
Now, other functions take a little thought. But remembering calculus, many of them are susceptible to iterative methods, so my general plan formed: formulate a recursive calculation which converges to the real answer, and terminate the recursion when the value is within epsilon of the answer.<\/p>\n
For example, to compute sqrt<\/code>, we can use the well-known Newton-Raphson method, where we start with a guess g0<\/sub> as an approximation to √x, then:<\/p>\n
[latex]g_{n+1} = \\frac{(g_n + \\frac{x}{g_n})}{2}[\/latex]<\/p>\n
First, we need a function to determine termination, which I put in a detail<\/code> namespace:<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr bool feq(T x, T y)\r\n  {\r\n    return abs(x - y) <=\r\n      std::numeric_limits::epsilon();\r\n  }\r\n}\r\n<\/pre>\n
And then sqrt<\/code> is straightforward to express (using the value itself as the initial guess), with a driver function handling the domain of the argument and the throw<\/code> pattern, and delegating the recursion to a function in the detail<\/code> namespace.<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T sqrt(T x, T guess)\r\n  {\r\n    return feq(guess, (guess + x\/guess)\/T{2}) ? guess :\r\n      sqrt(x, (guess + x\/guess)\/T{2});\r\n  }\r\n}\r\ntemplate \r\nconstexpr T sqrt(T x)\r\n{\r\n  return x == 0 ? 0 :\r\n    x > 0 ? detail::sqrt(x, x) :\r\n    throw err::sqrt_domain_error;\r\n}\r\n<\/pre>\n
This ends up being a useful pattern for many functions. For example, cbrt<\/code> is practically identical; only the equational details change. And once we have sqrt<\/code>, hypot<\/code> is trivial of course.<\/p>\n
Next I tackled trigonometric functions. Remembering Taylor series and turning to that Oracle of All Things Maths, Wolfram Alpha<\/a>, we can find the Taylor series expansions:<\/p>\n
[latex]sin(x) = \\sum_{k=0}^{\\infty} \\frac{(-1)^k x^{1+2k}}{(1+2k)!}[\/latex]
\nand
\n[latex]cos(x) = \\sum_{k=0}^{\\infty} \\frac{(-1)^k x^{2k}}{(2k)!}[\/latex]<\/p>\n
These are basically the same formula, and with some massaging, they can share a common recursive function with slightly different initial constants:<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T trig_series(T x, T sum, T n, int i, int s, T t)\r\n  {\r\n    return feq(sum, sum + t*s\/n) ?\r\n      sum :\r\n      trig_series(x, sum + t*s\/n, n*i*(i+1), i+2, -s, t*x*x);\r\n  }\r\n}\r\ntemplate \r\nconstexpr T sin(T x)\r\n{\r\n  return true ?\r\n    detail::trig_series(x, x, T{6}, 4, -1, x*x*x) :\r\n    throw err::sin_runtime_error;\r\n}\r\ntemplate \r\nconstexpr T cos(T x)\r\n{\r\n  return true ?\r\n    detail::trig_series(x, T{1}, T{2}, 3, -1, x*x) :\r\n    throw err::cos_runtime_error;\r\n}\r\n<\/pre>\n
Once we have sin<\/code> and cos<\/code>, tan<\/code> is trivial of course. And exp<\/code> has almost the same series as sin<\/code> and cos<\/code> - easier, actually:<\/p>\n
[latex]e^x = \\sum_{k=0}^{\\infty} \\frac{x^k}{k!}[\/latex]<\/p>\n
Which yields in code:<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T exp(T x, T sum, T n, int i, T t)\r\n  {\r\n    return feq(sum, sum + t\/n) ?\r\n      sum :\r\n      exp(x, sum + t\/n, n * i, i+1, t * x);\r\n  }\r\n}\r\ntemplate \r\nconstexpr T exp(T x)\r\n{\r\n  return true ? detail::exp(x, T{1}, T{1}, 2, x) :\r\n    throw err::exp_runtime_error;\r\n}\r\n<\/pre>\n
No peeking at the floating-point details<\/strong><\/p>\n
Now, how about floor<\/code>, ceil<\/code>, trunc<\/code> and round<\/code>? They are all variations on a theme: solving one basically means solving all, with only slight differences for things like negative numbers. I thought about this for a half hour - there wasn't an obvious solution to me at first. We can't trivially cast to an integral type, because there is no integral type big enough to hold a floating point value in general. And in the name of strict portability, we can't be sure that the floating point format is IEEE 754 standard (although I think I am prepared to accept this as a limitation). So even if we could fiddle with the representation, portability goes out the window. Finally, constexpr<\/code> (including C++14 constexpr<\/code>) forbids reinterpret_cast<\/code> and other tricks like access through union<\/code> members.<\/p>\n
So how to compute floor<\/code>? After a little thought I hit upon an approach. In fact, I am relying on the IEEE 754 standard in a small way - the fact that any integer power of 2 is exactly representable as floating point. With that I formed the plan:<\/p>\n
    \n
  1. Start with a guess of zero.<\/li>\n
  2. Start with an increment: 2(std::numeric_limits::max_exponent<\/code> - 1)<\/sup>.<\/li>\n
  3. If (guess + increment) is larger than x, halve the increment.<\/li>\n
  4. Otherwise, add the increment to the guess, and that's the new guess.<\/li>\n
  5. Stop when the increment is less than 2.<\/li>\n<\/ol>\n
    This gives a binary-search like approach to finding the correct value for floor<\/code>. And the code is easy to write. First we need a function to raise a floating point value to an integer power, in the standard O(log n) way:<\/p>\n
    \r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T ipow(T x, int n)\r\n  {\r\n    return (n == 0) ? T{1} :\r\n      n == 1 ? x :\r\n      n > 1 ? ((n & 1) ? x * ipow(x, n-1) :\r\n                         ipow(x, n\/2) * ipow(x, n\/2)) :\r\n      T{1} \/ ipow(x, -n);\r\n  }\r\n}\r\n<\/pre>\n
    Then the actual floor<\/code> function, which uses ceil<\/code> for the case of negative numbers:<\/p>\n
    \r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T floor(T x, T guess, T inc)\r\n  {\r\n    return guess + inc <= x ? floor(x, guess + inc, inc) :\r\n      inc <= T{1} ? guess : floor(x, guess, inc\/T{2});\r\n  }\r\n}\r\nconstexpr float floor(float x)\r\n{\r\n  return x < 0 ? -ceil(-x) :\r\n    x >= 0 ? detail::floor(\r\n        x, 0.0f,\r\n        detail::ipow(2.0f,\r\n            std::numeric_limits::max_exponent-1)) :\r\n    throw err::floor_runtime_error;\r\n}\r\n<\/pre>\n
    The other functions ceil<\/code>, trunc<\/code> and round<\/code> are very similar.<\/p>\n
    That's some deep recursion<\/strong><\/p>\n
    Now, there is a snag with this plan. It works well enough for float<\/code>, but when we try it with double<\/code>, it falls down. Appendix B of the C++ standard (Implementation quantities) recommends that compilers offer a minimum recursive depth of 512 for constexpr<\/code> function invocations. And on my machine, std::numeric_limits<float>::max_exponent<\/code> is 128. But for double<\/code>, it's 1024. And for long double<\/code>, it's 16384. Since by halving the increment each time, we're basically doing a linear search on the exponent, 512 recursions isn't enough. I'm prepared to go with what the standard recommends in the sense that I don't mind massaging switches for compilers that don't follow the recommendations, but I'd rather not mess around with compiler switches for those that do. So how can I get around this?<\/p>\n
    Well, I'm using a binary search to cut down the increment. What if I increase the fanout so it's not a binary search, but a trinary search or more? How big of a fanout do I need? What I need is:<\/p>\n
    max_recursion_depth > max_exponent\/x, where 2x<\/sup> is the fanout<\/p>\n
    So to deal with double<\/code>, x = 3 and we need a fanout of 8. An octary search? Well, it's easy enough to code, even if it is a nested-ternary-operator-from-hell:<\/p>\n
    \r\ntemplate \r\nconstexpr T floor8(T x, T guess, T inc)\r\n{\r\n  return\r\n    inc < T{8} ? floor(x, guess, inc) :\r\n    guess + inc <= x ? floor8(x, guess + inc, inc) :\r\n    guess + (inc\/T{8})*T{7} <= x ? floor8(x, guess + (inc\/T{8})*T{7}, inc\/T{8}) :\r\n    guess + (inc\/T{8})*T{6} <= x ? floor8(x, guess + (inc\/T{8})*T{6}, inc\/T{8}) :\r\n    guess + (inc\/T{8})*T{5} <= x ? floor8(x, guess + (inc\/T{8})*T{5}, inc\/T{8}) :\r\n    guess + (inc\/T{8})*T{4} <= x ? floor8(x, guess + (inc\/T{8})*T{4}, inc\/T{8}) :\r\n    guess + (inc\/T{8})*T{3} <= x ? floor8(x, guess + (inc\/T{8})*T{3}, inc\/T{8}) :\r\n    guess + (inc\/T{8})*T{2} <= x ? floor8(x, guess + (inc\/T{8})*T{2}, inc\/T{8}) :\r\n    guess + inc\/T{8} <= x ? floor8(x, guess + inc\/T{8}, inc\/T{8}) :\r\n    floor8(x, guess, inc\/T{8});\r\n}\r\n<\/pre>\n
    (Apologies for the width, but I really can't make it much more readable.) For the base case, I revert to the regular binary search implementation of floor<\/code>. What about for long double<\/code>? Well, since max_exponent is 16384, x = 33 and we need a fanout of 233<\/sup>. That's not happening! For long double<\/code>, I have no choice but to revert to the C++14 constexpr<\/code> rules:<\/p>\n
    \r\nconstexpr long double floor(long double x)\r\n{\r\n  if (x < 0.0) return -ceil(-x);\r\n  long double inc = detail::ipow(\r\n      2.0l,\r\n      std::numeric_limits::max_exponent - 1);\r\n  long double guess = 0.0l;\r\n  for (;;)\r\n  {\r\n    while (guess + inc > x)\r\n    {\r\n      inc \/= 2.0l;\r\n      if (inc < 1.0l)\r\n        return guess;\r\n    }\r\n    guess += inc;\r\n  }\r\n  throw err::floor_runtime_error;\r\n}\r\n<\/pre>\n
    And that's enough of such stuff for one blog post. Next time<\/a>, I'll go through the rest of the maths functions I implemented, with help from Leonhard Euler and Wikipedia.<\/p>\n","protected":false},"excerpt":{"rendered":"
    (Start at the beginning of the series – and all the source can be found in my github repo) To ease into constexpr programming I decided to tackle some floating-point maths functions. Disclaimer: I’m not a mathematician and this code has not been rigorously tested for numeric stability or convergence…<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[22],"tags":[],"class_list":["post-1211","post","type-post","status-publish","format-standard","hentry","category-cpp"],"_links":{"self":[{"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/posts\/1211","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=1211"}],"version-history":[{"count":39,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/posts\/1211\/revisions"}],"predecessor-version":[{"id":1303,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=\/wp\/v2\/posts\/1211\/revisions\/1303"}],"wp:attachment":[{"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=1211"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=1211"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.elbeno.com\/blog\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=1211"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}