AuroraMiddleware/glibc
1
0
Fork 0
mirror of https://sourceware.org/git/glibc.git synced 2024-11-25 14:30:06 +00:00
Code Issues Packages Projects Releases Wiki Activity
3eab466e56
glibc/sysdeps/ieee754/dbl-64/s_tan.c
Wilco Dijkstra 476d692e8a math: Remove slow paths in tan [BZ #15267] 
Remove slow paths in tan. Add ULP annotations. Merge 'number' into 'mynumber'.
Remove unused entries from tan constants.

Reviewed-By: Paul Zimmermann <Paul.Zimmermann@inria.fr>
2021-03-11 14:26:36 +00:00

376 lines
8.0 KiB
C
Raw Blame History

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001-2021 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 <https://www.gnu.org/licenses/>.
*/
/*********************************************************************/
/* MODULE_NAME: utan.c */
/* */
/* FUNCTIONS: utan */
/* */
/* FILES NEEDED:dla.h endian.h mydefs.h utan.h */
/* branred.c */
/* utan.tbl */
/* */
/*********************************************************************/
#include <errno.h>
#include <float.h>
#include "endian.h"
#include <dla.h>
#include "mydefs.h"
#include <math.h>
#include <math_private.h>
#include <fenv_private.h>
#include <math-underflow.h>
#include <libm-alias-double.h>
#include <fenv.h>
#ifndef SECTION
# define SECTION
#endif
/* tan with max ULP of ~0.619 based on random sampling. */
double
SECTION
__tan (double x)
{
#include "utan.h"
#include "utan.tbl"
int ux, i, n;
double a, da, a2, b, db, c, dc, fi, gi, pz,
s, sy, t, t1, t2, t3, t4, w, x2, xn, y, ya,
yya, z0, z, z2;
mynumber num, v;
double retval;
int __branred (double, double *, double *);
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)
{
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;
y = x + t2;
retval = y;
/* Max ULP is 0.504. */
goto ret;
}
/* (III) The case 0.0608 < abs(x) <= 0.787 */
if (w <= g3.d)
{
i = ((int) (mfftnhf.d + 256 * 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);
y = fi + t2;
retval = (s * y);
/* Max ULP is 0.60. */
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;
}
/* (VI) The case 0.787 < abs(x) <= 25, 0 < 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)
{
/* -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
y = c + dc;
retval = (-y);
/* Max ULP is 0.506. */
goto ret;
}
else
{
/* tan */
y = a + t2;
retval = y;
/* Max ULP is 0.506. */
goto ret;
}
}
/* (VII) The case 0.787 < abs(x) <= 25, 0.0608 < abs(y) <= 0.787 */
i = ((int) (mfftnhf.d + 256 * 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);
y = gi - t2;
retval = (-sy * y);
/* Max ULP is 0.62. */
goto ret;
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
y = fi + t2;
retval = (sy * y);
/* Max ULP is 0.62. */
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;
}
/* (VIII) The case 25 < abs(x) <= 1e8, 0 < 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)
{
/* -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
y = c + dc;
retval = (-y);
/* Max ULP is 0.506. */
goto ret;
}
else
{
/* tan */
y = a + t2;
retval = y;
/* Max ULP is 0.506. */
goto ret;
}
}
/* (IX) The case 25 < abs(x) <= 1e8, 0.0608 < abs(y) <= 0.787 */
i = ((int) (mfftnhf.d + 256 * 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);
y = gi - t2;
retval = (-sy * y);
/* Max ULP is 0.62. */
goto ret;
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
y = fi + t2;
retval = (sy * y);
/* Max ULP is 0.62. */
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;
}
/* (X) The case 1e8 < abs(x) < 2**1024, 0 < 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)
{
/* -cot */
EADD (a, t2, b, db);
DIV2 (1.0, 0.0, b, db, c, dc, t1, t2, t3, t4);
y = c + dc;
retval = (-y);
/* Max ULP is 0.506. */
goto ret;
}
else
{
/* tan */
y = a + t2;
retval = y;
/* Max ULP is 0.507. */
goto ret;
}
}
/* (XI) The case 1e8 < abs(x) < 2**1024, 0.0608 < abs(y) <= 0.787 */
i = ((int) (mfftnhf.d + 256 * 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);
y = gi - t2;
retval = (-sy * y);
/* Max ULP is 0.62. */
goto ret;
}
else
{
/* tan */
t2 = pz * (gi + fi) / (gi - pz);
y = fi + t2;
retval = (sy * y);
/* Max ULP is 0.62. */
goto ret;
}
ret:
return retval;
}
#ifndef __tan
libm_alias_double (__tan, tan)
#endif
Reference in New Issue View Git Blame Copy Permalink