4f889260c3
Mark Murray graciously provided access to an aarch64 system to test the ld128 implementations. This patch address * Misuses of copysignl() in sinpil() and tanpil(). * Redo the splitting of argument 'x' into an integer part and remainder. The remainder must satify 0 <= r < 1. * Update the reduction of the integer part to something that can easily be seen as even or odd, e.g., sin(pi*x) = (-1)^n*sin(pi*r) with n <= 2^112 and we an reduce n by subtracting integer powers of 2. * In s_cospil.c, fix typos where 'x' is used where 'ax', the remainder, is required. * In tanpil(), fix the use of an uninitialized variable, ax = fabsl(ax), ax should be x in fabsl(). One item of note, in the limited tested on aarch64, the max ULP for sinpil() and cospil() were less than 1.1 ULP, which is higher that the desired max ULP less than 1. This was traced to the kernel for cosl() in the fundamental interval [0,pi/4]. The coefficients in the minmax polynomial likely need refinement. PR: 218514 MFC after: 1 week
124 lines
3.2 KiB
C
124 lines
3.2 KiB
C
/*-
|
|
* Copyright (c) 2017-2021 Steven G. Kargl
|
|
* All rights reserved.
|
|
*
|
|
* Redistribution and use in source and binary forms, with or without
|
|
* modification, are permitted provided that the following conditions
|
|
* are met:
|
|
* 1. Redistributions of source code must retain the above copyright
|
|
* notice unmodified, this list of conditions, and the following
|
|
* disclaimer.
|
|
* 2. Redistributions in binary form must reproduce the above copyright
|
|
* notice, this list of conditions and the following disclaimer in the
|
|
* documentation and/or other materials provided with the distribution.
|
|
*
|
|
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
|
|
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
|
|
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
|
|
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
|
|
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
|
|
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
|
|
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
|
|
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
|
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
|
|
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
*/
|
|
|
|
/*
|
|
* See ../src/s_tanpi.c for implementation details.
|
|
*/
|
|
|
|
#include "math.h"
|
|
#include "math_private.h"
|
|
|
|
/*
|
|
* pi_hi contains the leading 56 bits of a 169 bit approximation for pi.
|
|
*/
|
|
static const long double
|
|
pi_hi = 3.14159265358979322702026593105983920e+00L,
|
|
pi_lo = 1.14423774522196636802434264184180742e-17L;
|
|
|
|
static inline long double
|
|
__kernel_tanpil(long double x)
|
|
{
|
|
long double hi, lo, t;
|
|
|
|
if (x < 0.25) {
|
|
hi = (double)x;
|
|
lo = x - hi;
|
|
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
|
|
hi *= pi_hi;
|
|
_2sumF(hi, lo);
|
|
t = __kernel_tanl(hi, lo, -1);
|
|
} else if (x > 0.25) {
|
|
x = 0.5 - x;
|
|
hi = (double)x;
|
|
lo = x - hi;
|
|
lo = lo * (pi_lo + pi_hi) + hi * pi_lo;
|
|
hi *= pi_hi;
|
|
_2sumF(hi, lo);
|
|
t = - __kernel_tanl(hi, lo, 1);
|
|
} else
|
|
t = 1;
|
|
|
|
return (t);
|
|
}
|
|
|
|
volatile static const double vzero = 0;
|
|
|
|
long double
|
|
tanpil(long double x)
|
|
{
|
|
long double ax, hi, lo, xf, t;
|
|
uint32_t ix;
|
|
|
|
ax = fabsl(x);
|
|
|
|
if (ax < 1) {
|
|
if (ax < 0.5) {
|
|
if (ax < 0x1p-60) {
|
|
if (x == 0)
|
|
return (x);
|
|
hi = (double)x;
|
|
hi *= 0x1p113L;
|
|
lo = x * 0x1p113L - hi;
|
|
t = (pi_lo + pi_hi) * lo + pi_lo * lo +
|
|
pi_hi * hi;
|
|
return (t * 0x1p-113L);
|
|
}
|
|
t = __kernel_tanpil(ax);
|
|
} else if (ax == 0.5)
|
|
return ((ax - ax) / (ax - ax));
|
|
else
|
|
t = -__kernel_tanpil(1 - ax);
|
|
return (x < 0 ? -t : t);
|
|
}
|
|
|
|
if (ix < 0x1p112) {
|
|
/* Split x = n + r with 0 <= r < 1. */
|
|
xf = (ax + 0x1p112L) - 0x1p112L; /* Integer part */
|
|
ax -= xf; /* Remainder */
|
|
if (ax < 0) {
|
|
ax += 1;
|
|
xf -= 1;
|
|
}
|
|
|
|
if (ax < 0.5)
|
|
t = ax == 0 ? 0 : __kernel_tanpil(ax);
|
|
else if (ax == 0.5)
|
|
return ((ax - ax) / (ax - ax));
|
|
else
|
|
t = -__kernel_tanpil(1 - ax);
|
|
return (x < 0 ? -t : t);
|
|
}
|
|
|
|
/* x = +-inf or nan. */
|
|
if (isinf(x) || isnan(x))
|
|
return (vzero / vzero);
|
|
|
|
/*
|
|
* |x| >= 0x1p112 is always an integer, so return +-0.
|
|
*/
|
|
return (copysignl(0, x));
|
|
}
|