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