glibc/sysdeps/ieee754/dbl-64/s_sin.c
Siddhesh Poyarekar 9d84d0e51d Use fabs(x) instead of branching on signedness of input to sin and cos
The sin and cos code is inconsistent about its use of fabs to get the
absolute value of X where in some places it conditionalizes the code
while in others it uses fabs.  fabs seems to be a better candidate in
most cases because it avoids a branch.  Similarly there is an attempt
to make it easier for the compiler to emit conditional assignment
instructions (like fcsel on aarch64) where it can, by isolating
conditional assignment constructs from the rest of the expression.

A further benefit of this change is to identify common constructs
across functions and consolidate them in future patches.

	* sysdeps/ieee754/dbl-64/s_sin.c (do_cos_slow): Use ternary
	instead of if/else.
	(do_sin_slow): Likewise.
	(do_sincos_1): Use fabs instead of if/else.
	(do_sincos_2): Likewise.
	(__sin): Likewise.
	(__cos): Likewise.
	(slow2): Likewise.
	(sloww): Likewise.
	(sloww1): Likewise.  Drop argument M.
	(sloww2): Use fabs instead of if/else.
	(bsloww): Likewise.
	(bsloww1): Likewise.
	(bsloww2): Likewise.
2016-08-30 13:01:59 +05:30

969 lines
28 KiB
C

/*
* IBM Accurate Mathematical Library
* written by International Business Machines Corp.
* Copyright (C) 2001-2016 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:usncs.c */
/* */
/* FUNCTIONS: usin */
/* ucos */
/* slow */
/* slow1 */
/* slow2 */
/* sloww */
/* sloww1 */
/* sloww2 */
/* bsloww */
/* bsloww1 */
/* bsloww2 */
/* cslow2 */
/* FILES NEEDED: dla.h endian.h mpa.h mydefs.h usncs.h */
/* branred.c sincos32.c dosincos.c mpa.c */
/* sincos.tbl */
/* */
/* An ultimate sin and routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) or cos(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 "mydefs.h"
#include "usncs.h"
#include "MathLib.h"
#include <math.h>
#include <math_private.h>
#include <fenv.h>
/* Helper macros to compute sin of the input values. */
#define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
#define POLYNOMIAL(xx) (POLYNOMIAL2 (xx) + s1)
/* The computed polynomial is a variation of the Taylor series expansion for
sin(a):
a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
on. The result is returned to LHS and correction in COR. */
#define TAYLOR_SIN(xx, a, da, cor) \
({ \
double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \
double res = (a) + t; \
(cor) = ((a) - res) + t; \
res; \
})
/* This is again a variation of the Taylor series expansion with the term
x^3/3! expanded into the following for better accuracy:
bb * x ^ 3 + 3 * aa * x * x1 * x2 + aa * x1 ^ 3 + aa * x2 ^ 3
The correction term is dx and bb + aa = -1/3!
*/
#define TAYLOR_SLOW(x0, dx, cor) \
({ \
static const double th2_36 = 206158430208.0; /* 1.5*2**37 */ \
double xx = (x0) * (x0); \
double x1 = ((x0) + th2_36) - th2_36; \
double y = aa * x1 * x1 * x1; \
double r = (x0) + y; \
double x2 = ((x0) - x1) + (dx); \
double t = (((POLYNOMIAL2 (xx) + bb) * xx + 3.0 * aa * x1 * x2) \
* (x0) + aa * x2 * x2 * x2 + (dx)); \
t = (((x0) - r) + y) + t; \
double res = r + t; \
(cor) = (r - res) + t; \
res; \
})
#define SINCOS_TABLE_LOOKUP(u, sn, ssn, cs, ccs) \
({ \
int4 k = u.i[LOW_HALF] << 2; \
sn = __sincostab.x[k]; \
ssn = __sincostab.x[k + 1]; \
cs = __sincostab.x[k + 2]; \
ccs = __sincostab.x[k + 3]; \
})
#ifndef SECTION
# define SECTION
#endif
extern const union
{
int4 i[880];
double x[440];
} __sincostab attribute_hidden;
static const double
sn3 = -1.66666666666664880952546298448555E-01,
sn5 = 8.33333214285722277379541354343671E-03,
cs2 = 4.99999999999999999999950396842453E-01,
cs4 = -4.16666666666664434524222570944589E-02,
cs6 = 1.38888874007937613028114285595617E-03;
static const double t22 = 0x1.8p22;
void __dubsin (double x, double dx, double w[]);
void __docos (double x, double dx, double w[]);
double __mpsin (double x, double dx, bool reduce_range);
double __mpcos (double x, double dx, bool reduce_range);
static double slow (double x);
static double slow1 (double x);
static double slow2 (double x);
static double sloww (double x, double dx, double orig, int n);
static double sloww1 (double x, double dx, double orig, int n);
static double sloww2 (double x, double dx, double orig, int n);
static double bsloww (double x, double dx, double orig, int n);
static double bsloww1 (double x, double dx, double orig, int n);
static double bsloww2 (double x, double dx, double orig, int n);
int __branred (double x, double *a, double *aa);
static double cslow2 (double x);
/* Given a number partitioned into U and X such that U is an index into the
sin/cos table, this macro computes the cosine of the number by combining
the sin and cos of X (as computed by a variation of the Taylor series) with
the values looked up from the sin/cos table to get the result in RES and a
correction value in COR. */
static double
do_cos (mynumber u, double x, double *corp)
{
double xx, s, sn, ssn, c, cs, ccs, res, cor;
xx = x * x;
s = x + x * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
cor = (ccs - s * ssn - cs * c) - sn * s;
res = cs + cor;
cor = (cs - res) + cor;
*corp = cor;
return res;
}
/* A more precise variant of DO_COS where the number is partitioned into U, X
and DX. EPS is the adjustment to the correction COR. */
static double
do_cos_slow (mynumber u, double x, double dx, double eps, double *corp)
{
double xx, y, x1, x2, e1, e2, res, cor;
double s, sn, ssn, c, cs, ccs;
xx = x * x;
s = x * xx * (sn3 + xx * sn5);
c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
x1 = (x + t22) - t22;
x2 = (x - x1) + dx;
e1 = (sn + t22) - t22;
e2 = (sn - e1) + ssn;
cor = (ccs - cs * c - e1 * x2 - e2 * x) - sn * s;
y = cs - e1 * x1;
cor = cor + ((cs - y) - e1 * x1);
res = y + cor;
cor = (y - res) + cor;
cor = 1.0005 * cor + ((cor > 0) ? eps : -eps);
*corp = cor;
return res;
}
/* Given a number partitioned into U and X and DX such that U is an index into
the sin/cos table, this macro computes the sine of the number by combining
the sin and cos of X (as computed by a variation of the Taylor series) with
the values looked up from the sin/cos table to get the result in RES and a
correction value in COR. */
static double
do_sin (mynumber u, double x, double dx, double *corp)
{
double xx, s, sn, ssn, c, cs, ccs, cor, res;
xx = x * x;
s = x + (dx + x * xx * (sn3 + xx * sn5));
c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
cor = (ssn + s * ccs - sn * c) + cs * s;
res = sn + cor;
cor = (sn - res) + cor;
*corp = cor;
return res;
}
/* A more precise variant of res = do_sin where the number is partitioned into U, X
and DX. EPS is the adjustment to the correction COR. */
static double
do_sin_slow (mynumber u, double x, double dx, double eps, double *corp)
{
double xx, y, x1, x2, c1, c2, res, cor;
double s, sn, ssn, c, cs, ccs;
xx = x * x;
s = x * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
x1 = (x + t22) - t22;
x2 = (x - x1) + dx;
c1 = (cs + t22) - t22;
c2 = (cs - c1) + ccs;
cor = (ssn + s * ccs + cs * s + c2 * x + c1 * x2 - sn * x * dx) - sn * c;
y = sn + c1 * x1;
cor = cor + ((sn - y) + c1 * x1);
res = y + cor;
cor = (y - res) + cor;
cor = 1.0005 * cor + ((cor > 0) ? eps : -eps);
*corp = cor;
return res;
}
/* Reduce range of X and compute sin of a + da. K is the amount by which to
rotate the quadrants. This allows us to use the same routine to compute cos
by simply rotating the quadrants by 1. */
static inline double
__always_inline
reduce_and_compute (double x, unsigned int k)
{
double retval = 0, a, da;
unsigned int n = __branred (x, &a, &da);
k = (n + k) % 4;
switch (k)
{
case 2:
a = -a;
da = -da;
/* Fall through. */
case 0:
if (a * a < 0.01588)
retval = bsloww (a, da, x, n);
else
retval = bsloww1 (a, da, x, n);
break;
case 1:
case 3:
retval = bsloww2 (a, da, x, n);
break;
}
return retval;
}
static inline int4
__always_inline
reduce_sincos_1 (double x, double *a, double *da)
{
mynumber v;
double t = (x * hpinv + toint);
double xn = t - toint;
v.x = t;
double y = (x - xn * mp1) - xn * mp2;
int4 n = v.i[LOW_HALF] & 3;
double db = xn * mp3;
double b = y - db;
db = (y - b) - db;
*a = b;
*da = db;
return n;
}
/* Compute sin (A + DA). cos can be computed by shifting the quadrant N
clockwise. */
static double
__always_inline
do_sincos_1 (double a, double da, double x, int4 n, int4 k)
{
double xx, retval, res, cor, y;
mynumber u;
double eps = fabs (x) * 1.2e-30;
int k1 = (n + k) & 3;
switch (k1)
{ /* quarter of unit circle */
case 2:
a = -a;
da = -da;
/* Fall through. */
case 0:
xx = a * a;
if (xx < 0.01588)
{
/* Taylor series. */
res = TAYLOR_SIN (xx, a, da, cor);
cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
retval = (res == res + cor) ? res : sloww (a, da, x, k);
}
else
{
double db = (a > 0 ? da : -da);
u.x = big + fabs (a);
y = fabs (a) - (u.x - big);
res = do_sin (u, y, db, &cor);
cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: sloww1 (a, da, x, k));
}
break;
case 1:
case 3:
{
double db = (a > 0 ? da : -da);
u.x = big + fabs (a);
y = fabs (a) - (u.x - big) + db;
res = do_cos (u, y, &cor);
cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
retval = ((res == res + cor) ? ((k1 & 2) ? -res : res)
: sloww2 (a, da, x, n));
break;
}
}
return retval;
}
static inline int4
__always_inline
reduce_sincos_2 (double x, double *a, double *da)
{
mynumber v;
double t = (x * hpinv + toint);
double xn = t - toint;
v.x = t;
double xn1 = (xn + 8.0e22) - 8.0e22;
double xn2 = xn - xn1;
double y = ((((x - xn1 * mp1) - xn1 * mp2) - xn2 * mp1) - xn2 * mp2);
int4 n = v.i[LOW_HALF] & 3;
double db = xn1 * pp3;
t = y - db;
db = (y - t) - db;
db = (db - xn2 * pp3) - xn * pp4;
double b = t + db;
db = (t - b) + db;
*a = b;
*da = db;
return n;
}
/* Compute sin (A + DA). cos can be computed by shifting the quadrant N
clockwise. */
static double
__always_inline
do_sincos_2 (double a, double da, double x, int4 n, int4 k)
{
double res, retval, cor, xx;
mynumber u;
double eps = 1.0e-24;
k = (n + k) & 3;
switch (k)
{
case 2:
a = -a;
da = -da;
/* Fall through. */
case 0:
xx = a * a;
if (xx < 0.01588)
{
/* Taylor series. */
res = TAYLOR_SIN (xx, a, da, cor);
cor = (cor > 0) ? 1.02 * cor + eps : 1.02 * cor - eps;
retval = (res == res + cor) ? res : bsloww (a, da, x, n);
}
else
{
double db = (a > 0 ? da : -da);
u.x = big + fabs (a);
double y = fabs (a) - (u.x - big);
res = do_sin (u, y, db, &cor);
cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: bsloww1 (a, da, x, n));
}
break;
case 1:
case 3:
{
double db = (a > 0 ? da : -da);
u.x = big + fabs (a);
double y = fabs (a) - (u.x - big) + db;
res = do_cos (u, y, &cor);
cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
retval = ((res == res + cor) ? ((n & 2) ? -res : res)
: bsloww2 (a, da, x, n));
break;
}
}
return retval;
}
/*******************************************************************/
/* An ultimate sin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) */
/*******************************************************************/
#ifdef IN_SINCOS
static double
#else
double
SECTION
#endif
__sin (double x)
{
double xx, res, t, cor, y, s, c, sn, ssn, cs, ccs;
mynumber u;
int4 k, m;
double retval = 0;
#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
#endif
u.x = x;
m = u.i[HIGH_HALF];
k = 0x7fffffff & m; /* no sign */
if (k < 0x3e500000) /* if x->0 =>sin(x)=x */
{
math_check_force_underflow (x);
retval = x;
}
/*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
else if (k < 0x3fd00000)
{
xx = x * x;
/* Taylor series. */
t = POLYNOMIAL (xx) * (xx * x);
res = x + t;
cor = (x - res) + t;
retval = (res == res + 1.07 * cor) ? res : slow (x);
} /* else if (k < 0x3fd00000) */
/*---------------------------- 0.25<|x|< 0.855469---------------------- */
else if (k < 0x3feb6000)
{
u.x = big + fabs (x);
y = fabs (x) - (u.x - big);
y = (x > 0 ? y : -y);
xx = y * y;
s = y + y * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
if (m <= 0)
{
sn = -sn;
ssn = -ssn;
}
cor = (ssn + s * ccs - sn * c) + cs * s;
res = sn + cor;
cor = (sn - res) + cor;
retval = (res == res + 1.096 * cor) ? res : slow1 (x);
} /* else if (k < 0x3feb6000) */
/*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
else if (k < 0x400368fd)
{
t = hp0 - fabs (x);
u.x = big + fabs (t);
y = fabs (t) - (u.x - big);
y = ((t >= 0) ? hp1 : -hp1) + y;
res = do_cos (u, y, &cor);
retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x);
} /* else if (k < 0x400368fd) */
#ifndef IN_SINCOS
/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
else if (k < 0x419921FB)
{
double a, da;
int4 n = reduce_sincos_1 (x, &a, &da);
retval = do_sincos_1 (a, da, x, n, 0);
} /* else if (k < 0x419921FB ) */
/*---------------------105414350 <|x|< 281474976710656 --------------------*/
else if (k < 0x42F00000)
{
double a, da;
int4 n = reduce_sincos_2 (x, &a, &da);
retval = do_sincos_2 (a, da, x, n, 0);
} /* else if (k < 0x42F00000 ) */
/* -----------------281474976710656 <|x| <2^1024----------------------------*/
else if (k < 0x7ff00000)
retval = reduce_and_compute (x, 0);
/*--------------------- |x| > 2^1024 ----------------------------------*/
else
{
if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
__set_errno (EDOM);
retval = x / x;
}
#endif
return retval;
}
/*******************************************************************/
/* An ultimate cos routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of cos(x) */
/*******************************************************************/
#ifdef IN_SINCOS
static double
#else
double
SECTION
#endif
__cos (double x)
{
double y, xx, res, cor, a, da;
mynumber u;
int4 k, m;
double retval = 0;
#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
#endif
u.x = x;
m = u.i[HIGH_HALF];
k = 0x7fffffff & m;
/* |x|<2^-27 => cos(x)=1 */
if (k < 0x3e400000)
retval = 1.0;
else if (k < 0x3feb6000)
{ /* 2^-27 < |x| < 0.855469 */
y = fabs (x);
u.x = big + y;
y = y - (u.x - big);
res = do_cos (u, y, &cor);
retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
} /* else if (k < 0x3feb6000) */
else if (k < 0x400368fd)
{ /* 0.855469 <|x|<2.426265 */ ;
y = hp0 - fabs (x);
a = y + hp1;
da = (y - a) + hp1;
xx = a * a;
if (xx < 0.01588)
{
res = TAYLOR_SIN (xx, a, da, cor);
cor = (cor > 0) ? 1.02 * cor + 1.0e-31 : 1.02 * cor - 1.0e-31;
retval = (res == res + cor) ? res : sloww (a, da, x, 1);
}
else
{
double db = (a > 0 ? da : -da);
u.x = big + fabs (a);
y = fabs (a) - (u.x - big);
res = do_sin (u, y, db, &cor);
cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: sloww1 (a, da, x, 1));
}
} /* else if (k < 0x400368fd) */
#ifndef IN_SINCOS
else if (k < 0x419921FB)
{ /* 2.426265<|x|< 105414350 */
double a, da;
int4 n = reduce_sincos_1 (x, &a, &da);
retval = do_sincos_1 (a, da, x, n, 1);
} /* else if (k < 0x419921FB ) */
else if (k < 0x42F00000)
{
double a, da;
int4 n = reduce_sincos_2 (x, &a, &da);
retval = do_sincos_2 (a, da, x, n, 1);
} /* else if (k < 0x42F00000 ) */
/* 281474976710656 <|x| <2^1024 */
else if (k < 0x7ff00000)
retval = reduce_and_compute (x, 1);
else
{
if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
__set_errno (EDOM);
retval = x / x; /* |x| > 2^1024 */
}
#endif
return retval;
}
/************************************************************************/
/* Routine compute sin(x) for 2^-26 < |x|< 0.25 by Taylor with more */
/* precision and if still doesn't accurate enough by mpsin or dubsin */
/************************************************************************/
static double
SECTION
slow (double x)
{
double res, cor, w[2];
res = TAYLOR_SLOW (x, 0, cor);
if (res == res + 1.0007 * cor)
return res;
__dubsin (fabs (x), 0, w);
if (w[0] == w[0] + 1.000000001 * w[1])
return (x > 0) ? w[0] : -w[0];
return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
}
/*******************************************************************************/
/* Routine compute sin(x) for 0.25<|x|< 0.855469 by __sincostab.tbl and Taylor */
/* and if result still doesn't accurate enough by mpsin or dubsin */
/*******************************************************************************/
static double
SECTION
slow1 (double x)
{
mynumber u;
double w[2], y, cor, res;
y = fabs (x);
u.x = big + y;
y = y - (u.x - big);
res = do_sin_slow (u, y, 0, 0, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
__dubsin (fabs (x), 0, w);
if (w[0] == w[0] + 1.000000005 * w[1])
return (x > 0) ? w[0] : -w[0];
return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
}
/**************************************************************************/
/* Routine compute sin(x) for 0.855469 <|x|<2.426265 by __sincostab.tbl */
/* and if result still doesn't accurate enough by mpsin or dubsin */
/**************************************************************************/
static double
SECTION
slow2 (double x)
{
mynumber u;
double w[2], y, y1, y2, cor, res, del;
double t = hp0 - fabs (x);
u.x = big + fabs (t);
y = fabs (t) - (u.x - big);
del = (t >= 0) ? hp1 : -hp1;
res = do_cos_slow (u, y, del, 0, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
y = fabs (x) - hp0;
y1 = y - hp1;
y2 = (y - y1) - hp1;
__docos (y1, y2, w);
if (w[0] == w[0] + 1.000000005 * w[1])
return (x > 0) ? w[0] : -w[0];
return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
}
/***************************************************************************/
/* Routine compute sin(x+dx) (Double-Length number) where x is small enough*/
/* to use Taylor series around zero and (x+dx) */
/* in first or third quarter of unit circle.Routine receive also */
/* (right argument) the original value of x for computing error of */
/* result.And if result not accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
static double
SECTION
sloww (double x, double dx, double orig, int k)
{
double y, t, res, cor, w[2], a, da, xn;
mynumber v;
int4 n;
res = TAYLOR_SLOW (x, dx, cor);
double eps = fabs (orig) * 3.1e-30;
cor = 1.0005 * cor + ((cor > 0) ? eps : -eps);
if (res == res + cor)
return res;
a = fabs (x);
da = (x > 0) ? dx : -dx;
__dubsin (a, da, w);
eps = fabs (orig) * 1.1e-30;
cor = 1.000000001 * w[1] + ((w[1] > 0) ? eps : -eps);
if (w[0] == w[0] + cor)
return (x > 0) ? w[0] : -w[0];
t = (orig * hpinv + toint);
xn = t - toint;
v.x = t;
y = (orig - xn * mp1) - xn * mp2;
n = (v.i[LOW_HALF] + k) & 3;
da = xn * pp3;
t = y - da;
da = (y - t) - da;
y = xn * pp4;
a = t - y;
da = ((t - a) - y) + da;
if (n & 2)
{
a = -a;
da = -da;
}
x = fabs (a);
dx = (a > 0) ? da : -da;
__dubsin (x, dx, w);
eps = fabs (orig) * 1.1e-40;
cor = 1.000000001 * w[1] + ((w[1] > 0) ? eps : -eps);
if (w[0] == w[0] + cor)
return (a > 0) ? w[0] : -w[0];
return k ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
}
/***************************************************************************/
/* Routine compute sin(x+dx) (Double-Length number) where x in first or */
/* third quarter of unit circle.Routine receive also (right argument) the */
/* original value of x for computing error of result.And if result not */
/* accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
static double
SECTION
sloww1 (double x, double dx, double orig, int k)
{
mynumber u;
double w[2], y, cor, res;
u.x = big + fabs (x);
y = fabs (x) - (u.x - big);
dx = (x > 0 ? dx : -dx);
res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
__dubsin (fabs (x), dx, w);
double eps = 1.1e-30 * fabs (orig);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? eps : -eps);
if (w[0] == w[0] + cor)
return (x > 0) ? w[0] : -w[0];
return (k == 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
}
/***************************************************************************/
/* Routine compute sin(x+dx) (Double-Length number) where x in second or */
/* fourth quarter of unit circle.Routine receive also the original value */
/* and quarter(n= 1or 3)of x for computing error of result.And if result not*/
/* accurate enough routine calls mpsin1 or dubsin */
/***************************************************************************/
static double
SECTION
sloww2 (double x, double dx, double orig, int n)
{
mynumber u;
double w[2], y, cor, res;
u.x = big + fabs (x);
y = fabs (x) - (u.x - big);
dx = (x > 0 ? dx : -dx);
res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
if (res == res + cor)
return (n & 2) ? -res : res;
__docos (fabs (x), dx, w);
double eps = 1.1e-30 * fabs (orig);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? eps : -eps);
if (w[0] == w[0] + cor)
return (n & 2) ? -w[0] : w[0];
return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
}
/***************************************************************************/
/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
/* is small enough to use Taylor series around zero and (x+dx) */
/* in first or third quarter of unit circle.Routine receive also */
/* (right argument) the original value of x for computing error of */
/* result.And if result not accurate enough routine calls other routines */
/***************************************************************************/
static double
SECTION
bsloww (double x, double dx, double orig, int n)
{
double res, cor, w[2], a, da;
res = TAYLOR_SLOW (x, dx, cor);
cor = 1.0005 * cor + ((cor > 0) ? 1.1e-24 : -1.1e-24);
if (res == res + cor)
return res;
a = fabs (x);
da = (x > 0) ? dx : -dx;
__dubsin (a, da, w);
cor = 1.000000001 * w[1] + ((w[1] > 0) ? 1.1e-24 : -1.1e-24);
if (w[0] == w[0] + cor)
return (x > 0) ? w[0] : -w[0];
return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
}
/***************************************************************************/
/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
/* in first or third quarter of unit circle.Routine receive also */
/* (right argument) the original value of x for computing error of result.*/
/* And if result not accurate enough routine calls other routines */
/***************************************************************************/
static double
SECTION
bsloww1 (double x, double dx, double orig, int n)
{
mynumber u;
double w[2], y, cor, res;
y = fabs (x);
u.x = big + y;
y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
res = do_sin_slow (u, y, dx, 1.1e-24, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
__dubsin (fabs (x), dx, w);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? 1.1e-24 : -1.1e-24);
if (w[0] == w[0] + cor)
return (x > 0) ? w[0] : -w[0];
return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
}
/***************************************************************************/
/* Routine compute sin(x+dx) or cos(x+dx) (Double-Length number) where x */
/* in second or fourth quarter of unit circle.Routine receive also the */
/* original value and quarter(n= 1or 3)of x for computing error of result. */
/* And if result not accurate enough routine calls other routines */
/***************************************************************************/
static double
SECTION
bsloww2 (double x, double dx, double orig, int n)
{
mynumber u;
double w[2], y, cor, res;
y = fabs (x);
u.x = big + y;
y = y - (u.x - big);
dx = (x > 0) ? dx : -dx;
res = do_cos_slow (u, y, dx, 1.1e-24, &cor);
if (res == res + cor)
return (n & 2) ? -res : res;
__docos (fabs (x), dx, w);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? 1.1e-24 : -1.1e-24);
if (w[0] == w[0] + cor)
return (n & 2) ? -w[0] : w[0];
return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
}
/************************************************************************/
/* Routine compute cos(x) for 2^-27 < |x|< 0.25 by Taylor with more */
/* precision and if still doesn't accurate enough by mpcos or docos */
/************************************************************************/
static double
SECTION
cslow2 (double x)
{
mynumber u;
double w[2], y, cor, res;
y = fabs (x);
u.x = big + y;
y = y - (u.x - big);
res = do_cos_slow (u, y, 0, 0, &cor);
if (res == res + cor)
return res;
y = fabs (x);
__docos (y, 0, w);
if (w[0] == w[0] + 1.000000005 * w[1])
return w[0];
return __mpcos (x, 0, false);
}
#ifndef __cos
weak_alias (__cos, cos)
# ifdef NO_LONG_DOUBLE
strong_alias (__cos, __cosl)
weak_alias (__cos, cosl)
# endif
#endif
#ifndef __sin
weak_alias (__sin, sin)
# ifdef NO_LONG_DOUBLE
strong_alias (__sin, __sinl)
weak_alias (__sin, sinl)
# endif
#endif