libc/AOR_v20.02/math/v_cos.c

   1 /*
   2  * Double-precision vector cos 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 HalfPi v_f64 (0x1.921fb54442d18p+0)
  35 #define Pi1 v_f64 (0x1.921fb54442d18p+1)
  36 #define Pi2 v_f64 (0x1.1a62633145c06p-53)
  37 #define Pi3 v_f64 (0x1.c1cd129024e09p-106)
  38 #define Shift v_f64 (0x1.8p52)
  39 #define RangeVal v_f64 (0x1p23)
  40 #define AbsMask v_u64 (0x7fffffffffffffff)
  41
  42 VPCS_ATTR
  43 __attribute__ ((noinline)) static v_f64_t
  44 specialcase (v_f64_t x, v_f64_t y, v_u64_t cmp)
  45 {
  46   return v_call_f64 (cos, x, y, cmp);
  47 }
  48
  49 VPCS_ATTR
  50 v_f64_t
  51 V_NAME(cos) (v_f64_t x)
  52 {
  53   v_f64_t n, r, r2, y;
  54   v_u64_t odd, cmp;
  55
  56   r = v_as_f64_u64 (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/2)/pi) - 0.5.  */
  60   n = v_fma_f64 (InvPi, r + HalfPi, Shift);
  61   odd = v_as_u64_f64 (n) << 63;
  62   n -= Shift;
  63   n -= v_f64 (0.5);
  64
  65   /* r = |x| - n*pi  (range reduction into -pi/2 .. pi/2).  */
  66   r = v_fma_f64 (-Pi1, n, r);
  67   r = v_fma_f64 (-Pi2, n, r);
  68   r = v_fma_f64 (-Pi3, n, r);
  69
  70   /* sin(r) poly approx.  */
  71   r2 = r * r;
  72   y = v_fma_f64 (C7, r2, C6);
  73   y = v_fma_f64 (y, r2, C5);
  74   y = v_fma_f64 (y, r2, C4);
  75   y = v_fma_f64 (y, r2, C3);
  76   y = v_fma_f64 (y, r2, C2);
  77   y = v_fma_f64 (y, r2, C1);
  78   y = v_fma_f64 (y * r2, r, r);
  79
  80   /* sign.  */
  81   y = v_as_f64_u64 (v_as_u64_f64 (y) ^ odd);
  82
  83   if (unlikely (v_any_u64 (cmp)))
  84     return specialcase (x, y, cmp);
  85   return y;
  86 }
  87 VPCS_ALIAS
  88 #endif