|
Packit Service |
82fcde |
/*
|
|
Packit Service |
82fcde |
* ====================================================
|
|
Packit Service |
82fcde |
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
Packit Service |
82fcde |
*
|
|
Packit Service |
82fcde |
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
Packit Service |
82fcde |
* Permission to use, copy, modify, and distribute this
|
|
Packit Service |
82fcde |
* software is freely granted, provided that this notice
|
|
Packit Service |
82fcde |
* is preserved.
|
|
Packit Service |
82fcde |
* ====================================================
|
|
Packit Service |
82fcde |
*/
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
/*
|
|
Packit Service |
82fcde |
Long double expansions are
|
|
Packit Service |
82fcde |
Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
|
|
Packit Service |
82fcde |
and are incorporated herein by permission of the author. The author
|
|
Packit Service |
82fcde |
reserves the right to distribute this material elsewhere under different
|
|
Packit Service |
82fcde |
copying permissions. These modifications are distributed here under
|
|
Packit Service |
82fcde |
the following terms:
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
This library is free software; you can redistribute it and/or
|
|
Packit Service |
82fcde |
modify it under the terms of the GNU Lesser General Public
|
|
Packit Service |
82fcde |
License as published by the Free Software Foundation; either
|
|
Packit Service |
82fcde |
version 2.1 of the License, or (at your option) any later version.
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
This library is distributed in the hope that it will be useful,
|
|
Packit Service |
82fcde |
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit Service |
82fcde |
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit Service |
82fcde |
Lesser General Public License for more details.
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
You should have received a copy of the GNU Lesser General Public
|
|
Packit Service |
82fcde |
License along with this library; if not, see
|
|
Packit Service |
82fcde |
<http://www.gnu.org/licenses/>. */
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
/* __kernel_tanl( x, y, k )
|
|
Packit Service |
82fcde |
* kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
|
|
Packit Service |
82fcde |
* Input x is assumed to be bounded by ~pi/4 in magnitude.
|
|
Packit Service |
82fcde |
* Input y is the tail of x.
|
|
Packit Service |
82fcde |
* Input k indicates whether tan (if k=1) or
|
|
Packit Service |
82fcde |
* -1/tan (if k= -1) is returned.
|
|
Packit Service |
82fcde |
*
|
|
Packit Service |
82fcde |
* Algorithm
|
|
Packit Service |
82fcde |
* 1. Since tan(-x) = -tan(x), we need only to consider positive x.
|
|
Packit Service |
82fcde |
* 2. if x < 2^-33, return x with inexact if x!=0.
|
|
Packit Service |
82fcde |
* 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
|
|
Packit Service |
82fcde |
* on [0,0.67433].
|
|
Packit Service |
82fcde |
*
|
|
Packit Service |
82fcde |
* Note: tan(x+y) = tan(x) + tan'(x)*y
|
|
Packit Service |
82fcde |
* ~ tan(x) + (1+x*x)*y
|
|
Packit Service |
82fcde |
* Therefore, for better accuracy in computing tan(x+y), let
|
|
Packit Service |
82fcde |
* r = x^3 * R(x^2)
|
|
Packit Service |
82fcde |
* then
|
|
Packit Service |
82fcde |
* tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
|
|
Packit Service |
82fcde |
*
|
|
Packit Service |
82fcde |
* 4. For x in [0.67433,pi/4], let y = pi/4 - x, then
|
|
Packit Service |
82fcde |
* tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
|
|
Packit Service |
82fcde |
* = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
|
|
Packit Service |
82fcde |
*/
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
#include <float.h>
|
|
Packit Service |
82fcde |
#include <math.h>
|
|
Packit Service |
82fcde |
#include <math_private.h>
|
|
Packit Service |
82fcde |
#include <math-underflow.h>
|
|
Packit Service |
82fcde |
#include <libc-diag.h>
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
static const long double
|
|
Packit Service |
82fcde |
one = 1.0L,
|
|
Packit Service |
82fcde |
pio4hi = 0xc.90fdaa22168c235p-4L,
|
|
Packit Service |
82fcde |
pio4lo = -0x3.b399d747f23e32ecp-68L,
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
/* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
|
|
Packit Service |
82fcde |
0 <= x <= 0.6743316650390625
|
|
Packit Service |
82fcde |
Peak relative error 8.0e-36 */
|
|
Packit Service |
82fcde |
TH = 3.333333333333333333333333333333333333333E-1L,
|
|
Packit Service |
82fcde |
T0 = -1.813014711743583437742363284336855889393E7L,
|
|
Packit Service |
82fcde |
T1 = 1.320767960008972224312740075083259247618E6L,
|
|
Packit Service |
82fcde |
T2 = -2.626775478255838182468651821863299023956E4L,
|
|
Packit Service |
82fcde |
T3 = 1.764573356488504935415411383687150199315E2L,
|
|
Packit Service |
82fcde |
T4 = -3.333267763822178690794678978979803526092E-1L,
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
U0 = -1.359761033807687578306772463253710042010E8L,
|
|
Packit Service |
82fcde |
U1 = 6.494370630656893175666729313065113194784E7L,
|
|
Packit Service |
82fcde |
U2 = -4.180787672237927475505536849168729386782E6L,
|
|
Packit Service |
82fcde |
U3 = 8.031643765106170040139966622980914621521E4L,
|
|
Packit Service |
82fcde |
U4 = -5.323131271912475695157127875560667378597E2L;
|
|
Packit Service |
82fcde |
/* 1.000000000000000000000000000000000000000E0 */
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
long double
|
|
Packit Service |
82fcde |
__kernel_tanl (long double x, long double y, int iy)
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
long double z, r, v, w, s;
|
|
Packit Service |
82fcde |
long double absx = fabsl (x);
|
|
Packit Service |
82fcde |
int sign;
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
if (absx < 0x1p-33)
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
if ((int) x == 0)
|
|
Packit Service |
82fcde |
{ /* generate inexact */
|
|
Packit Service |
82fcde |
if (x == 0 && iy == -1)
|
|
Packit Service |
82fcde |
return one / fabsl (x);
|
|
Packit Service |
82fcde |
else if (iy == 1)
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
math_check_force_underflow_nonneg (absx);
|
|
Packit Service |
82fcde |
return x;
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
else
|
|
Packit Service |
82fcde |
return -one / x;
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
if (absx >= 0.6743316650390625L)
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
if (signbit (x))
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
x = -x;
|
|
Packit Service |
82fcde |
y = -y;
|
|
Packit Service |
82fcde |
sign = -1;
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
else
|
|
Packit Service |
82fcde |
sign = 1;
|
|
Packit Service |
82fcde |
z = pio4hi - x;
|
|
Packit Service |
82fcde |
w = pio4lo - y;
|
|
Packit Service |
82fcde |
x = z + w;
|
|
Packit Service |
82fcde |
y = 0.0;
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
z = x * x;
|
|
Packit Service |
82fcde |
r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
|
|
Packit Service |
82fcde |
v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
|
|
Packit Service |
82fcde |
r = r / v;
|
|
Packit Service |
82fcde |
|
|
Packit Service |
82fcde |
s = z * x;
|
|
Packit Service |
82fcde |
r = y + z * (s * r + y);
|
|
Packit Service |
82fcde |
r += TH * s;
|
|
Packit Service |
82fcde |
w = x + r;
|
|
Packit Service |
82fcde |
if (absx >= 0.6743316650390625L)
|
|
Packit Service |
82fcde |
{
|
|
Packit Service |
82fcde |
v = (long double) iy;
|
|
Packit Service |
82fcde |
w = (v - 2.0 * (x - (w * w / (w + v) - r)));
|
|
Packit Service |
82fcde |
/* SIGN is set for arguments that reach this code, but not
|
|
Packit Service |
82fcde |
otherwise, resulting in warnings that it may be used
|
|
Packit Service |
82fcde |
uninitialized although in the cases where it is used it has
|
|
Packit Service |
82fcde |
always been set. */
|
|
Packit Service |
82fcde |
DIAG_PUSH_NEEDS_COMMENT;
|
|
Packit Service |
82fcde |
DIAG_IGNORE_NEEDS_COMMENT (4.8, "-Wmaybe-uninitialized");
|
|
Packit Service |
82fcde |
if (sign < 0)
|
|
Packit Service |
82fcde |
w = -w;
|
|
Packit Service |
82fcde |
DIAG_POP_NEEDS_COMMENT;
|
|
Packit Service |
82fcde |
return w;
|
|
Packit Service |
82fcde |
}
|
|
Packit Service |
82fcde |
if (iy == 1)
|
|
Packit Service |
82fcde |
return w;
|
|
Packit Service |
82fcde |
else
|
|
Packit Service |
82fcde |
return -1.0 / (x + r);
|
|
Packit Service |
82fcde |
}
|