{"id":1245,"date":"2015-10-13T22:23:36","date_gmt":"2015-10-14T05:23:36","guid":{"rendered":"http:\/\/www.elbeno.com\/blog\/?p=1245"},"modified":"2015-10-15T08:54:16","modified_gmt":"2015-10-15T15:54:16","slug":"more-constexpr-floating-point-computation","status":"publish","type":"post","link":"https:\/\/www.elbeno.com\/blog\/?p=1245","title":{"rendered":"More constexpr floating-point computation"},"content":{"rendered":"

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

In the last post<\/a>, I covered my first forays into writing C++11-style constexpr<\/code> floating point math functions. Once I’d done some trig functions, exp<\/code>, and floor<\/code> and friends, it seemed like the rest would be a piece of cake. Not entirely…<\/p>\n

But there was some easy stuff. With trunc<\/code> sorted out, fmod<\/code> was just:<\/p>\n

\r\nconstexpr float fmod(float x, float y)\r\n{\r\n  return y != 0 ? x - trunc(x\/y)*y :\r\n    throw err::fmod_domain_error;\r\n}\r\n<\/pre>\n
And remainder<\/code> was very similar. Also, fmax<\/code>, fmin<\/code> and fdim<\/code> were trivial: I won’t bore you with the code.<\/p>\n
Euler was quite good at maths<\/strong><\/p>\n
Now I come from a games background, and there are a few maths functions that get used all the time in games. I already have sqrt<\/code>. What I want now is atan<\/code> and atan2<\/code>. To warm up, I implemented the infinite series calculations for asin<\/code> and acos<\/code>, each of which are only defined for an input |x| <= 1. That was fairly straightforward, broadly in line with what I'd done before for the series expansions of sin<\/code> and cos<\/code>.<\/p>\n
Inverting tan<\/code> was another matter: but a search revealed a way to do it<\/a> using asin<\/code>:<\/p>\n
[latex]arctan(x) = arcsin(\\frac{x}{\\sqrt{x^2+1}})[\/latex]<\/p>\n
Which I went with for a while. But then I spotted<\/a> the magic words:<\/p>\n
“Leonhard Euler found a more efficient series for the arctangent”<\/p>\n
[latex]arctan(z) = \\frac{z}{1+z^2} \\sum_{n=0}^{\\infty} \\prod_{k=1}^n \\frac{2kz^2}{(2k+1)(1+z^2)}[\/latex]<\/p>\n
More efficient (which I take here to mean more quickly converging)? And discovered by Euler? That’s good enough for me.<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T atan_term(T x2, int k)\r\n  {\r\n    return (T{2}*static_cast(k)*x2)\r\n      \/ ((T{2}*static_cast(k)+T{1}) * (T{1}+x2));\r\n  }\r\n  template \r\n  constexpr T atan_product(T x, int k)\r\n  {\r\n    return k == 1 ? atan_term(x*x, k) :\r\n      atan_term(x*x, k) * atan_product(x, k-1);\r\n  }\r\n  template \r\n  constexpr T atan_sum(T x, T sum, int n)\r\n  {\r\n    return sum + atan_product(x, n) == sum ?\r\n      sum :\r\n      atan_sum(x, sum + atan_product(x, n), n+1);\r\n  }\r\n}\r\ntemplate \r\nconstexpr T atan(T x)\r\n{\r\n  return true ?\r\n    x \/ (T{1} + x*x) * detail::atan_sum(x, T{1}, 1) :\r\n    throw err::atan_runtime_error;\r\n}\r\n<\/pre>\n
And atan2<\/code> follows from atan<\/code> relatively easily.<\/p>\n
It’s big, it’s heavy, it’s wood<\/strong><\/p>\n
Looking over what I’d implemented so far, the next important function (and one that would serve as a building block for a few others) was log<\/code>. So I programmed it using my friend the Taylor series expansion.<\/p>\n
Hm. It turns out that a naive implementation of log<\/code> is quite slow to converge. (Yes, what a surprise.) OK, a bit of searching yields fruit<\/a>: “… use Newton’s method to invert the exponentional function … the iteration simplifies to:”<\/p>\n
[latex]y_{n+1} = y_n + 2 \\frac{x – e^{y_n}}{x + e^{y_n}}[\/latex]<\/p>\n
“… which has cubic convergence to ln(x).”<\/p>\n
Result! I coded it up.<\/p>\n
\r\nnamespace detail\r\n{\r\n  template \r\n  constexpr T log_iter(T x, T y)\r\n  {\r\n    return y + T{2} * (x - cx::exp(y)) \/ (x + cx::exp(y));\r\n  }\r\n  template \r\n  constexpr T log(T x, T y)\r\n  {\r\n    return feq(y, log_iter(x, y)) ? y : log(x, log_iter(x, y));\r\n  }\r\n}\r\n<\/pre>\n
It worked well. A week or so later, an article about C++11 constexpr log<\/a> surfaced serendipitously, and on reading it I realised my implementation wasn’t very well tested, so I added some tests. The tests revealed numerical instability for larger input values, so I took a tip from the article and used a couple of identities to constrain the input:<\/p>\n
[latex]ln(x) = ln(\\frac{x}{e}) + 1[\/latex]
\nand
\n[latex]ln(x) = ln(e.x) – 1[\/latex]<\/p>\n
Coding these expressions to constrain the input proper to detail::log<\/code> seemed to solve the instability issues.<\/p>\n
Are we done yet?<\/strong><\/p>\n
Now I was nearing the end of what I felt was a quorum of functions. To wit:<\/p>\n