1  /* Double-precision vector (SVE) sin function.
       2  
       3     Copyright (C) 2023 Free Software Foundation, Inc.
       4     This file is part of the GNU C Library.
       5  
       6     The GNU C Library is free software; you can redistribute it and/or
       7     modify it under the terms of the GNU Lesser General Public
       8     License as published by the Free Software Foundation; either
       9     version 2.1 of the License, or (at your option) any later version.
      10  
      11     The GNU C Library is distributed in the hope that it will be useful,
      12     but WITHOUT ANY WARRANTY; without even the implied warranty of
      13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      14     Lesser General Public License for more details.
      15  
      16     You should have received a copy of the GNU Lesser General Public
      17     License along with the GNU C Library; if not, see
      18     <https://www.gnu.org/licenses/>.  */
      19  
      20  #include "sv_math.h"
      21  
      22  static const struct data
      23  {
      24    double inv_pi, half_pi, inv_pi_over_2, pi_over_2_1, pi_over_2_2, pi_over_2_3,
      25        shift;
      26  } data = {
      27    /* Polynomial coefficients are hard-wired in the FTMAD instruction.  */
      28    .inv_pi = 0x1.45f306dc9c883p-2,
      29    .half_pi = 0x1.921fb54442d18p+0,
      30    .inv_pi_over_2 = 0x1.45f306dc9c882p-1,
      31    .pi_over_2_1 = 0x1.921fb50000000p+0,
      32    .pi_over_2_2 = 0x1.110b460000000p-26,
      33    .pi_over_2_3 = 0x1.1a62633145c07p-54,
      34    .shift = 0x1.8p52
      35  };
      36  
      37  #define RangeVal 0x4160000000000000 /* asuint64 (0x1p23).  */
      38  
      39  static svfloat64_t NOINLINE
      40  special_case (svfloat64_t x, svfloat64_t y, svbool_t cmp)
      41  {
      42    return sv_call_f64 (sin, x, y, cmp);
      43  }
      44  
      45  /* A fast SVE implementation of sin based on trigonometric
      46     instructions (FTMAD, FTSSEL, FTSMUL).
      47     Maximum observed error in 2.52 ULP:
      48     SV_NAME_D1 (sin)(0x1.2d2b00df69661p+19) got 0x1.10ace8f3e786bp-40
      49  					  want 0x1.10ace8f3e7868p-40.  */
      50  svfloat64_t SV_NAME_D1 (sin) (svfloat64_t x, const svbool_t pg)
      51  {
      52    const struct data *d = ptr_barrier (&data);
      53  
      54    svfloat64_t r = svabs_f64_x (pg, x);
      55    svuint64_t sign
      56        = sveor_u64_x (pg, svreinterpret_u64_f64 (x), svreinterpret_u64_f64 (r));
      57    svbool_t cmp = svcmpge_n_u64 (pg, svreinterpret_u64_f64 (r), RangeVal);
      58  
      59    /* Load first two pio2-related constants to one vector.  */
      60    svfloat64_t invpio2_and_pio2_1
      61        = svld1rq_f64 (svptrue_b64 (), &d->inv_pi_over_2);
      62  
      63    /* n = rint(|x|/(pi/2)).  */
      64    svfloat64_t q = svmla_lane_f64 (sv_f64 (d->shift), r, invpio2_and_pio2_1, 0);
      65    svfloat64_t n = svsub_n_f64_x (pg, q, d->shift);
      66  
      67    /* r = |x| - n*(pi/2)  (range reduction into -pi/4 .. pi/4).  */
      68    r = svmls_lane_f64 (r, n, invpio2_and_pio2_1, 1);
      69    r = svmls_n_f64_x (pg, r, n, d->pi_over_2_2);
      70    r = svmls_n_f64_x (pg, r, n, d->pi_over_2_3);
      71  
      72    /* Final multiplicative factor: 1.0 or x depending on bit #0 of q.  */
      73    svfloat64_t f = svtssel_f64 (r, svreinterpret_u64_f64 (q));
      74  
      75    /* sin(r) poly approx.  */
      76    svfloat64_t r2 = svtsmul_f64 (r, svreinterpret_u64_f64 (q));
      77    svfloat64_t y = sv_f64 (0.0);
      78    y = svtmad_f64 (y, r2, 7);
      79    y = svtmad_f64 (y, r2, 6);
      80    y = svtmad_f64 (y, r2, 5);
      81    y = svtmad_f64 (y, r2, 4);
      82    y = svtmad_f64 (y, r2, 3);
      83    y = svtmad_f64 (y, r2, 2);
      84    y = svtmad_f64 (y, r2, 1);
      85    y = svtmad_f64 (y, r2, 0);
      86  
      87    /* Apply factor.  */
      88    y = svmul_f64_x (pg, f, y);
      89  
      90    /* sign = y^sign.  */
      91    y = svreinterpret_f64_u64 (
      92        sveor_u64_x (pg, svreinterpret_u64_f64 (y), sign));
      93  
      94    if (__glibc_unlikely (svptest_any (pg, cmp)))
      95      return special_case (x, y, cmp);
      96    return y;
      97  }