[PATCH v4 2/4] aarch64: Add vector implementations of sin routines

[Szabolcs Nagy]

The 06/28/2023 12:19, Joe Ramsay via Libc-alpha wrote:
> +++ b/sysdeps/aarch64/fpu/sin_advsimd.c
> @@ -0,0 +1,106 @@
> +/* Double-precision vector (Advanced SIMD) sin function.
...
> +float64x2_t VPCS_ATTR V_NAME_D1 (sin) (float64x2_t x)
> +{
> +  const struct data *d = ptr_barrier (&data);
> +  float64x2_t n, r, r2, r3, r4, y, t1, t2, t3;
> +  uint64x2_t odd, cmp, eqz;
> +
> +#if WANT_SIMD_EXCEPT
> +  /* Detect |x| <= TinyBound or |x| >= RangeVal. If fenv exceptions are to be
> +     triggered correctly, set any special lanes to 1 (which is neutral w.r.t.
> +     fenv). These lanes will be fixed by special-case handler later.  */
> +  uint64x2_t ir = vreinterpretq_u64_f64 (vabsq_f64 (x));
> +  cmp = vcgeq_u64 (vsubq_u64 (ir, TinyBound), Thresh);
> +  r = vbslq_f64 (cmp, vreinterpretq_f64_u64 (cmp), x);
> +#else
> +  r = x;
> +  cmp = vcageq_f64 (d->range_val, x);
> +  cmp = vceqzq_u64 (cmp); /* cmp = ~cmp.  */
> +#endif
> +  eqz = vceqzq_f64 (x);
> +
> +  /* n = rint(|x|/pi).  */
> +  n = vfmaq_f64 (d->shift, d->inv_pi, r);
> +  odd = vshlq_n_u64 (vreinterpretq_u64_f64 (n), 63);
> +  n = vsubq_f64 (n, d->shift);
> +
> +  /* r = |x| - n*pi  (range reduction into -pi/2 .. pi/2).  */
> +  r = vfmsq_f64 (r, d->pi_1, n);
> +  r = vfmsq_f64 (r, d->pi_2, n);
> +  r = vfmsq_f64 (r, d->pi_3, n);
> +
> +  /* sin(r) poly approx.  */
> +  r2 = vmulq_f64 (r, r);
> +  r3 = vmulq_f64 (r2, r);
> +  r4 = vmulq_f64 (r2, r2);
> +
> +  t1 = vfmaq_f64 (C (4), C (5), r2);
> +  t2 = vfmaq_f64 (C (2), C (3), r2);
> +  t3 = vfmaq_f64 (C (0), C (1), r2);
> +
> +  y = vfmaq_f64 (t1, C (6), r4);
> +  y = vfmaq_f64 (t2, y, r4);
> +  y = vfmaq_f64 (t3, y, r4);
> +  y = vfmaq_f64 (r, y, r3);
> +
> +  /* Sign of 0 is discarded by polynomial, so copy it back here.  */
> +  if (__glibc_unlikely (v_any_u64 (eqz)))
> +    y = vbslq_f64 (eqz, x, y);

this check is for dealing with the sin(-0.0) case.

but -ffast-math can already break the sign of 0 and libmvec
symbols are supposed to be used with -ffast-math (at least via
gcc auto vectorizer).

so do we need to provide quality guarantees beyond -ffast-math
in libmvec? or can we ignore math tests that check for the sign
of 0 and change the code accordingly.

the wiki does not seem to cover this under "C99 compliance":
https://sourceware.org/glibc/wiki/libmvec


> +
> +  if (__glibc_unlikely (v_any_u64 (cmp)))
> +    return special_case (x, y, odd, cmp);
> +  return vreinterpretq_f64_u64 (veorq_u64 (vreinterpretq_u64_f64 (y), odd));
> +}