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