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:

[latex]arctan(x) = arcsin(\frac{x}{\sqrt{x^2+1}})[/latex]

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

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

[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]

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

namespace detail
{
  template 
  constexpr T atan_term(T x2, int k)
  {
    return (T{2}*static_cast(k)*x2)
      / ((T{2}*static_cast(k)+T{1}) * (T{1}+x2));
  }
  template 
  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 
  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 
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:”

[latex]y_{n+1} = y_n + 2 \frac{x – e^{y_n}}{x + e^{y_n}}[/latex]

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

Result! I coded it up.

namespace detail
{
  template 
  constexpr T log_iter(T x, T y)
  {
    return y + T{2} * (x - cx::exp(y)) / (x + cx::exp(y));
  }
  template 
  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:

[latex]ln(x) = ln(\frac{x}{e}) + 1[/latex]
and
[latex]ln(x) = ln(e.x) – 1[/latex]

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:

[latex]x^y = (e^{ln(x)})^y = e^{ln(x).y}[/latex]

And for the general case, pow was solved.

template 
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.