AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-12-25 20:21:07 +00:00
Code Issues Packages Projects Releases Wiki Activity
9be7530cc0
glibc/sysdeps/ieee754/dbl-64/s_tan.c
Joseph Myers 38722448c6 Use libm_alias_double for dbl-64 atan, tan. 
This patch makes the dbl-64 atan and tan implementations use
libm_alias_double, removing the corresponding ldbl-opt wrappers.

Tested for x86_64, and with build-many-glibcs.py.  Installed stripped
shared libraries are unchanged on non-ldbl-opt platforms.  For
ldbl-opt configurations, the patch has the effect of causing
compat_symbol to define atanl and tanl in terms of __atan and __tan
instead of in terms of atan and tan, which is enough to change the
installed stripped libm.so.

	* sysdeps/ieee754/dbl-64/s_atan.c: Include <libm-alias-double.h>.
	(atan): Define using libm_alias_double.
	* sysdeps/ieee754/dbl-64/s_tan.c: Include <libm-alias-double.h>.
	(tan): Define using libm_alias_double.
	* sysdeps/ieee754/ldbl-opt/s_atan.c: Remove file.
	* sysdeps/ieee754/ldbl-opt/s_tan.c: Likewise.
2017-10-02 23:16:56 +00:00

850 lines
22 KiB
C
Raw Blame History

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001-2017 Free Software Foundation, Inc.
*
* This program 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.
*
* This program 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 this program; if not, see <http://www.gnu.org/licenses/>.
*/
/*********************************************************************/
/* MODULE_NAME: utan.c */
/* */
/* FUNCTIONS: utan */
/* tanMp */
/* */
/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h utan.h */
/* branred.c sincos32.c mptan.c */
/* utan.tbl */
/* */
/* An ultimate tan routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of tan(x). */
/* Assumption: Machine arithmetic operations are performed in */
/* round to nearest mode of IEEE 754 standard. */
/* */
/*********************************************************************/
#include <errno.h>
#include <float.h>
#include "endian.h"
#include <dla.h>
#include "mpa.h"
#include "MathLib.h"
#include <math.h>
#include <math_private.h>
#include <libm-alias-double.h>
#include <fenv.h>
#include <stap-probe.h>
#ifndef SECTION
# define SECTION
#endif
static double tanMp (double);
void __mptan (double, mp_no *, int);
double
SECTION
__tan (double x)
{
#include "utan.h"
#include "utan.tbl"
int ux, i, n;
double a, da, a2, b, db, c, dc, c1, cc1, c2, cc2, c3, cc3, fi, ffi, gi, pz,
s, sy, t, t1, t2, t3, t4, t7, t8, t9, t10, w, x2, xn, xx2, y, ya,
yya, z0, z, zz, z2, zz2;
#ifndef DLA_FMS
double t5, t6;
#endif
int p;
number num, v;
mp_no mpa, mpt1, mpt2;
double retval;
int __branred (double, double *, double *);
int __mpranred (double, mp_no *, int);
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
/* x=+-INF, x=NaN */
num.d = x;
ux = num.i[HIGH_HALF];
if ((ux & 0x7ff00000) == 0x7ff00000)
{
if ((ux & 0x7fffffff) == 0x7ff00000)
__set_errno (EDOM);
retval = x - x;
goto ret;
}
w = (x < 0.0) ? -x : x;
/* (I) The case abs(x) <= 1.259e-8 */
if (w <= g1.d)
{
math_check_force_underflow_nonneg (w);
retval = x;
goto ret;
}
/* (II) The case 1.259e-8 < abs(x) <= 0.0608 */
if (w <= g2.d)
{
/* First stage */
x2 = x * x;
t2 = d9.d + x2 * d11.d;
t2 = d7.d + x2 * t2;
t2 = d5.d + x2 * t2;
t2 = d3.d + x2 * t2;
t2 *= x * x2;
if ((y = x + (t2 - u1.d * t2)) == x + (t2 + u1.d * t2))
{
retval = y;
goto ret;
}
/* Second stage */
c1 = a25.d + x2 * a27.d;
c1 = a23.d + x2 * c1;
c1 = a21.d + x2 * c1;
c1 = a19.d + x2 * c1;
c1 = a17.d + x2 * c1;
c1 = a15.d + x2 * c1;
c1 *= x2;
EMULV (x, x, x2, xx2, t1, t2, t3, t4, t5);
ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (x, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (x, 0.0, c2, cc2, c1, cc1, t1, t2);
if ((y = c1 + (cc1 - u2.d * c1)) == c1 + (cc1 + u2.d * c1))
{
retval = y;
goto ret;
}
retval = tanMp (x);
goto ret;
}
/* (III) The case 0.0608 < abs(x) <= 0.787 */
if (w <= g3.d)
{
/* First stage */
i = ((int) (mfftnhf.d + TWO8 * w));
z = w - xfg[i][0].d;
z2 = z * z;
s = (x < 0.0) ? -1 : 1;
pz = z + z * z2 * (e0.d + z2 * e1.d);
fi = xfg[i][1].d;
gi = xfg[i][2].d;
t2 = pz * (gi + fi) / (gi - pz);
if ((y = fi + (t2 - fi * u3.d)) == fi + (t2 + fi * u3.d))
{
retval = (s * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = fi * ua3.d + t3 * ub3.d;
if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
{
retval = (s * y);
goto ret;
}
/* Second stage */
ffi = xfg[i][3].d;
c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
EMULV (z, z, z2, zz2, t1, t2, t3, t4, t5);
ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (z, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (z, 0.0, c2, cc2, c1, cc1, t1, t2);
ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8);
SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u4.d * c3)) == c3 + (cc3 + u4.d * c3))
{
retval = (s * y);
goto ret;
}
retval = tanMp (x);
goto ret;
}
/* (---) The case 0.787 < abs(x) <= 25 */
if (w <= g4.d)
{
/* Range reduction by algorithm i */
t = (x * hpinv.d + toint.d);
xn = t - toint.d;
v.d = t;
t1 = (x - xn * mp1.d) - xn * mp2.d;
n = v.i[LOW_HALF] & 0x00000001;
da = xn * mp3.d;
a = t1 - da;
da = (t1 - a) - da;
if (a < 0.0)
{
ya = -a;
yya = -da;
sy = -1;
}
else
{
ya = a;
yya = da;
sy = 1;
}
/* (IV),(V) The case 0.787 < abs(x) <= 25, abs(y) <= 1e-7 */
if (ya <= gy1.d)
{
retval = tanMp (x);
goto ret;
}
/* (VI) The case 0.787 < abs(x) <= 25, 1e-7 < abs(y) <= 0.0608 */
if (ya <= gy2.d)
{
a2 = a * a;
t2 = d9.d + a2 * d11.d;
t2 = d7.d + a2 * t2;
t2 = d5.d + a2 * t2;
t2 = d3.d + a2 * t2;
t2 = da + a * a2 * t2;
if (n)
{
/* First stage -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4, t5, t6, t7, t8,
t9, t10);
if ((y = c + (dc - u6.d * c)) == c + (dc + u6.d * c))
{
retval = (-y);
goto ret;
}
}
else
{
/* First stage tan */
if ((y = a + (t2 - u5.d * a)) == a + (t2 + u5.d * a))
{
retval = y;
goto ret;
}
}
/* Second stage */
/* Range reduction by algorithm ii */
t = (x * hpinv.d + toint.d);
xn = t - toint.d;
v.d = t;
t1 = (x - xn * mp1.d) - xn * mp2.d;
n = v.i[LOW_HALF] & 0x00000001;
da = xn * pp3.d;
t = t1 - da;
da = (t1 - t) - da;
t1 = xn * pp4.d;
a = t - t1;
da = ((t - a) - t1) + da;
/* Second stage */
EADD (a, da, t1, t2);
a = t1;
da = t2;
MUL2 (a, da, a, da, x2, xx2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = a25.d + x2 * a27.d;
c1 = a23.d + x2 * c1;
c1 = a21.d + x2 * c1;
c1 = a19.d + x2 * c1;
c1 = a17.d + x2 * c1;
c1 = a15.d + x2 * c1;
c1 *= x2;
ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (a, da, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
if (n)
{
/* Second stage -cot */
DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7,
t8, t9, t10);
if ((y = c2 + (cc2 - u8.d * c2)) == c2 + (cc2 + u8.d * c2))
{
retval = (-y);
goto ret;
}
}
else
{
/* Second stage tan */
if ((y = c1 + (cc1 - u7.d * c1)) == c1 + (cc1 + u7.d * c1))
{
retval = y;
goto ret;
}
}
retval = tanMp (x);
goto ret;
}
/* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
/* First stage */
i = ((int) (mfftnhf.d + TWO8 * ya));
z = (z0 = (ya - xfg[i][0].d)) + yya;
z2 = z * z;
pz = z + z * z2 * (e0.d + z2 * e1.d);
fi = xfg[i][1].d;
gi = xfg[i][2].d;
if (n)
{
/* -cot */
t2 = pz * (fi + gi) / (fi + pz);
if ((y = gi - (t2 - gi * u10.d)) == gi - (t2 + gi * u10.d))
{
retval = (-sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = gi * ua10.d + t3 * ub10.d;
if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
if ((y = fi + (t2 - fi * u9.d)) == fi + (t2 + fi * u9.d))
{
retval = (sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = fi * ua9.d + t3 * ub9.d;
if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
{
retval = (sy * y);
goto ret;
}
}
/* Second stage */
ffi = xfg[i][3].d;
EADD (z0, yya, z, zz)
MUL2 (z, zz, z, zz, z2, zz2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8);
SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
if (n)
{
/* -cot */
DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u12.d * c3)) == c3 + (cc3 + u12.d * c3))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u11.d * c3)) == c3 + (cc3 + u11.d * c3))
{
retval = (sy * y);
goto ret;
}
}
retval = tanMp (x);
goto ret;
}
/* (---) The case 25 < abs(x) <= 1e8 */
if (w <= g5.d)
{
/* Range reduction by algorithm ii */
t = (x * hpinv.d + toint.d);
xn = t - toint.d;
v.d = t;
t1 = (x - xn * mp1.d) - xn * mp2.d;
n = v.i[LOW_HALF] & 0x00000001;
da = xn * pp3.d;
t = t1 - da;
da = (t1 - t) - da;
t1 = xn * pp4.d;
a = t - t1;
da = ((t - a) - t1) + da;
EADD (a, da, t1, t2);
a = t1;
da = t2;
if (a < 0.0)
{
ya = -a;
yya = -da;
sy = -1;
}
else
{
ya = a;
yya = da;
sy = 1;
}
/* (+++) The case 25 < abs(x) <= 1e8, abs(y) <= 1e-7 */
if (ya <= gy1.d)
{
retval = tanMp (x);
goto ret;
}
/* (VIII) The case 25 < abs(x) <= 1e8, 1e-7 < abs(y) <= 0.0608 */
if (ya <= gy2.d)
{
a2 = a * a;
t2 = d9.d + a2 * d11.d;
t2 = d7.d + a2 * t2;
t2 = d5.d + a2 * t2;
t2 = d3.d + a2 * t2;
t2 = da + a * a2 * t2;
if (n)
{
/* First stage -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4, t5, t6, t7, t8,
t9, t10);
if ((y = c + (dc - u14.d * c)) == c + (dc + u14.d * c))
{
retval = (-y);
goto ret;
}
}
else
{
/* First stage tan */
if ((y = a + (t2 - u13.d * a)) == a + (t2 + u13.d * a))
{
retval = y;
goto ret;
}
}
/* Second stage */
MUL2 (a, da, a, da, x2, xx2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = a25.d + x2 * a27.d;
c1 = a23.d + x2 * c1;
c1 = a21.d + x2 * c1;
c1 = a19.d + x2 * c1;
c1 = a17.d + x2 * c1;
c1 = a15.d + x2 * c1;
c1 *= x2;
ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (a, da, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
if (n)
{
/* Second stage -cot */
DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7,
t8, t9, t10);
if ((y = c2 + (cc2 - u16.d * c2)) == c2 + (cc2 + u16.d * c2))
{
retval = (-y);
goto ret;
}
}
else
{
/* Second stage tan */
if ((y = c1 + (cc1 - u15.d * c1)) == c1 + (cc1 + u15.d * c1))
{
retval = (y);
goto ret;
}
}
retval = tanMp (x);
goto ret;
}
/* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
/* First stage */
i = ((int) (mfftnhf.d + TWO8 * ya));
z = (z0 = (ya - xfg[i][0].d)) + yya;
z2 = z * z;
pz = z + z * z2 * (e0.d + z2 * e1.d);
fi = xfg[i][1].d;
gi = xfg[i][2].d;
if (n)
{
/* -cot */
t2 = pz * (fi + gi) / (fi + pz);
if ((y = gi - (t2 - gi * u18.d)) == gi - (t2 + gi * u18.d))
{
retval = (-sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = gi * ua18.d + t3 * ub18.d;
if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
if ((y = fi + (t2 - fi * u17.d)) == fi + (t2 + fi * u17.d))
{
retval = (sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = fi * ua17.d + t3 * ub17.d;
if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
{
retval = (sy * y);
goto ret;
}
}
/* Second stage */
ffi = xfg[i][3].d;
EADD (z0, yya, z, zz);
MUL2 (z, zz, z, zz, z2, zz2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8);
SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
if (n)
{
/* -cot */
DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u20.d * c3)) == c3 + (cc3 + u20.d * c3))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u19.d * c3)) == c3 + (cc3 + u19.d * c3))
{
retval = (sy * y);
goto ret;
}
}
retval = tanMp (x);
goto ret;
}
/* (---) The case 1e8 < abs(x) < 2**1024 */
/* Range reduction by algorithm iii */
n = (__branred (x, &a, &da)) & 0x00000001;
EADD (a, da, t1, t2);
a = t1;
da = t2;
if (a < 0.0)
{
ya = -a;
yya = -da;
sy = -1;
}
else
{
ya = a;
yya = da;
sy = 1;
}
/* (+++) The case 1e8 < abs(x) < 2**1024, abs(y) <= 1e-7 */
if (ya <= gy1.d)
{
retval = tanMp (x);
goto ret;
}
/* (X) The case 1e8 < abs(x) < 2**1024, 1e-7 < abs(y) <= 0.0608 */
if (ya <= gy2.d)
{
a2 = a * a;
t2 = d9.d + a2 * d11.d;
t2 = d7.d + a2 * t2;
t2 = d5.d + a2 * t2;
t2 = d3.d + a2 * t2;
t2 = da + a * a2 * t2;
if (n)
{
/* First stage -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c + (dc - u22.d * c)) == c + (dc + u22.d * c))
{
retval = (-y);
goto ret;
}
}
else
{
/* First stage tan */
if ((y = a + (t2 - u21.d * a)) == a + (t2 + u21.d * a))
{
retval = y;
goto ret;
}
}
/* Second stage */
/* Reduction by algorithm iv */
p = 10;
n = (__mpranred (x, &mpa, p)) & 0x00000001;
__mp_dbl (&mpa, &a, p);
__dbl_mp (a, &mpt1, p);
__sub (&mpa, &mpt1, &mpt2, p);
__mp_dbl (&mpt2, &da, p);
MUL2 (a, da, a, da, x2, xx2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = a25.d + x2 * a27.d;
c1 = a23.d + x2 * c1;
c1 = a21.d + x2 * c1;
c1 = a19.d + x2 * c1;
c1 = a17.d + x2 * c1;
c1 = a15.d + x2 * c1;
c1 *= x2;
ADD2 (a13.d, aa13.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a11.d, aa11.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a9.d, aa9.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a7.d, aa7.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a5.d, aa5.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (x2, xx2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (a, da, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a, da, c2, cc2, c1, cc1, t1, t2);
if (n)
{
/* Second stage -cot */
DIV2 (1.0, 0.0, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8,
t9, t10);
if ((y = c2 + (cc2 - u24.d * c2)) == c2 + (cc2 + u24.d * c2))
{
retval = (-y);
goto ret;
}
}
else
{
/* Second stage tan */
if ((y = c1 + (cc1 - u23.d * c1)) == c1 + (cc1 + u23.d * c1))
{
retval = y;
goto ret;
}
}
retval = tanMp (x);
goto ret;
}
/* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
/* First stage */
i = ((int) (mfftnhf.d + TWO8 * ya));
z = (z0 = (ya - xfg[i][0].d)) + yya;
z2 = z * z;
pz = z + z * z2 * (e0.d + z2 * e1.d);
fi = xfg[i][1].d;
gi = xfg[i][2].d;
if (n)
{
/* -cot */
t2 = pz * (fi + gi) / (fi + pz);
if ((y = gi - (t2 - gi * u26.d)) == gi - (t2 + gi * u26.d))
{
retval = (-sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = gi * ua26.d + t3 * ub26.d;
if ((y = gi - (t2 - t4)) == gi - (t2 + t4))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
if ((y = fi + (t2 - fi * u25.d)) == fi + (t2 + fi * u25.d))
{
retval = (sy * y);
goto ret;
}
t3 = (t2 < 0.0) ? -t2 : t2;
t4 = fi * ua25.d + t3 * ub25.d;
if ((y = fi + (t2 - t4)) == fi + (t2 + t4))
{
retval = (sy * y);
goto ret;
}
}
/* Second stage */
ffi = xfg[i][3].d;
EADD (z0, yya, z, zz);
MUL2 (z, zz, z, zz, z2, zz2, t1, t2, t3, t4, t5, t6, t7, t8);
c1 = z2 * (a7.d + z2 * (a9.d + z2 * a11.d));
ADD2 (a5.d, aa5.d, c1, 0.0, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (a3.d, aa3.d, c1, cc1, c2, cc2, t1, t2);
MUL2 (z2, zz2, c2, cc2, c1, cc1, t1, t2, t3, t4, t5, t6, t7, t8);
MUL2 (z, zz, c1, cc1, c2, cc2, t1, t2, t3, t4, t5, t6, t7, t8);
ADD2 (z, zz, c2, cc2, c1, cc1, t1, t2);
ADD2 (fi, ffi, c1, cc1, c2, cc2, t1, t2);
MUL2 (fi, ffi, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8);
SUB2 (1.0, 0.0, c3, cc3, c1, cc1, t1, t2);
if (n)
{
/* -cot */
DIV2 (c1, cc1, c2, cc2, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u28.d * c3)) == c3 + (cc3 + u28.d * c3))
{
retval = (-sy * y);
goto ret;
}
}
else
{
/* tan */
DIV2 (c2, cc2, c1, cc1, c3, cc3, t1, t2, t3, t4, t5, t6, t7, t8, t9,
t10);
if ((y = c3 + (cc3 - u27.d * c3)) == c3 + (cc3 + u27.d * c3))
{
retval = (sy * y);
goto ret;
}
}
retval = tanMp (x);
goto ret;
ret:
return retval;
}
/* multiple precision stage */
/* Convert x to multi precision number,compute tan(x) by mptan() routine */
/* and converts result back to double */
static double
SECTION
tanMp (double x)
{
int p;
double y;
mp_no mpy;
p = 32;
__mptan (x, &mpy, p);
__mp_dbl (&mpy, &y, p);
LIBC_PROBE (slowtan, 2, &x, &y);
return y;
}
#ifndef __tan
libm_alias_double (__tan, tan)
#endif
Reference in New Issue View Git Blame Copy Permalink