[libc] Provide Arm Optimized Routines for the LLVM libc project.
[lldb.git] / libc / AOR_v20.02 / math / sincosf.h
1 /*
2  * Header for sinf, cosf and sincosf.
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 <stdint.h>
10 #include <math.h>
11 #include "math_config.h"
12
13 /* 2PI * 2^-64.  */
14 static const double pi63 = 0x1.921FB54442D18p-62;
15 /* PI / 4.  */
16 static const double pio4 = 0x1.921FB54442D18p-1;
17
18 /* The constants and polynomials for sine and cosine.  */
19 typedef struct
20 {
21   double sign[4];               /* Sign of sine in quadrants 0..3.  */
22   double hpi_inv;               /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS).  */
23   double hpi;                   /* PI / 2.  */
24   double c0, c1, c2, c3, c4;    /* Cosine polynomial.  */
25   double s1, s2, s3;            /* Sine polynomial.  */
26 } sincos_t;
27
28 /* Polynomial data (the cosine polynomial is negated in the 2nd entry).  */
29 extern const sincos_t __sincosf_table[2] HIDDEN;
30
31 /* Table with 4/PI to 192 bit precision.  */
32 extern const uint32_t __inv_pio4[] HIDDEN;
33
34 /* Top 12 bits of the float representation with the sign bit cleared.  */
35 static inline uint32_t
36 abstop12 (float x)
37 {
38   return (asuint (x) >> 20) & 0x7ff;
39 }
40
41 /* Compute the sine and cosine of inputs X and X2 (X squared), using the
42    polynomial P and store the results in SINP and COSP.  N is the quadrant,
43    if odd the cosine and sine polynomials are swapped.  */
44 static inline void
45 sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp,
46               float *cosp)
47 {
48   double x3, x4, x5, x6, s, c, c1, c2, s1;
49
50   x4 = x2 * x2;
51   x3 = x2 * x;
52   c2 = p->c3 + x2 * p->c4;
53   s1 = p->s2 + x2 * p->s3;
54
55   /* Swap sin/cos result based on quadrant.  */
56   float *tmp = (n & 1 ? cosp : sinp);
57   cosp = (n & 1 ? sinp : cosp);
58   sinp = tmp;
59
60   c1 = p->c0 + x2 * p->c1;
61   x5 = x3 * x2;
62   x6 = x4 * x2;
63
64   s = x + x3 * p->s1;
65   c = c1 + x4 * p->c2;
66
67   *sinp = s + x5 * s1;
68   *cosp = c + x6 * c2;
69 }
70
71 /* Return the sine of inputs X and X2 (X squared) using the polynomial P.
72    N is the quadrant, and if odd the cosine polynomial is used.  */
73 static inline float
74 sinf_poly (double x, double x2, const sincos_t *p, int n)
75 {
76   double x3, x4, x6, x7, s, c, c1, c2, s1;
77
78   if ((n & 1) == 0)
79     {
80       x3 = x * x2;
81       s1 = p->s2 + x2 * p->s3;
82
83       x7 = x3 * x2;
84       s = x + x3 * p->s1;
85
86       return s + x7 * s1;
87     }
88   else
89     {
90       x4 = x2 * x2;
91       c2 = p->c3 + x2 * p->c4;
92       c1 = p->c0 + x2 * p->c1;
93
94       x6 = x4 * x2;
95       c = c1 + x4 * p->c2;
96
97       return c + x6 * c2;
98     }
99 }
100
101 /* Fast range reduction using single multiply-subtract.  Return the modulo of
102    X as a value between -PI/4 and PI/4 and store the quadrant in NP.
103    The values for PI/2 and 2/PI are accessed via P.  Since PI/2 as a double
104    is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
105    the result is accurate for |X| <= 120.0.  */
106 static inline double
107 reduce_fast (double x, const sincos_t *p, int *np)
108 {
109   double r;
110 #if TOINT_INTRINSICS
111   /* Use fast round and lround instructions when available.  */
112   r = x * p->hpi_inv;
113   *np = converttoint (r);
114   return x - roundtoint (r) * p->hpi;
115 #else
116   /* Use scaled float to int conversion with explicit rounding.
117      hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
118      This avoids inaccuracies introduced by truncating negative values.  */
119   r = x * p->hpi_inv;
120   int n = ((int32_t)r + 0x800000) >> 24;
121   *np = n;
122   return x - n * p->hpi;
123 #endif
124 }
125
126 /* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
127    XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
128    Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
129    Reduction uses a table of 4/PI with 192 bits of precision.  A 32x96->128 bit
130    multiply computes the exact 2.62-bit fixed-point modulo.  Since the result
131    can have at most 29 leading zeros after the binary point, the double
132    precision result is accurate to 33 bits.  */
133 static inline double
134 reduce_large (uint32_t xi, int *np)
135 {
136   const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
137   int shift = (xi >> 23) & 7;
138   uint64_t n, res0, res1, res2;
139
140   xi = (xi & 0xffffff) | 0x800000;
141   xi <<= shift;
142
143   res0 = xi * arr[0];
144   res1 = (uint64_t)xi * arr[4];
145   res2 = (uint64_t)xi * arr[8];
146   res0 = (res2 >> 32) | (res0 << 32);
147   res0 += res1;
148
149   n = (res0 + (1ULL << 61)) >> 62;
150   res0 -= n << 62;
151   double x = (int64_t)res0;
152   *np = n;
153   return x * pi63;
154 }