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