libc/AOR_v20.02/math/v_sin.c

   1 /*
   2  * Double-precision vector sin function.
   3  *
   4  * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   5  * See https://llvm.org/LICENSE.txt for license information.
   6  * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   7  */
   8
   9 #include "mathlib.h"
  10 #include "v_math.h"
  11 #if V_SUPPORTED
  12
  13 static const double Poly[] = {
  14 /* worst-case error is 3.5 ulp.
  15    abs error: 0x1.be222a58p-53 in [-pi/2, pi/2].  */
  16 -0x1.9f4a9c8b21dc9p-41,
  17  0x1.60e88a10163f2p-33,
  18 -0x1.ae6361b7254e7p-26,
  19  0x1.71de382e8d62bp-19,
  20 -0x1.a01a019aeb4ffp-13,
  21  0x1.111111110b25ep-7,
  22 -0x1.55555555554c3p-3,
  23 };
  24
  25 #define C7 v_f64 (Poly[0])
  26 #define C6 v_f64 (Poly[1])
  27 #define C5 v_f64 (Poly[2])
  28 #define C4 v_f64 (Poly[3])
  29 #define C3 v_f64 (Poly[4])
  30 #define C2 v_f64 (Poly[5])
  31 #define C1 v_f64 (Poly[6])
  32
  33 #define InvPi v_f64 (0x1.45f306dc9c883p-2)
  34 #define Pi1 v_f64 (0x1.921fb54442d18p+1)
  35 #define Pi2 v_f64 (0x1.1a62633145c06p-53)
  36 #define Pi3 v_f64 (0x1.c1cd129024e09p-106)
  37 #define Shift v_f64 (0x1.8p52)
  38 #define RangeVal v_f64 (0x1p23)
  39 #define AbsMask v_u64 (0x7fffffffffffffff)
  40
  41 VPCS_ATTR
  42 __attribute__ ((noinline)) static v_f64_t
  43 specialcase (v_f64_t x, v_f64_t y, v_u64_t cmp)
  44 {
  45   return v_call_f64 (sin, x, y, cmp);
  46 }
  47
  48 VPCS_ATTR
  49 v_f64_t
  50 V_NAME(sin) (v_f64_t x)
  51 {
  52   v_f64_t n, r, r2, y;
  53   v_u64_t sign, odd, cmp;
  54
  55   r = v_as_f64_u64 (v_as_u64_f64 (x) & AbsMask);
  56   sign = v_as_u64_f64 (x) & ~AbsMask;
  57   cmp = v_cond_u64 (v_as_u64_f64 (r) >= v_as_u64_f64 (RangeVal));
  58
  59   /* n = rint(|x|/pi).  */
  60   n = v_fma_f64 (InvPi, r, Shift);
  61   odd = v_as_u64_f64 (n) << 63;
  62   n -= Shift;
  63
  64   /* r = |x| - n*pi  (range reduction into -pi/2 .. pi/2).  */
  65   r = v_fma_f64 (-Pi1, n, r);
  66   r = v_fma_f64 (-Pi2, n, r);
  67   r = v_fma_f64 (-Pi3, n, r);
  68
  69   /* sin(r) poly approx.  */
  70   r2 = r * r;
  71   y = v_fma_f64 (C7, r2, C6);
  72   y = v_fma_f64 (y, r2, C5);
  73   y = v_fma_f64 (y, r2, C4);
  74   y = v_fma_f64 (y, r2, C3);
  75   y = v_fma_f64 (y, r2, C2);
  76   y = v_fma_f64 (y, r2, C1);
  77   y = v_fma_f64 (y * r2, r, r);
  78
  79   /* sign.  */
  80   y = v_as_f64_u64 (v_as_u64_f64 (y) ^ sign ^ odd);
  81
  82   if (unlikely (v_any_u64 (cmp)))
  83     return specialcase (x, y, cmp);
  84   return y;
  85 }
  86 VPCS_ALIAS
  87 #endif