xref: /freebsd/contrib/arm-optimized-routines/math/aarch64/sincospi_4u.c (revision f3087bef11543b42e0d69b708f367097a4118d24)
1 /*
2  * Double-precision scalar sincospi function.
3  *
4  * Copyright (c) 2024, Arm Limited.
5  * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6  */
7 
8 #include "mathlib.h"
9 #include "math_config.h"
10 #include "test_sig.h"
11 #include "test_defs.h"
12 #include "poly_scalar_f64.h"
13 
14 /* Taylor series coefficents for sin(pi * x).
15    C2 coefficient (orginally ~=5.16771278) has been split into two parts:
16    C2_hi = 4, C2_lo = C2 - C2_hi (~=1.16771278)
17    This change in magnitude reduces floating point rounding errors.
18    C2_hi is then reintroduced after the polynomial approxmation.  */
19 const static struct sincospi_data
20 {
21   double poly[10];
22 } sincospi_data = {
23   /* Taylor series coefficents for sin(pi * x).  */
24   .poly = { 0x1.921fb54442d184p1, -0x1.2aef39896f94bp0, 0x1.466bc6775ab16p1,
25 	    -0x1.32d2cce62dc33p-1, 0x1.507834891188ep-4, -0x1.e30750a28c88ep-8,
26 	    0x1.e8f48308acda4p-12, -0x1.6fc0032b3c29fp-16,
27 	    0x1.af86ae521260bp-21, -0x1.012a9870eeb7dp-25 },
28 };
29 
30 /* Top 12 bits of a double (sign and exponent bits).  */
31 static inline uint64_t
abstop12(double x)32 abstop12 (double x)
33 {
34   return (asuint64 (x) >> 52) & 0x7ff;
35 }
36 
37 /* Triages special cases into 4 categories:
38      -1 or +1 if iy represents half an integer
39        -1 if round(y) is odd.
40        +1 if round(y) is even.
41      -2 or +2 if iy represents and integer.
42        -2 if iy is odd.
43        +2 if iy is even.
44    The argument is the bit representation of a positive non-zero
45    finite floating-point value which is either a half or an integer.  */
46 static inline int
checkint(uint64_t iy)47 checkint (uint64_t iy)
48 {
49   int e = iy >> 52;
50   if (e > 0x3ff + 52)
51     return 2;
52   if (iy & ((1ULL << (0x3ff + 52 - e)) - 1))
53     {
54       if ((iy - 1) & 2)
55 	return -1;
56       else
57 	return 1;
58     }
59   if (iy & (1 << (0x3ff + 52 - e)))
60     return -2;
61   return 2;
62 }
63 
64 /* Approximation for scalar double-precision sincospi(x).
65    Maximum error for sin: 3.46 ULP:
66       sincospif_sin(0x1.3d8a067cd8961p+14) got 0x1.ffe609a279008p-1 want
67    0x1.ffe609a27900cp-1.
68    Maximum error for cos: 3.66 ULP:
69       sincospif_cos(0x1.a0ec6997557eep-24) got 0x1.ffffffffffe59p-1 want
70    0x1.ffffffffffe5dp-1.  */
71 void
arm_math_sincospi(double x,double * out_sin,double * out_cos)72 arm_math_sincospi (double x, double *out_sin, double *out_cos)
73 {
74   const struct sincospi_data *d = ptr_barrier (&sincospi_data);
75   uint64_t sign = asuint64 (x) & 0x8000000000000000;
76 
77   if (likely (abstop12 (x) < abstop12 (0x1p51)))
78     {
79       /* ax = |x| - n (range reduction into -1/2 .. 1/2).  */
80       double ar_s = x - rint (x);
81 
82       /* We know that cospi(x) = sinpi(0.5 - x)
83 	 range reduction and offset into sinpi range -1/2 .. 1/2
84 	 ax = 0.5 - |x - rint(x)|.  */
85       double ar_c = 0.5 - fabs (ar_s);
86 
87       /* ss = sin(pi * ax).  */
88       double ar2_s = ar_s * ar_s;
89       double ar2_c = ar_c * ar_c;
90       double ar4_s = ar2_s * ar2_s;
91       double ar4_c = ar2_c * ar2_c;
92 
93       uint64_t cc_sign = ((uint64_t) llrint (x)) << 63;
94       uint64_t ss_sign = cc_sign;
95       if (ar_s == 0)
96 	ss_sign = sign;
97 
98       double ss = pw_horner_9_f64 (ar2_s, ar4_s, d->poly);
99       double cc = pw_horner_9_f64 (ar2_c, ar4_c, d->poly);
100 
101       /* As all values are reduced to -1/2 .. 1/2, the result of cos(x)
102 	 always be positive, therefore, the sign must be introduced
103 	 based upon if x rounds to odd or even. For sin(x) the sign is
104 	 copied from x.  */
105       *out_sin
106 	  = asdouble (asuint64 (fma (-4 * ar2_s, ar_s, ss * ar_s)) ^ ss_sign);
107       *out_cos
108 	  = asdouble (asuint64 (fma (-4 * ar2_c, ar_c, cc * ar_c)) ^ cc_sign);
109     }
110   else
111     {
112       /* When abs(x) > 0x1p51, the x will be either
113 	    - Half integer (relevant if abs(x) in [0x1p51, 0x1p52])
114 	    - Odd integer  (relevant if abs(x) in [0x1p52, 0x1p53])
115 	    - Even integer (relevant if abs(x) in [0x1p53, inf])
116 	    - Inf or NaN.  */
117       if (abstop12 (x) >= 0x7ff)
118 	{
119 	  double inv_result = __math_invalid (x);
120 	  *out_sin = inv_result;
121 	  *out_cos = inv_result;
122 	  return;
123 	}
124       else
125 	{
126 	  uint64_t ax = asuint64 (x) & 0x7fffffffffffffff;
127 	  int m = checkint (ax);
128 	  /* The case where ax is half integer.  */
129 	  if (m & 1)
130 	    {
131 	      *out_sin = sign ? -m : m;
132 	      *out_cos = 0;
133 	      return;
134 	    }
135 	  /* The case where ax is integer.  */
136 	  else
137 	    {
138 	      *out_sin = asdouble (sign);
139 	      *out_cos = m >> 1;
140 	      return;
141 	    }
142 	}
143     }
144 }
145 
146 #if WANT_TRIGPI_TESTS
147 TEST_DISABLE_FENV (arm_math_sincospi_sin)
148 TEST_DISABLE_FENV (arm_math_sincospi_cos)
149 TEST_ULP (arm_math_sincospi_sin, 2.96)
150 TEST_ULP (arm_math_sincospi_cos, 3.16)
151 #  define SINCOS_INTERVAL(lo, hi, n)                                          \
152     TEST_SYM_INTERVAL (arm_math_sincospi_sin, lo, hi, n)                      \
153     TEST_SYM_INTERVAL (arm_math_sincospi_cos, lo, hi, n)
154 SINCOS_INTERVAL (0, 0x1p-63, 10000)
155 SINCOS_INTERVAL (0x1p-63, 0.5, 50000)
156 SINCOS_INTERVAL (0.5, 0x1p51, 50000)
157 SINCOS_INTERVAL (0x1p51, inf, 10000)
158 #endif
159