summaryrefslogtreecommitdiffstatshomepage
path: root/src/lib/softfloat/fsincos.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/softfloat/fsincos.c')
-rwxr-xr-xsrc/lib/softfloat/fsincos.c611
1 files changed, 305 insertions, 306 deletions
diff --git a/src/lib/softfloat/fsincos.c b/src/lib/softfloat/fsincos.c
index fb148cc8fbe..84c67c47270 100755
--- a/src/lib/softfloat/fsincos.c
+++ b/src/lib/softfloat/fsincos.c
@@ -47,7 +47,7 @@ static const floatx80 floatx80_default_nan = packFloatx80(0, 0xffff, U64(0xfffff
INLINE bits64 extractFloatx80Frac( floatx80 a )
{
- return a.low;
+ return a.low;
}
@@ -59,7 +59,7 @@ INLINE bits64 extractFloatx80Frac( floatx80 a )
INLINE int32 extractFloatx80Exp( floatx80 a )
{
- return a.high & 0x7FFF;
+ return a.high & 0x7FFF;
}
@@ -71,7 +71,7 @@ INLINE int32 extractFloatx80Exp( floatx80 a )
INLINE flag extractFloatx80Sign( floatx80 a )
{
- return a.high>>15;
+ return a.high>>15;
}
@@ -83,12 +83,12 @@ INLINE flag extractFloatx80Sign( floatx80 a )
INLINE floatx80 propagateFloatx80NaNOneArg(floatx80 a)
{
- if (floatx80_is_signaling_nan(a))
- float_raise(float_flag_invalid);
+ if (floatx80_is_signaling_nan(a))
+ float_raise(float_flag_invalid);
- a.low |= U64(0xC000000000000000);
+ a.low |= U64(0xC000000000000000);
- return a;
+ return a;
}
/*----------------------------------------------------------------------------
@@ -100,62 +100,62 @@ INLINE floatx80 propagateFloatx80NaNOneArg(floatx80 a)
void normalizeFloatx80Subnormal(UINT64 aSig, INT32 *zExpPtr, UINT64 *zSigPtr)
{
- int shiftCount = countLeadingZeros64(aSig);
- *zSigPtr = aSig<<shiftCount;
- *zExpPtr = 1 - shiftCount;
+ int shiftCount = countLeadingZeros64(aSig);
+ *zSigPtr = aSig<<shiftCount;
+ *zExpPtr = 1 - shiftCount;
}
/* reduce trigonometric function argument using 128-bit precision
M_PI approximation */
static UINT64 argument_reduction_kernel(UINT64 aSig0, int Exp, UINT64 *zSig0, UINT64 *zSig1)
{
- UINT64 term0, term1, term2;
- UINT64 aSig1 = 0;
-
- shortShift128Left(aSig1, aSig0, Exp, &aSig1, &aSig0);
- UINT64 q = estimateDiv128To64(aSig1, aSig0, FLOAT_PI_HI);
- mul128By64To192(FLOAT_PI_HI, FLOAT_PI_LO, q, &term0, &term1, &term2);
- sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
- while ((INT64)(*zSig1) < 0) {
- --q;
- add192(*zSig1, *zSig0, term2, 0, FLOAT_PI_HI, FLOAT_PI_LO, zSig1, zSig0, &term2);
- }
- *zSig1 = term2;
- return q;
+ UINT64 term0, term1, term2;
+ UINT64 aSig1 = 0;
+
+ shortShift128Left(aSig1, aSig0, Exp, &aSig1, &aSig0);
+ UINT64 q = estimateDiv128To64(aSig1, aSig0, FLOAT_PI_HI);
+ mul128By64To192(FLOAT_PI_HI, FLOAT_PI_LO, q, &term0, &term1, &term2);
+ sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
+ while ((INT64)(*zSig1) < 0) {
+ --q;
+ add192(*zSig1, *zSig0, term2, 0, FLOAT_PI_HI, FLOAT_PI_LO, zSig1, zSig0, &term2);
+ }
+ *zSig1 = term2;
+ return q;
}
static int reduce_trig_arg(int expDiff, int &zSign, UINT64 &aSig0, UINT64 &aSig1)
{
- UINT64 term0, term1, q = 0;
-
- if (expDiff < 0) {
- shift128Right(aSig0, 0, 1, &aSig0, &aSig1);
- expDiff = 0;
- }
- if (expDiff > 0) {
- q = argument_reduction_kernel(aSig0, expDiff, &aSig0, &aSig1);
- }
- else {
- if (FLOAT_PI_HI <= aSig0) {
- aSig0 -= FLOAT_PI_HI;
- q = 1;
- }
- }
-
- shift128Right(FLOAT_PI_HI, FLOAT_PI_LO, 1, &term0, &term1);
- if (! lt128(aSig0, aSig1, term0, term1))
- {
- int lt = lt128(term0, term1, aSig0, aSig1);
- int eq = eq128(aSig0, aSig1, term0, term1);
-
- if ((eq && (q & 1)) || lt) {
- zSign = !zSign;
- ++q;
- }
- if (lt) sub128(FLOAT_PI_HI, FLOAT_PI_LO, aSig0, aSig1, &aSig0, &aSig1);
- }
-
- return (int)(q & 3);
+ UINT64 term0, term1, q = 0;
+
+ if (expDiff < 0) {
+ shift128Right(aSig0, 0, 1, &aSig0, &aSig1);
+ expDiff = 0;
+ }
+ if (expDiff > 0) {
+ q = argument_reduction_kernel(aSig0, expDiff, &aSig0, &aSig1);
+ }
+ else {
+ if (FLOAT_PI_HI <= aSig0) {
+ aSig0 -= FLOAT_PI_HI;
+ q = 1;
+ }
+ }
+
+ shift128Right(FLOAT_PI_HI, FLOAT_PI_LO, 1, &term0, &term1);
+ if (! lt128(aSig0, aSig1, term0, term1))
+ {
+ int lt = lt128(term0, term1, aSig0, aSig1);
+ int eq = eq128(aSig0, aSig1, term0, term1);
+
+ if ((eq && (q & 1)) || lt) {
+ zSign = !zSign;
+ ++q;
+ }
+ if (lt) sub128(FLOAT_PI_HI, FLOAT_PI_LO, aSig0, aSig1, &aSig0, &aSig1);
+ }
+
+ return (int)(q & 3);
}
#define SIN_ARR_SIZE 11
@@ -163,32 +163,32 @@ static int reduce_trig_arg(int expDiff, int &zSign, UINT64 &aSig0, UINT64 &aSig1
static float128 sin_arr[SIN_ARR_SIZE] =
{
- PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
- PACK_FLOAT_128(0xbffc555555555555, 0x5555555555555555), /* 3 */
- PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
- PACK_FLOAT_128(0xbff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
- PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
- PACK_FLOAT_128(0xbfe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
- PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
- PACK_FLOAT_128(0xbfd6ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 15 */
- PACK_FLOAT_128(0x3fce952c77030ad4, 0xa6b2605197771b00), /* 17 */
- PACK_FLOAT_128(0xbfc62f49b4681415, 0x724ca1ec3b7b9675), /* 19 */
- PACK_FLOAT_128(0x3fbd71b8ef6dcf57, 0x18bef146fcee6e45) /* 21 */
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
+ PACK_FLOAT_128(0xbffc555555555555, 0x5555555555555555), /* 3 */
+ PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
+ PACK_FLOAT_128(0xbff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
+ PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
+ PACK_FLOAT_128(0xbfe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
+ PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
+ PACK_FLOAT_128(0xbfd6ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 15 */
+ PACK_FLOAT_128(0x3fce952c77030ad4, 0xa6b2605197771b00), /* 17 */
+ PACK_FLOAT_128(0xbfc62f49b4681415, 0x724ca1ec3b7b9675), /* 19 */
+ PACK_FLOAT_128(0x3fbd71b8ef6dcf57, 0x18bef146fcee6e45) /* 21 */
};
static float128 cos_arr[COS_ARR_SIZE] =
{
- PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 0 */
- PACK_FLOAT_128(0xbffe000000000000, 0x0000000000000000), /* 2 */
- PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
- PACK_FLOAT_128(0xbff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
- PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
- PACK_FLOAT_128(0xbfe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
- PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
- PACK_FLOAT_128(0xbfda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
- PACK_FLOAT_128(0x3fd2ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 16 */
- PACK_FLOAT_128(0xbfca6827863b97d9, 0x77bb004886a2c2ab), /* 18 */
- PACK_FLOAT_128(0x3fc1e542ba402022, 0x507a9cad2bf8f0bb) /* 20 */
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 0 */
+ PACK_FLOAT_128(0xbffe000000000000, 0x0000000000000000), /* 2 */
+ PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
+ PACK_FLOAT_128(0xbff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
+ PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
+ PACK_FLOAT_128(0xbfe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
+ PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
+ PACK_FLOAT_128(0xbfda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
+ PACK_FLOAT_128(0x3fd2ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 16 */
+ PACK_FLOAT_128(0xbfca6827863b97d9, 0x77bb004886a2c2ab), /* 18 */
+ PACK_FLOAT_128(0x3fc1e542ba402022, 0x507a9cad2bf8f0bb) /* 20 */
};
extern float128 OddPoly (float128 x, float128 *arr, unsigned n);
@@ -196,27 +196,27 @@ extern float128 OddPoly (float128 x, float128 *arr, unsigned n);
/* 0 <= x <= pi/4 */
INLINE float128 poly_sin(float128 x)
{
- // 3 5 7 9 11 13 15
- // x x x x x x x
- // sin (x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- =
- // 3! 5! 7! 9! 11! 13! 15!
- //
- // 2 4 6 8 10 12 14
- // x x x x x x x
- // = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- ] =
- // 3! 5! 7! 9! 11! 13! 15!
- //
- // 3 3
- // -- 4k -- 4k+2
- // p(x) = > C * x > 0 q(x) = > C * x < 0
- // -- 2k -- 2k+1
- // k=0 k=0
- //
- // 2
- // sin(x) ~ x * [ p(x) + x * q(x) ]
- //
-
- return OddPoly(x, sin_arr, SIN_ARR_SIZE);
+ // 3 5 7 9 11 13 15
+ // x x x x x x x
+ // sin (x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- =
+ // 3! 5! 7! 9! 11! 13! 15!
+ //
+ // 2 4 6 8 10 12 14
+ // x x x x x x x
+ // = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- ] =
+ // 3! 5! 7! 9! 11! 13! 15!
+ //
+ // 3 3
+ // -- 4k -- 4k+2
+ // p(x) = > C * x > 0 q(x) = > C * x < 0
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // 2
+ // sin(x) ~ x * [ p(x) + x * q(x) ]
+ //
+
+ return OddPoly(x, sin_arr, SIN_ARR_SIZE);
}
extern float128 EvenPoly(float128 x, float128 *arr, unsigned n);
@@ -224,53 +224,53 @@ extern float128 EvenPoly(float128 x, float128 *arr, unsigned n);
/* 0 <= x <= pi/4 */
INLINE float128 poly_cos(float128 x)
{
- // 2 4 6 8 10 12 14
- // x x x x x x x
- // cos (x) ~ 1 - --- + --- - --- + --- - ---- + ---- - ----
- // 2! 4! 6! 8! 10! 12! 14!
- //
- // 3 3
- // -- 4k -- 4k+2
- // p(x) = > C * x > 0 q(x) = > C * x < 0
- // -- 2k -- 2k+1
- // k=0 k=0
- //
- // 2
- // cos(x) ~ [ p(x) + x * q(x) ]
- //
-
- return EvenPoly(x, cos_arr, COS_ARR_SIZE);
+ // 2 4 6 8 10 12 14
+ // x x x x x x x
+ // cos (x) ~ 1 - --- + --- - --- + --- - ---- + ---- - ----
+ // 2! 4! 6! 8! 10! 12! 14!
+ //
+ // 3 3
+ // -- 4k -- 4k+2
+ // p(x) = > C * x > 0 q(x) = > C * x < 0
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // 2
+ // cos(x) ~ [ p(x) + x * q(x) ]
+ //
+
+ return EvenPoly(x, cos_arr, COS_ARR_SIZE);
}
INLINE void sincos_invalid(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
{
- if (sin_a) *sin_a = a;
- if (cos_a) *cos_a = a;
+ if (sin_a) *sin_a = a;
+ if (cos_a) *cos_a = a;
}
INLINE void sincos_tiny_argument(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
{
- if (sin_a) *sin_a = a;
- if (cos_a) *cos_a = floatx80_one;
+ if (sin_a) *sin_a = a;
+ if (cos_a) *cos_a = floatx80_one;
}
static floatx80 sincos_approximation(int neg, float128 r, UINT64 quotient)
{
- if (quotient & 0x1) {
- r = poly_cos(r);
- neg = 0;
- } else {
- r = poly_sin(r);
- }
+ if (quotient & 0x1) {
+ r = poly_cos(r);
+ neg = 0;
+ } else {
+ r = poly_sin(r);
+ }
- floatx80 result = float128_to_floatx80(r);
- if (quotient & 0x2)
- neg = ! neg;
+ floatx80 result = float128_to_floatx80(r);
+ if (quotient & 0x2)
+ neg = ! neg;
- if (neg)
- result = floatx80_chs(result);
+ if (neg)
+ result = floatx80_chs(result);
- return result;
+ return result;
}
// =================================================
@@ -295,92 +295,92 @@ static floatx80 sincos_approximation(int neg, float128 r, UINT64 quotient)
int sf_fsincos(floatx80 a, floatx80 *sin_a, floatx80 *cos_a)
{
- UINT64 aSig0, aSig1 = 0;
- INT32 aExp, zExp, expDiff;
- int aSign, zSign;
- int q = 0;
-
- aSig0 = extractFloatx80Frac(a);
- aExp = extractFloatx80Exp(a);
- aSign = extractFloatx80Sign(a);
-
- /* invalid argument */
- if (aExp == 0x7FFF) {
- if ((UINT64) (aSig0<<1)) {
- sincos_invalid(sin_a, cos_a, propagateFloatx80NaNOneArg(a));
- return 0;
- }
-
- float_raise(float_flag_invalid);
- sincos_invalid(sin_a, cos_a, floatx80_default_nan);
- return 0;
- }
-
- if (aExp == 0) {
- if (aSig0 == 0) {
- sincos_tiny_argument(sin_a, cos_a, a);
- return 0;
- }
+ UINT64 aSig0, aSig1 = 0;
+ INT32 aExp, zExp, expDiff;
+ int aSign, zSign;
+ int q = 0;
+
+ aSig0 = extractFloatx80Frac(a);
+ aExp = extractFloatx80Exp(a);
+ aSign = extractFloatx80Sign(a);
+
+ /* invalid argument */
+ if (aExp == 0x7FFF) {
+ if ((UINT64) (aSig0<<1)) {
+ sincos_invalid(sin_a, cos_a, propagateFloatx80NaNOneArg(a));
+ return 0;
+ }
+
+ float_raise(float_flag_invalid);
+ sincos_invalid(sin_a, cos_a, floatx80_default_nan);
+ return 0;
+ }
+
+ if (aExp == 0) {
+ if (aSig0 == 0) {
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
// float_raise(float_flag_denormal);
- /* handle pseudo denormals */
- if (! (aSig0 & U64(0x8000000000000000)))
- {
- float_raise(float_flag_inexact);
- if (sin_a)
- float_raise(float_flag_underflow);
- sincos_tiny_argument(sin_a, cos_a, a);
- return 0;
- }
-
- normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
- }
-
- zSign = aSign;
- zExp = EXP_BIAS;
- expDiff = aExp - zExp;
-
- /* argument is out-of-range */
- if (expDiff >= 63)
- return -1;
-
- float_raise(float_flag_inexact);
-
- if (expDiff < -1) { // doesn't require reduction
- if (expDiff <= -68) {
- a = packFloatx80(aSign, aExp, aSig0);
- sincos_tiny_argument(sin_a, cos_a, a);
- return 0;
- }
- zExp = aExp;
- }
- else {
- q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
- }
-
- /* **************************** */
- /* argument reduction completed */
- /* **************************** */
-
- /* using float128 for approximation */
- float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
-
- if (aSign) q = -q;
- if (sin_a) *sin_a = sincos_approximation(zSign, r, q);
- if (cos_a) *cos_a = sincos_approximation(zSign, r, q+1);
-
- return 0;
+ /* handle pseudo denormals */
+ if (! (aSig0 & U64(0x8000000000000000)))
+ {
+ float_raise(float_flag_inexact);
+ if (sin_a)
+ float_raise(float_flag_underflow);
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
+
+ normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
+ }
+
+ zSign = aSign;
+ zExp = EXP_BIAS;
+ expDiff = aExp - zExp;
+
+ /* argument is out-of-range */
+ if (expDiff >= 63)
+ return -1;
+
+ float_raise(float_flag_inexact);
+
+ if (expDiff < -1) { // doesn't require reduction
+ if (expDiff <= -68) {
+ a = packFloatx80(aSign, aExp, aSig0);
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
+ zExp = aExp;
+ }
+ else {
+ q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+ }
+
+ /* **************************** */
+ /* argument reduction completed */
+ /* **************************** */
+
+ /* using float128 for approximation */
+ float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
+
+ if (aSign) q = -q;
+ if (sin_a) *sin_a = sincos_approximation(zSign, r, q);
+ if (cos_a) *cos_a = sincos_approximation(zSign, r, q+1);
+
+ return 0;
}
int floatx80_fsin(floatx80 &a)
{
- return sf_fsincos(a, &a, 0);
+ return sf_fsincos(a, &a, 0);
}
int floatx80_fcos(floatx80 &a)
{
- return sf_fsincos(a, 0, &a);
+ return sf_fsincos(a, 0, &a);
}
// =================================================
@@ -412,83 +412,83 @@ int floatx80_fcos(floatx80 &a)
int floatx80_ftan(floatx80 &a)
{
- UINT64 aSig0, aSig1 = 0;
- INT32 aExp, zExp, expDiff;
- int aSign, zSign;
- int q = 0;
-
- aSig0 = extractFloatx80Frac(a);
- aExp = extractFloatx80Exp(a);
- aSign = extractFloatx80Sign(a);
-
- /* invalid argument */
- if (aExp == 0x7FFF) {
- if ((UINT64) (aSig0<<1))
- {
- a = propagateFloatx80NaNOneArg(a);
- return 0;
- }
-
- float_raise(float_flag_invalid);
- a = floatx80_default_nan;
- return 0;
- }
-
- if (aExp == 0) {
- if (aSig0 == 0) return 0;
+ UINT64 aSig0, aSig1 = 0;
+ INT32 aExp, zExp, expDiff;
+ int aSign, zSign;
+ int q = 0;
+
+ aSig0 = extractFloatx80Frac(a);
+ aExp = extractFloatx80Exp(a);
+ aSign = extractFloatx80Sign(a);
+
+ /* invalid argument */
+ if (aExp == 0x7FFF) {
+ if ((UINT64) (aSig0<<1))
+ {
+ a = propagateFloatx80NaNOneArg(a);
+ return 0;
+ }
+
+ float_raise(float_flag_invalid);
+ a = floatx80_default_nan;
+ return 0;
+ }
+
+ if (aExp == 0) {
+ if (aSig0 == 0) return 0;
// float_raise(float_flag_denormal);
- /* handle pseudo denormals */
- if (! (aSig0 & U64(0x8000000000000000)))
- {
- float_raise(float_flag_inexact | float_flag_underflow);
- return 0;
- }
- normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
- }
-
- zSign = aSign;
- zExp = EXP_BIAS;
- expDiff = aExp - zExp;
-
- /* argument is out-of-range */
- if (expDiff >= 63)
- return -1;
-
- float_raise(float_flag_inexact);
-
- if (expDiff < -1) { // doesn't require reduction
- if (expDiff <= -68) {
- a = packFloatx80(aSign, aExp, aSig0);
- return 0;
- }
- zExp = aExp;
- }
- else {
- q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
- }
-
- /* **************************** */
- /* argument reduction completed */
- /* **************************** */
-
- /* using float128 for approximation */
- float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
-
- float128 sin_r = poly_sin(r);
- float128 cos_r = poly_cos(r);
-
- if (q & 0x1) {
- r = float128_div(cos_r, sin_r);
- zSign = ! zSign;
- } else {
- r = float128_div(sin_r, cos_r);
- }
-
- a = float128_to_floatx80(r);
- if (zSign)
- a = floatx80_chs(a);
-
- return 0;
+ /* handle pseudo denormals */
+ if (! (aSig0 & U64(0x8000000000000000)))
+ {
+ float_raise(float_flag_inexact | float_flag_underflow);
+ return 0;
+ }
+ normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
+ }
+
+ zSign = aSign;
+ zExp = EXP_BIAS;
+ expDiff = aExp - zExp;
+
+ /* argument is out-of-range */
+ if (expDiff >= 63)
+ return -1;
+
+ float_raise(float_flag_inexact);
+
+ if (expDiff < -1) { // doesn't require reduction
+ if (expDiff <= -68) {
+ a = packFloatx80(aSign, aExp, aSig0);
+ return 0;
+ }
+ zExp = aExp;
+ }
+ else {
+ q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+ }
+
+ /* **************************** */
+ /* argument reduction completed */
+ /* **************************** */
+
+ /* using float128 for approximation */
+ float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
+
+ float128 sin_r = poly_sin(r);
+ float128 cos_r = poly_cos(r);
+
+ if (q & 0x1) {
+ r = float128_div(cos_r, sin_r);
+ zSign = ! zSign;
+ } else {
+ r = float128_div(sin_r, cos_r);
+ }
+
+ a = float128_to_floatx80(r);
+ if (zSign)
+ a = floatx80_chs(a);
+
+ return 0;
}
// 2 3 4 n
@@ -504,30 +504,30 @@ int floatx80_ftan(floatx80 &a)
float128 EvalPoly(float128 x, float128 *arr, unsigned n)
{
- float128 x2 = float128_mul(x, x);
- unsigned i;
-
- assert(n > 1);
-
- float128 r1 = arr[--n];
- i = n;
- while(i >= 2) {
- r1 = float128_mul(r1, x2);
- i -= 2;
- r1 = float128_add(r1, arr[i]);
- }
- if (i) r1 = float128_mul(r1, x);
-
- float128 r2 = arr[--n];
- i = n;
- while(i >= 2) {
- r2 = float128_mul(r2, x2);
- i -= 2;
- r2 = float128_add(r2, arr[i]);
- }
- if (i) r2 = float128_mul(r2, x);
-
- return float128_add(r1, r2);
+ float128 x2 = float128_mul(x, x);
+ unsigned i;
+
+ assert(n > 1);
+
+ float128 r1 = arr[--n];
+ i = n;
+ while(i >= 2) {
+ r1 = float128_mul(r1, x2);
+ i -= 2;
+ r1 = float128_add(r1, arr[i]);
+ }
+ if (i) r1 = float128_mul(r1, x);
+
+ float128 r2 = arr[--n];
+ i = n;
+ while(i >= 2) {
+ r2 = float128_mul(r2, x2);
+ i -= 2;
+ r2 = float128_add(r2, arr[i]);
+ }
+ if (i) r2 = float128_mul(r2, x);
+
+ return float128_add(r1, r2);
}
// 2 4 6 8 2n
@@ -544,7 +544,7 @@ float128 EvalPoly(float128 x, float128 *arr, unsigned n)
float128 EvenPoly(float128 x, float128 *arr, unsigned n)
{
- return EvalPoly(float128_mul(x, x), arr, n);
+ return EvalPoly(float128_mul(x, x), arr, n);
}
// 3 5 7 9 2n+1
@@ -564,6 +564,5 @@ float128 EvenPoly(float128 x, float128 *arr, unsigned n)
float128 OddPoly(float128 x, float128 *arr, unsigned n)
{
- return float128_mul(x, EvenPoly(x, arr, n));
+ return float128_mul(x, EvenPoly(x, arr, n));
}
-