1/* Compute cosine of argument.
2 Copyright (C) 2017-2018 Free Software Foundation, Inc.
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
17 <http://www.gnu.org/licenses/>. */
18
19#include <errno.h>
20#include <math.h>
21#include <math_private.h>
22#include <libm-alias-float.h>
23#include "s_sincosf.h"
24
25#ifndef COSF
26# define COSF_FUNC __cosf
27#else
28# define COSF_FUNC COSF
29#endif
30
31float
32COSF_FUNC (float x)
33{
34 double theta = x;
35 double abstheta = fabs (theta);
36 if (isless (abstheta, M_PI_4))
37 {
38 double cx;
39 if (abstheta >= 0x1p-5)
40 {
41 const double theta2 = theta * theta;
42 /* Chebyshev polynomial of the form for cos:
43 * 1 + x^2 (C0 + x^2 (C1 + x^2 (C2 + x^2 (C3 + x^2 * C4)))). */
44 cx = C3 + theta2 * C4;
45 cx = C2 + theta2 * cx;
46 cx = C1 + theta2 * cx;
47 cx = C0 + theta2 * cx;
48 cx = 1. + theta2 * cx;
49 return cx;
50 }
51 else if (abstheta >= 0x1p-27)
52 {
53 /* A simpler Chebyshev approximation is close enough for this range:
54 * 1 + x^2 (CC0 + x^3 * CC1). */
55 const double theta2 = theta * theta;
56 cx = CC0 + theta * theta2 * CC1;
57 cx = 1.0 + theta2 * cx;
58 return cx;
59 }
60 else
61 {
62 /* For small enough |theta|, this is close enough. */
63 return 1.0 - abstheta;
64 }
65 }
66 else /* |theta| >= Pi/4. */
67 {
68 if (isless (abstheta, 9 * M_PI_4))
69 {
70 /* There are cases where FE_UPWARD rounding mode can
71 produce a result of abstheta * inv_PI_4 == 9,
72 where abstheta < 9pi/4, so the domain for
73 pio2_table must go to 5 (9 / 2 + 1). */
74 unsigned int n = (abstheta * inv_PI_4) + 1;
75 theta = abstheta - pio2_table[n / 2];
76 return reduced_cos (theta, n);
77 }
78 else if (isless (abstheta, INFINITY))
79 {
80 if (abstheta < 0x1p+23)
81 {
82 unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
83 double x = n / 2;
84 theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
85 /* Argument reduction needed. */
86 return reduced_cos (theta, n);
87 }
88 else /* |theta| >= 2^23. */
89 {
90 x = fabsf (x);
91 int exponent;
92 GET_FLOAT_WORD (exponent, x);
93 exponent = (exponent >> FLOAT_EXPONENT_SHIFT)
94 - FLOAT_EXPONENT_BIAS;
95 exponent += 3;
96 exponent /= 28;
97 double a = invpio4_table[exponent] * x;
98 double b = invpio4_table[exponent + 1] * x;
99 double c = invpio4_table[exponent + 2] * x;
100 double d = invpio4_table[exponent + 3] * x;
101 uint64_t l = a;
102 l &= ~0x7;
103 a -= l;
104 double e = a + b;
105 l = e;
106 e = a - l;
107 if (l & 1)
108 {
109 e -= 1.0;
110 e += b;
111 e += c;
112 e += d;
113 e *= M_PI_4;
114 return reduced_cos (e, l + 1);
115 }
116 else
117 {
118 e += b;
119 e += c;
120 e += d;
121 if (e <= 1.0)
122 {
123 e *= M_PI_4;
124 return reduced_cos (e, l + 1);
125 }
126 else
127 {
128 l++;
129 e -= 2.0;
130 e *= M_PI_4;
131 return reduced_cos (e, l + 1);
132 }
133 }
134 }
135 }
136 else
137 {
138 int32_t ix;
139 GET_FLOAT_WORD (ix, abstheta);
140 /* cos(Inf or NaN) is NaN. */
141 if (ix == 0x7f800000) /* Inf. */
142 __set_errno (EDOM);
143 return x - x;
144 }
145 }
146}
147
148#ifndef COSF
149libm_alias_float (__cos, cos)
150#endif
151