1 //===-- Collection of utils for cosf/sinf/sincosf ---------------*- C++ -*-===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8
9 #ifndef LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H
10 #define LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H
11
12 #include "math_utils.h"
13
14 #include <stdint.h>
15
16 namespace __llvm_libc {
17
18 // 2PI * 2^-64.
19 static const double pi63 = as_double(0x3c1921fb54442d18);
20 // PI / 4.
21 static const double pio4 = as_double(0x3fe921fb54442d18);
22
23 // The constants and polynomials for sine and cosine.
24 typedef struct {
25 double sign[4]; // Sign of sine in quadrants 0..3.
26 double hpi_inv; // 2 / PI ( * 2^24 ).
27 double hpi; // PI / 2.
28 double c0, c1, c2, c3, c4; // Cosine polynomial.
29 double s1, s2, s3; // Sine polynomial.
30 } sincos_t;
31
32 // Polynomial data (the cosine polynomial is negated in the 2nd entry).
33 extern const sincos_t __sincosf_table[2];
34
35 // Table with 4/PI to 192 bit precision.
36 extern const uint32_t __inv_pio4[];
37
38 // Top 12 bits of the float representation with the sign bit cleared.
abstop12(float x)39 static inline uint32_t abstop12(float x) {
40 return (as_uint32_bits(x) >> 20) & 0x7ff;
41 }
42
43 // Compute the sine and cosine of inputs X and X2 (X squared), using the
44 // polynomial P and store the results in SINP and COSP. N is the quadrant,
45 // if odd the cosine and sine polynomials are swapped.
sincosf_poly(double x,double x2,const sincos_t * p,int n,float * sinp,float * cosp)46 static inline void sincosf_poly(double x, double x2, const sincos_t *p, int n,
47 float *sinp, float *cosp) {
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.
sinf_poly(double x,double x2,const sincos_t * p,int n)73 static inline float sinf_poly(double x, double x2, const sincos_t *p, int n) {
74 double x3, x4, x6, x7, s, c, c1, c2, s1;
75
76 if ((n & 1) == 0) {
77 x3 = x * x2;
78 s1 = p->s2 + x2 * p->s3;
79
80 x7 = x3 * x2;
81 s = x + x3 * p->s1;
82
83 return s + x7 * s1;
84 } else {
85 x4 = x2 * x2;
86 c2 = p->c3 + x2 * p->c4;
87 c1 = p->c0 + x2 * p->c1;
88
89 x6 = x4 * x2;
90 c = c1 + x4 * p->c2;
91
92 return c + x6 * c2;
93 }
94 }
95
96 // Fast range reduction using single multiply-subtract. Return the modulo of
97 // X as a value between -PI/4 and PI/4 and store the quadrant in NP.
98 // The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
99 // is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
100 // the result is accurate for |X| <= 120.0.
reduce_fast(double x,const sincos_t * p,int * np)101 static inline double reduce_fast(double x, const sincos_t *p, int *np) {
102 double r;
103 // Use scaled float to int conversion with explicit rounding.
104 // hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
105 // This avoids inaccuracies introduced by truncating negative values.
106 r = x * p->hpi_inv;
107 int n = ((int32_t)r + 0x800000) >> 24;
108 *np = n;
109 return x - n * p->hpi;
110 }
111
112 // Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
113 // XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
114 // Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
115 // Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit
116 // multiply computes the exact 2.62-bit fixed-point modulo. Since the result
117 // can have at most 29 leading zeros after the binary point, the double
118 // precision result is accurate to 33 bits.
reduce_large(uint32_t xi,int * np)119 static inline double reduce_large(uint32_t xi, int *np) {
120 const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
121 int shift = (xi >> 23) & 7;
122 uint64_t n, res0, res1, res2;
123
124 xi = (xi & 0xffffff) | 0x800000;
125 xi <<= shift;
126
127 res0 = xi * arr[0];
128 res1 = (uint64_t)xi * arr[4];
129 res2 = (uint64_t)xi * arr[8];
130 res0 = (res2 >> 32) | (res0 << 32);
131 res0 += res1;
132
133 n = (res0 + (1ULL << 61)) >> 62;
134 res0 -= n << 62;
135 double x = (int64_t)res0;
136 *np = n;
137 return x * pi63;
138 }
139
140 } // namespace __llvm_libc
141
142 #endif // LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H
143