## More constexpr floating-point computation

October 13th, 2015(Start at the beginning of the series – and all the source can be found in my github repo)

In the last post, I covered my first forays into writing C++11-style `constexpr`

floating point math functions. Once I’d done some trig functions, `exp`

, and `floor`

and friends, it seemed like the rest would be a piece of cake. Not entirely…

But there was some easy stuff. With `trunc`

sorted out, `fmod`

was just:

constexpr float fmod(float x, float y) { return y != 0 ? x - trunc(x/y)*y : throw err::fmod_domain_error; } |

And `remainder`

was very similar. Also, `fmax`

, `fmin`

and `fdim`

were trivial: I won’t bore you with the code.

**Euler was quite good at maths**

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`

. What I want now is `atan`

and `atan2`

. To warm up, I implemented the infinite series calculations for `asin`

and `acos`

, 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`

and `cos`

.

Inverting `tan`

was another matter: but a search revealed a way to do it using `asin`

:

Which I went with for a while. But then I spotted the magic words:

“Leonhard Euler found a more efficient series for the arctangent”

More efficient (which I take here to mean more quickly converging)? And discovered by Euler? That’s good enough for me.

namespace detail { template <typename T> constexpr T atan_term(T x2, int k) { return (T{2}*static_cast<T>(k)*x2) / ((T{2}*static_cast<T>(k)+T{1}) * (T{1}+x2)); } template <typename T> constexpr T atan_product(T x, int k) { return k == 1 ? atan_term(x*x, k) : atan_term(x*x, k) * atan_product(x, k-1); } template <typename T> constexpr T atan_sum(T x, T sum, int n) { return sum + atan_product(x, n) == sum ? sum : atan_sum(x, sum + atan_product(x, n), n+1); } } template <typename T> constexpr T atan(T x) { return true ? x / (T{1} + x*x) * detail::atan_sum(x, T{1}, 1) : throw err::atan_runtime_error; } |

And `atan2`

follows from `atan`

relatively easily.

**It’s big, it’s heavy, it’s wood**

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`

. So I programmed it using my friend the Taylor series expansion.

Hm. It turns out that a naive implementation of `log`

is quite slow to converge. (Yes, what a surprise.) OK, a bit of searching yields fruit: “… use Newton’s method to invert the exponentional function … the iteration simplifies to:”

“… which has cubic convergence to ln(x).”

Result! I coded it up.

namespace detail { template <typename T> constexpr T log_iter(T x, T y) { return y + T{2} * (x - cx::exp(y)) / (x + cx::exp(y)); } template <typename T> constexpr T log(T x, T y) { return feq(y, log_iter(x, y)) ? y : log(x, log_iter(x, y)); } } |

It worked well. A week or so later, an article about C++11 constexpr log 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:

and

Coding these expressions to constrain the input proper to `detail::log`

seemed to solve the instability issues.

**Are we done yet?**

Now I was nearing the end of what I felt was a quorum of functions. To wit:

`sqrt`

,`cbrt`

,`hypot`

`sin`

,`cos`

,`tan`

`floor`

,`ceil`

,`trunc`

,`round`

`asin`

,`acos`

,`atan`

,`atan2`

`fmod`

,`remainder`

`exp`

,`log`

,`log2`

,`log10`

With a decent implementation of `exp`

and `log`

it was easy to add the hyperbolic trig functions and their inverses. The only obviously missing function now was `pow`

. Recall that I already had `ipow`

for raising a floating-point value to an integral power, but I had been puzzling over a more general version of `pow`

to allow raising to a floating-point power for a while. And then I suddenly realised the blindingly obvious:

And for the general case, `pow`

was solved.

template <typename T> constexpr T pow(T x, T y) { return true ? exp(log(x)*y) : throw err::pow_runtime_error; } |

I had now reached a satisfying conclusion to maths functions and I was getting fairly well-versed in writing `constexpr`

code. Pausing only to implement `erf`

(because why not?), I turned my attention now to other `constexpr`

matters.