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. 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. 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. 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. 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. 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