sin/cos with constexpr

[Marcel Mueller]

Is there some work around to get a constexpr version of trigonometric
math functions like sin, cos?

I want to store precalculated twiddle factors for GPU accelerated FFT
calculation in the read only data segment. Basically they are all
constants of the form
static const float twiddles[] =
{ (float)cos(M_2PI*0/256), (float)sin(M_2PI*0/256),
(float)cos(M_2PI*1/256), (float)sin(M_2PI*1/256),
...
(float)cos(M_2PI*127/256), (float)sin(M_2PI*127/256),
};
Of course, this will not work since sin/cos are not constexpr.

My current work around is to create the twiddle table with the GPU
assembler (it can handle math functions in compile time constants) and
read it with the CPU just to feed them to the shader code.

Is there some native C++11 work around?


Marcel

[mark]

On 2015-09-13 21:24, Marcel Mueller wrote:
> Is there some work around to get a constexpr version of trigonometric
> math functions like sin, cos?

I don't think there is a 'sane' possibility. But, you if you don't care
about sanity, you can evaluate a taylor series as constexpr:
http://prosepoetrycode.potterpcs.net/2015/07/more-fun-with-constexpr/

GGC >= 4.7 has the cmath stuff as constexpr, but it's not even required
by C++14.

[Öö Tiib]

No. After we, programmers realized how powerful tool is constexpr
some compiler vendors made some things that weren't constexpr to be
constexpr (like 'std::fmod'). Some other compiler vendors who were
late with implementing constexpr whatsoever (only achieved non-fully
in VS 2015) used their power to forbid that.

Their excuse was usual pathetic nonsense about breaking legacy.
Those functions have side effect of setting 'errno' on case of some
domain errors. So one can make code that "breaks" because 'sin' or
'cos' is constexpr.

Long story short ... unless you don't mind to use table lookup
(your cosine and sine tables seem to be 128 entries long) there
are no standard ways. You either have to generate your twiddle
factor table with code generators or to use non-standard compiler
intrinsics that work compile time.

[K. Frank]

Hi Marcel!

On Sunday, September 13, 2015 at 3:25:10 PM UTC-4, Marcel Mueller wrote:
> Is there some work around to get a constexpr version of trigonometric
> math functions like sin, cos?

> ...

>
> Is there some native C++11 work around?
>
> Marcel

Just a side note -- I don't have an answer to your
work-around question.

Floating-point constexpr raises issues with
cross-compilation, especially with transcendental
functions.

This is because, respecting (legitimate) differences
across hardware, c++ does not require floating-point
computation to be exactly portable.

So there is really no way to ensure that, for example,
"constexpr cos (1.23)" is equal to "cos (1.23)" as
evaluated at run time. (Although they should be equal
up to some reasonable round-off error, where "reasonable"
is a quality-of-implementation issue).

This isn't really a show-stopper, but it does let one
more potential floating-point "gotcha" creep in that
we would then have to avoid tripping over.


Happy Floating-Point hacking!


K. Frank

[Marcel Mueller]

On 13.09.15 23.54, K. Frank wrote:
> Floating-point constexpr raises issues with
> cross-compilation, especially with transcendental
> functions.

> This isn't really a show-stopper, but it does let one
> more potential floating-point "gotcha" creep in that
> we would then have to avoid tripping over.

Since I need to crop to 32 bit IEEE float anyway to get GPU support,
this is no big deal. Probably any real life platform should have enough
accuracy in the double calculation to produce reasonable exact values in
single precision, even when cross compiling.

> Happy Floating-Point hacking!

Well, most of the work is not related to numerics but to cache
efficiency and things like that.


Marcel

[Marcel Mueller]

Hmm, I must have done something wrong...
gcc 4.8.2
I got "initializer element is not a constant expression".

I will check this later.
Since the code is platform dependent anyway I could live with the
restriction to gcc at this place.

Oops, I got it. This file is named .c - shame on me (not my code).


Marcel

[Christian Gollwitzer]

Am 13.09.15 um 22:13 schrieb mark:

> On 2015-09-13 21:24, Marcel Mueller wrote:
>> Is there some work around to get a constexpr version of trigonometric
>> math functions like sin, cos?
>
> I don't think there is a 'sane' possibility. But, you if you don't care
> about sanity, you can evaluate a taylor series as constexpr:
> http://prosepoetrycode.potterpcs.net/2015/07/more-fun-with-constexpr/
>

A comment to this: it's fun to see how templates can be abused to
compute the cosine series. However you should be aware that there might
be precision problems. Serious implementations of trigonometric
functions can be precise to the last digit in floating point, but they
don't use the Taylor expansion, instead a minmax approximation is used.
cephes on netlib contains many special functions as pure C code, which
could probably be converted to constexpr.

http://www.netlib.org/cephes/cmath.tgz

sin and cos in double precision are here:

https://github.com/jeremybarnes/cephes/blob/master/cmath/sin.c


Christian

[Nobody]

On Sun, 13 Sep 2015 22:13:21 +0200, mark wrote:

>> Is there some work around to get a constexpr version of trigonometric
>> math functions like sin, cos?
>
> I don't think there is a 'sane' possibility. But, you if you don't care
> about sanity, you can evaluate a taylor series as constexpr:

Another option is a recursive approach based upon the sum/difference
formulae and half-angle formulae.

Corollary: for angles of the form pi*i/(2^n) where i and n are
integers (aka Binary Angular Measurements or BAMs for short), the sine,
cosine and tangent of such angles aren't transcendental but constructible,
i.e. can be evaluated using only +, -, *, / and sqrt(). Thus it's possible
to determine the nearest representable value using only rational
arithmetic (so no "table-maker's dilemma").

[David Brown]

On 13/09/15 23:54, K. Frank wrote:
> Hi Marcel!
>
> On Sunday, September 13, 2015 at 3:25:10 PM UTC-4, Marcel Mueller wrote:
>> Is there some work around to get a constexpr version of trigonometric
>> math functions like sin, cos?
>> ...
>>
>> Is there some native C++11 work around?
>>
>> Marcel
>
> Just a side note -- I don't have an answer to your
> work-around question.
>
> Floating-point constexpr raises issues with
> cross-compilation, especially with transcendental
> functions.
>
> This is because, respecting (legitimate) differences
> across hardware, c++ does not require floating-point
> computation to be exactly portable.
>
> So there is really no way to ensure that, for example,
> "constexpr cos (1.23)" is equal to "cos (1.23)" as
> evaluated at run time. (Although they should be equal
> up to some reasonable round-off error, where "reasonable"
> is a quality-of-implementation issue).

There /is/ a way to handle this - the compiler can have libraries that
emulate the target details. gcc (since 4.5) has used a library to do
this, and will happily and confidently handle constant compile-time
evaluation of floating point operations, and including transcendental
ones, including cross-compilation.

To help matters, if you are looking for "fast and correct for normal
code" rather than "follow IEEE to the letter and the last bit", you can
use the "-ffast-math" flag to tell gcc not to worry /too/ much about
such details. Other compilers may have something similar.

[Christian Gollwitzer]

Am 14.09.15 um 09:14 schrieb Nobody:

> On Sun, 13 Sep 2015 22:13:21 +0200, mark wrote:
>
>>> Is there some work around to get a constexpr version of trigonometric
>>> math functions like sin, cos?
>>
>> I don't think there is a 'sane' possibility. But, you if you don't care
>> about sanity, you can evaluate a taylor series as constexpr:
>
> Another option is a recursive approach based upon the sum/difference
> formulae and half-angle formulae.

which is usually called the CORDIC algorithm. It is simple enough to be
implemented in hardware (pocket calculators), but inferior in speed if a
fast multiplier is available. Of course, for compile-time copmutation
speed is probably the least important

https://en.wikipedia.org/wiki/CORDIC

Christian

[Nobody]

On Mon, 14 Sep 2015 21:20:52 +0200, Christian Gollwitzer wrote:

>> Another option is a recursive approach based upon the sum/difference
>> formulae and half-angle formulae.
>
> which is usually called the CORDIC algorithm.

They're different algorithms. CORDIC is based upon adding or subtracting
angles of the form arctan(1/(2^n)); the half-angle approach uses angles of
the form pi/(2^n).

CORDIC is more efficient on simple hardware, as it only requires shift,
add/subtract, and a lookup table.

Half-angle doesn't require the lookup table, but it does require square
roots. But std::sqrt() doesn't appear to be constexpr either, so you'd
need to implement that as well.

The main advantage of half-angle over CORDIC or a Taylor series is that it
doesn't require a set of constants whose number and accuracy needs to
change according to the desired accuracy of the result. Although that
probably doesn't matter here.

[Nobody]

On Sun, 13 Sep 2015 22:13:21 +0200, mark wrote:

> On 2015-09-13 21:24, Marcel Mueller wrote:
>> Is there some work around to get a constexpr version of trigonometric
>> math functions like sin, cos?
>
> I don't think there is a 'sane' possibility. But, you if you don't care
> about sanity, you can evaluate a taylor series as constexpr:
> http://prosepoetrycode.potterpcs.net/2015/07/more-fun-with-constexpr/

FWIW, here's a complete implementation of such:

#include <cmath>

namespace {
static const int terms = 20;

constexpr double alternate(double x2, int i, int k, double xn, long long nf)
{
return (i > terms) ? 0.0 :
(k*xn/nf + alternate(x2, i+2, -k, xn*x2, nf*(i+1)*(i+2)));
}
}

namespace taylor {
constexpr double sin(double x)
{
return alternate(x*x, 1, 1, x, 1);
}

constexpr double cos(double x)
{
return 1.0 + alternate(x*x, 2, -1, x*x, 2);
}
}

It doesn't always agree with libc's sin/cos, but any difference appears to
be limited to the last bit. sin() differs in ~15% of cases, cos() in ~30%.