]>
Commit | Line | Data |
---|---|---|
984ae996 | 1 | /* Used by sinf, cosf and sincosf functions. |
6d7e8eda | 2 | Copyright (C) 2018-2023 Free Software Foundation, Inc. |
984ae996 RS |
3 | This file is part of the GNU C Library. |
4 | ||
5 | The GNU C Library is free software; you can redistribute it and/or | |
6 | modify it under the terms of the GNU Lesser General Public | |
7 | License as published by the Free Software Foundation; either | |
8 | version 2.1 of the License, or (at your option) any later version. | |
9 | ||
10 | The GNU C Library is distributed in the hope that it will be useful, | |
11 | but WITHOUT ANY WARRANTY; without even the implied warranty of | |
12 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
13 | Lesser General Public License for more details. | |
14 | ||
15 | You should have received a copy of the GNU Lesser General Public | |
16 | License along with the GNU C Library; if not, see | |
5a82c748 | 17 | <https://www.gnu.org/licenses/>. */ |
984ae996 | 18 | |
ea5c662c WD |
19 | #include <stdint.h> |
20 | #include <math.h> | |
21 | #include "math_config.h" | |
8700a785 | 22 | #include <sincosf_poly.h> |
ea5c662c | 23 | |
ea5c662c WD |
24 | /* 2PI * 2^-64. */ |
25 | static const double pi63 = 0x1.921FB54442D18p-62; | |
26 | /* PI / 4. */ | |
d3e4f5a1 | 27 | static const float pio4 = 0x1.921FB6p-1f; |
ea5c662c | 28 | |
ea5c662c WD |
29 | /* Polynomial data (the cosine polynomial is negated in the 2nd entry). */ |
30 | extern const sincos_t __sincosf_table[2] attribute_hidden; | |
31 | ||
32 | /* Table with 4/PI to 192 bit precision. */ | |
33 | extern const uint32_t __inv_pio4[] attribute_hidden; | |
34 | ||
35 | /* Top 12 bits of the float representation with the sign bit cleared. */ | |
36 | static inline uint32_t | |
37 | abstop12 (float x) | |
38 | { | |
39 | return (asuint (x) >> 20) & 0x7ff; | |
40 | } | |
41 | ||
ea5c662c WD |
42 | /* Fast range reduction using single multiply-subtract. Return the modulo of |
43 | X as a value between -PI/4 and PI/4 and store the quadrant in NP. | |
44 | The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double | |
45 | is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4, | |
46 | the result is accurate for |X| <= 120.0. */ | |
47 | static inline double | |
48 | reduce_fast (double x, const sincos_t *p, int *np) | |
49 | { | |
50 | double r; | |
51 | #if TOINT_INTRINSICS | |
52 | /* Use fast round and lround instructions when available. */ | |
53 | r = x * p->hpi_inv; | |
54 | *np = converttoint (r); | |
55 | return x - roundtoint (r) * p->hpi; | |
56 | #else | |
57 | /* Use scaled float to int conversion with explicit rounding. | |
58 | hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31. | |
59 | This avoids inaccuracies introduced by truncating negative values. */ | |
60 | r = x * p->hpi_inv; | |
61 | int n = ((int32_t)r + 0x800000) >> 24; | |
62 | *np = n; | |
63 | return x - n * p->hpi; | |
64 | #endif | |
65 | } | |
66 | ||
67 | /* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic. | |
68 | XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored). | |
69 | Return the modulo between -PI/4 and PI/4 and store the quadrant in NP. | |
70 | Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit | |
71 | multiply computes the exact 2.62-bit fixed-point modulo. Since the result | |
72 | can have at most 29 leading zeros after the binary point, the double | |
73 | precision result is accurate to 33 bits. */ | |
74 | static inline double | |
75 | reduce_large (uint32_t xi, int *np) | |
76 | { | |
77 | const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15]; | |
78 | int shift = (xi >> 23) & 7; | |
79 | uint64_t n, res0, res1, res2; | |
80 | ||
81 | xi = (xi & 0xffffff) | 0x800000; | |
82 | xi <<= shift; | |
83 | ||
84 | res0 = xi * arr[0]; | |
85 | res1 = (uint64_t)xi * arr[4]; | |
86 | res2 = (uint64_t)xi * arr[8]; | |
87 | res0 = (res2 >> 32) | (res0 << 32); | |
88 | res0 += res1; | |
89 | ||
90 | n = (res0 + (1ULL << 61)) >> 62; | |
91 | res0 -= n << 62; | |
92 | double x = (int64_t)res0; | |
93 | *np = n; | |
94 | return x * pi63; | |
95 | } |