mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-22 21:10:07 +00:00
96 lines
3.3 KiB
C
96 lines
3.3 KiB
C
/* Used by sinf, cosf and sincosf functions.
|
|
Copyright (C) 2018-2023 Free Software Foundation, Inc.
|
|
This file is part of the GNU C Library.
|
|
|
|
The GNU C Library is free software; you can redistribute it and/or
|
|
modify it under the terms of the GNU Lesser General Public
|
|
License as published by the Free Software Foundation; either
|
|
version 2.1 of the License, or (at your option) any later version.
|
|
|
|
The GNU C Library is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public
|
|
License along with the GNU C Library; if not, see
|
|
<https://www.gnu.org/licenses/>. */
|
|
|
|
#include <stdint.h>
|
|
#include <math.h>
|
|
#include "math_config.h"
|
|
#include <sincosf_poly.h>
|
|
|
|
/* 2PI * 2^-64. */
|
|
static const double pi63 = 0x1.921FB54442D18p-62;
|
|
/* PI / 4. */
|
|
static const float pio4 = 0x1.921FB6p-1f;
|
|
|
|
/* Polynomial data (the cosine polynomial is negated in the 2nd entry). */
|
|
extern const sincos_t __sincosf_table[2] attribute_hidden;
|
|
|
|
/* Table with 4/PI to 192 bit precision. */
|
|
extern const uint32_t __inv_pio4[] attribute_hidden;
|
|
|
|
/* Top 12 bits of the float representation with the sign bit cleared. */
|
|
static inline uint32_t
|
|
abstop12 (float x)
|
|
{
|
|
return (asuint (x) >> 20) & 0x7ff;
|
|
}
|
|
|
|
/* Fast range reduction using single multiply-subtract. Return the modulo of
|
|
X as a value between -PI/4 and PI/4 and store the quadrant in NP.
|
|
The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double
|
|
is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4,
|
|
the result is accurate for |X| <= 120.0. */
|
|
static inline double
|
|
reduce_fast (double x, const sincos_t *p, int *np)
|
|
{
|
|
double r;
|
|
#if TOINT_INTRINSICS
|
|
/* Use fast round and lround instructions when available. */
|
|
r = x * p->hpi_inv;
|
|
*np = converttoint (r);
|
|
return x - roundtoint (r) * p->hpi;
|
|
#else
|
|
/* Use scaled float to int conversion with explicit rounding.
|
|
hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31.
|
|
This avoids inaccuracies introduced by truncating negative values. */
|
|
r = x * p->hpi_inv;
|
|
int n = ((int32_t)r + 0x800000) >> 24;
|
|
*np = n;
|
|
return x - n * p->hpi;
|
|
#endif
|
|
}
|
|
|
|
/* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic.
|
|
XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored).
|
|
Return the modulo between -PI/4 and PI/4 and store the quadrant in NP.
|
|
Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit
|
|
multiply computes the exact 2.62-bit fixed-point modulo. Since the result
|
|
can have at most 29 leading zeros after the binary point, the double
|
|
precision result is accurate to 33 bits. */
|
|
static inline double
|
|
reduce_large (uint32_t xi, int *np)
|
|
{
|
|
const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15];
|
|
int shift = (xi >> 23) & 7;
|
|
uint64_t n, res0, res1, res2;
|
|
|
|
xi = (xi & 0xffffff) | 0x800000;
|
|
xi <<= shift;
|
|
|
|
res0 = xi * arr[0];
|
|
res1 = (uint64_t)xi * arr[4];
|
|
res2 = (uint64_t)xi * arr[8];
|
|
res0 = (res2 >> 32) | (res0 << 32);
|
|
res0 += res1;
|
|
|
|
n = (res0 + (1ULL << 61)) >> 62;
|
|
res0 -= n << 62;
|
|
double x = (int64_t)res0;
|
|
*np = n;
|
|
return x * pi63;
|
|
}
|