contrib/arm-optimized-routines/math/v_sin.c

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