## More constexpr floating-point computation

(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:

$arctan(x) = arcsin(\frac{x}{\sqrt{x^2+1}})$

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

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

$arctan(z) = \frac{z}{1+z^2} \sum_{n=0}^{\infty} \prod_{k=1}^n \frac{2kz^2}{(2k+1)(1+z^2)}$

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:”

$y_{n+1} = y_n + 2 \frac{x - e^{y_n}}{x + e^{y_n}}$

“… 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:

$ln(x) = ln(\frac{x}{e}) + 1$
and
$ln(x) = ln(e.x) - 1$

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:

$x^y = (e^{ln(x)})^y = e^{ln(x).y}$

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.