|
Packit |
6c4009 |
/* Compute cosine of argument.
|
|
Packit |
6c4009 |
Copyright (C) 2017-2018 Free Software Foundation, Inc.
|
|
Packit |
6c4009 |
This file is part of the GNU C Library.
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
The GNU C Library is free software; you can redistribute it and/or
|
|
Packit |
6c4009 |
modify it under the terms of the GNU Lesser General Public
|
|
Packit |
6c4009 |
License as published by the Free Software Foundation; either
|
|
Packit |
6c4009 |
version 2.1 of the License, or (at your option) any later version.
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
The GNU C Library is distributed in the hope that it will be useful,
|
|
Packit |
6c4009 |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
6c4009 |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
6c4009 |
Lesser General Public License for more details.
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
You should have received a copy of the GNU Lesser General Public
|
|
Packit |
6c4009 |
License along with the GNU C Library; if not, see
|
|
Packit |
6c4009 |
<http://www.gnu.org/licenses/>. */
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#include <errno.h>
|
|
Packit |
6c4009 |
#include <math.h>
|
|
Packit |
6c4009 |
#include <math_private.h>
|
|
Packit |
6c4009 |
#include <libm-alias-float.h>
|
|
Packit |
6c4009 |
#include "s_sincosf.h"
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#ifndef COSF
|
|
Packit |
6c4009 |
# define COSF_FUNC __cosf
|
|
Packit |
6c4009 |
#else
|
|
Packit |
6c4009 |
# define COSF_FUNC COSF
|
|
Packit |
6c4009 |
#endif
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
float
|
|
Packit |
6c4009 |
COSF_FUNC (float x)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
double theta = x;
|
|
Packit |
6c4009 |
double abstheta = fabs (theta);
|
|
Packit |
6c4009 |
if (isless (abstheta, M_PI_4))
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
double cx;
|
|
Packit |
6c4009 |
if (abstheta >= 0x1p-5)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
const double theta2 = theta * theta;
|
|
Packit |
6c4009 |
/* Chebyshev polynomial of the form for cos:
|
|
Packit |
6c4009 |
* 1 + x^2 (C0 + x^2 (C1 + x^2 (C2 + x^2 (C3 + x^2 * C4)))). */
|
|
Packit |
6c4009 |
cx = C3 + theta2 * C4;
|
|
Packit |
6c4009 |
cx = C2 + theta2 * cx;
|
|
Packit |
6c4009 |
cx = C1 + theta2 * cx;
|
|
Packit |
6c4009 |
cx = C0 + theta2 * cx;
|
|
Packit |
6c4009 |
cx = 1. + theta2 * cx;
|
|
Packit |
6c4009 |
return cx;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else if (abstheta >= 0x1p-27)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
/* A simpler Chebyshev approximation is close enough for this range:
|
|
Packit |
6c4009 |
* 1 + x^2 (CC0 + x^3 * CC1). */
|
|
Packit |
6c4009 |
const double theta2 = theta * theta;
|
|
Packit |
6c4009 |
cx = CC0 + theta * theta2 * CC1;
|
|
Packit |
6c4009 |
cx = 1.0 + theta2 * cx;
|
|
Packit |
6c4009 |
return cx;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
/* For small enough |theta|, this is close enough. */
|
|
Packit |
6c4009 |
return 1.0 - abstheta;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else /* |theta| >= Pi/4. */
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
if (isless (abstheta, 9 * M_PI_4))
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
/* There are cases where FE_UPWARD rounding mode can
|
|
Packit |
6c4009 |
produce a result of abstheta * inv_PI_4 == 9,
|
|
Packit |
6c4009 |
where abstheta < 9pi/4, so the domain for
|
|
Packit |
6c4009 |
pio2_table must go to 5 (9 / 2 + 1). */
|
|
Packit |
6c4009 |
unsigned int n = (abstheta * inv_PI_4) + 1;
|
|
Packit |
6c4009 |
theta = abstheta - pio2_table[n / 2];
|
|
Packit |
6c4009 |
return reduced_cos (theta, n);
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else if (isless (abstheta, INFINITY))
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
if (abstheta < 0x1p+23)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
|
|
Packit |
6c4009 |
double x = n / 2;
|
|
Packit |
6c4009 |
theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
|
|
Packit |
6c4009 |
/* Argument reduction needed. */
|
|
Packit |
6c4009 |
return reduced_cos (theta, n);
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else /* |theta| >= 2^23. */
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
x = fabsf (x);
|
|
Packit |
6c4009 |
int exponent;
|
|
Packit |
6c4009 |
GET_FLOAT_WORD (exponent, x);
|
|
Packit |
6c4009 |
exponent = (exponent >> FLOAT_EXPONENT_SHIFT)
|
|
Packit |
6c4009 |
- FLOAT_EXPONENT_BIAS;
|
|
Packit |
6c4009 |
exponent += 3;
|
|
Packit |
6c4009 |
exponent /= 28;
|
|
Packit |
6c4009 |
double a = invpio4_table[exponent] * x;
|
|
Packit |
6c4009 |
double b = invpio4_table[exponent + 1] * x;
|
|
Packit |
6c4009 |
double c = invpio4_table[exponent + 2] * x;
|
|
Packit |
6c4009 |
double d = invpio4_table[exponent + 3] * x;
|
|
Packit |
6c4009 |
uint64_t l = a;
|
|
Packit |
6c4009 |
l &= ~0x7;
|
|
Packit |
6c4009 |
a -= l;
|
|
Packit |
6c4009 |
double e = a + b;
|
|
Packit |
6c4009 |
l = e;
|
|
Packit |
6c4009 |
e = a - l;
|
|
Packit |
6c4009 |
if (l & 1)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
e -= 1.0;
|
|
Packit |
6c4009 |
e += b;
|
|
Packit |
6c4009 |
e += c;
|
|
Packit |
6c4009 |
e += d;
|
|
Packit |
6c4009 |
e *= M_PI_4;
|
|
Packit |
6c4009 |
return reduced_cos (e, l + 1);
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
e += b;
|
|
Packit |
6c4009 |
e += c;
|
|
Packit |
6c4009 |
e += d;
|
|
Packit |
6c4009 |
if (e <= 1.0)
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
e *= M_PI_4;
|
|
Packit |
6c4009 |
return reduced_cos (e, l + 1);
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
l++;
|
|
Packit |
6c4009 |
e -= 2.0;
|
|
Packit |
6c4009 |
e *= M_PI_4;
|
|
Packit |
6c4009 |
return reduced_cos (e, l + 1);
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
else
|
|
Packit |
6c4009 |
{
|
|
Packit |
6c4009 |
int32_t ix;
|
|
Packit |
6c4009 |
GET_FLOAT_WORD (ix, abstheta);
|
|
Packit |
6c4009 |
/* cos(Inf or NaN) is NaN. */
|
|
Packit |
6c4009 |
if (ix == 0x7f800000) /* Inf. */
|
|
Packit |
6c4009 |
__set_errno (EDOM);
|
|
Packit |
6c4009 |
return x - x;
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
}
|
|
Packit |
6c4009 |
|
|
Packit |
6c4009 |
#ifndef COSF
|
|
Packit |
6c4009 |
libm_alias_float (__cos, cos)
|
|
Packit |
6c4009 |
#endif
|