PPC64 libmvec sin, cos, sincos

PPC64 libmvec sin, cos, sincos

[GT]

Here is a suggestion to organize the implementation of the three functions
a little differently than was done for x86_64.

Because of the trigonometric identity cos (x) = sin( x + PI/2), selecting
either the sine or cosine polynomial series approximation should suffice to
compute all three functions.

Defining the series in a static inline function will let us call it with
input argument x or x +/- PI/2 for sine and cosine (depending on which of the
two series was chosen). For sincos, the function will be called twice: once
to obtain the cosines of the inputs and a second time to get the sines.

Doing that will allow us to maintain a single file which performs evaluations
for the three functions.

Any reason not to implement in this style?

Thanks.
Bert.

Re: PPC64 libmvec sin, cos, sincos

[Szabolcs Nagy]

On 03/03/2019 17:21, GT wrote:
> Here is a suggestion to organize the implementation of the three functions
> a little differently than was done for x86_64.
> 
> Because of the trigonometric identity cos (x) = sin( x + PI/2), selecting
> either the sine or cosine polynomial series approximation should suffice to
> compute all three functions.
> 
> Defining the series in a static inline function will let us call it with
> input argument x or x +/- PI/2 for sine and cosine (depending on which of the
> two series was chosen). For sincos, the function will be called twice: once
> to obtain the cosines of the inputs and a second time to get the sines.
> 
> Doing that will allow us to maintain a single file which performs evaluations
> for the three functions.
> 
> Any reason not to implement in this style?

i think it is better to compute

n = rint((x+pi/2)/pi) - 0.5
r = x - n*pi
sin(r)

(current x86_64 vector cos) instead of

n = rint((x+pi/2)/pi)
r = (x+pi/2) - n*pi
sin(r)

since -n*pi can be done with extra precision
(e.g. -n*pi1-n*pi2-n*pi3) but if you pass x+pi/2
down to a common implementation that may have a
large error (e.g. for x close to 0).

so the arg reduction may not be easy to share.
(and the final sign fixup)
but the polynomial can be the same i guess.