diff options
Diffstat (limited to 'arch/x86/math-emu/poly_tan.c')
| -rw-r--r-- | arch/x86/math-emu/poly_tan.c | 338 |
1 files changed, 164 insertions, 174 deletions
diff --git a/arch/x86/math-emu/poly_tan.c b/arch/x86/math-emu/poly_tan.c index 8df3e03..c0d181e 100644 --- a/arch/x86/math-emu/poly_tan.c +++ b/arch/x86/math-emu/poly_tan.c @@ -17,206 +17,196 @@ #include "control_w.h" #include "poly.h" - #define HiPOWERop 3 /* odd poly, positive terms */ -static const unsigned long long oddplterm[HiPOWERop] = -{ - 0x0000000000000000LL, - 0x0051a1cf08fca228LL, - 0x0000000071284ff7LL +static const unsigned long long oddplterm[HiPOWERop] = { + 0x0000000000000000LL, + 0x0051a1cf08fca228LL, + 0x0000000071284ff7LL }; #define HiPOWERon 2 /* odd poly, negative terms */ -static const unsigned long long oddnegterm[HiPOWERon] = -{ - 0x1291a9a184244e80LL, - 0x0000583245819c21LL +static const unsigned long long oddnegterm[HiPOWERon] = { + 0x1291a9a184244e80LL, + 0x0000583245819c21LL }; #define HiPOWERep 2 /* even poly, positive terms */ -static const unsigned long long evenplterm[HiPOWERep] = -{ - 0x0e848884b539e888LL, - 0x00003c7f18b887daLL +static const unsigned long long evenplterm[HiPOWERep] = { + 0x0e848884b539e888LL, + 0x00003c7f18b887daLL }; #define HiPOWERen 2 /* even poly, negative terms */ -static const unsigned long long evennegterm[HiPOWERen] = -{ - 0xf1f0200fd51569ccLL, - 0x003afb46105c4432LL +static const unsigned long long evennegterm[HiPOWERen] = { + 0xf1f0200fd51569ccLL, + 0x003afb46105c4432LL }; static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL; - /*--- poly_tan() ------------------------------------------------------------+ | | +---------------------------------------------------------------------------*/ -void poly_tan(FPU_REG *st0_ptr) +void poly_tan(FPU_REG * st0_ptr) { - long int exponent; - int invert; - Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, - argSignif, fix_up; - unsigned long adj; + long int exponent; + int invert; + Xsig argSq, argSqSq, accumulatoro, accumulatore, accum, + argSignif, fix_up; + unsigned long adj; - exponent = exponent(st0_ptr); + exponent = exponent(st0_ptr); #ifdef PARANOID - if ( signnegative(st0_ptr) ) /* Can't hack a number < 0.0 */ - { arith_invalid(0); return; } /* Need a positive number */ + if (signnegative(st0_ptr)) { /* Can't hack a number < 0.0 */ + arith_invalid(0); + return; + } /* Need a positive number */ #endif /* PARANOID */ - /* Split the problem into two domains, smaller and larger than pi/4 */ - if ( (exponent == 0) || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2)) ) - { - /* The argument is greater than (approx) pi/4 */ - invert = 1; - accum.lsw = 0; - XSIG_LL(accum) = significand(st0_ptr); - - if ( exponent == 0 ) - { - /* The argument is >= 1.0 */ - /* Put the binary point at the left. */ - XSIG_LL(accum) <<= 1; - } - /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */ - XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum); - /* This is a special case which arises due to rounding. */ - if ( XSIG_LL(accum) == 0xffffffffffffffffLL ) - { - FPU_settag0(TAG_Valid); - significand(st0_ptr) = 0x8a51e04daabda360LL; - setexponent16(st0_ptr, (0x41 + EXTENDED_Ebias) | SIGN_Negative); - return; + /* Split the problem into two domains, smaller and larger than pi/4 */ + if ((exponent == 0) + || ((exponent == -1) && (st0_ptr->sigh > 0xc90fdaa2))) { + /* The argument is greater than (approx) pi/4 */ + invert = 1; + accum.lsw = 0; + XSIG_LL(accum) = significand(st0_ptr); + + if (exponent == 0) { + /* The argument is >= 1.0 */ + /* Put the binary point at the left. */ + XSIG_LL(accum) <<= 1; + } + /* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */ + XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum); + /* This is a special case which arises due to rounding. */ + if (XSIG_LL(accum) == 0xffffffffffffffffLL) { + FPU_settag0(TAG_Valid); + significand(st0_ptr) = 0x8a51e04daabda360LL; + setexponent16(st0_ptr, + (0x41 + EXTENDED_Ebias) | SIGN_Negative); + return; + } + + argSignif.lsw = accum.lsw; + XSIG_LL(argSignif) = XSIG_LL(accum); + exponent = -1 + norm_Xsig(&argSignif); + } else { + invert = 0; + argSignif.lsw = 0; + XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr); + + if (exponent < -1) { + /* shift the argument right by the required places */ + if (FPU_shrx(&XSIG_LL(accum), -1 - exponent) >= + 0x80000000U) + XSIG_LL(accum)++; /* round up */ + } } - argSignif.lsw = accum.lsw; - XSIG_LL(argSignif) = XSIG_LL(accum); - exponent = -1 + norm_Xsig(&argSignif); - } - else - { - invert = 0; - argSignif.lsw = 0; - XSIG_LL(accum) = XSIG_LL(argSignif) = significand(st0_ptr); - - if ( exponent < -1 ) - { - /* shift the argument right by the required places */ - if ( FPU_shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U ) - XSIG_LL(accum) ++; /* round up */ - } - } - - XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw; - mul_Xsig_Xsig(&argSq, &argSq); - XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw; - mul_Xsig_Xsig(&argSqSq, &argSqSq); - - /* Compute the negative terms for the numerator polynomial */ - accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; - polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1); - mul_Xsig_Xsig(&accumulatoro, &argSq); - negate_Xsig(&accumulatoro); - /* Add the positive terms */ - polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1); - - - /* Compute the positive terms for the denominator polynomial */ - accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0; - polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1); - mul_Xsig_Xsig(&accumulatore, &argSq); - negate_Xsig(&accumulatore); - /* Add the negative terms */ - polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1); - /* Multiply by arg^2 */ - mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); - mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); - /* de-normalize and divide by 2 */ - shr_Xsig(&accumulatore, -2*(1+exponent) + 1); - negate_Xsig(&accumulatore); /* This does 1 - accumulator */ - - /* Now find the ratio. */ - if ( accumulatore.msw == 0 ) - { - /* accumulatoro must contain 1.0 here, (actually, 0) but it - really doesn't matter what value we use because it will - have negligible effect in later calculations - */ - XSIG_LL(accum) = 0x8000000000000000LL; - accum.lsw = 0; - } - else - { - div_Xsig(&accumulatoro, &accumulatore, &accum); - } - - /* Multiply by 1/3 * arg^3 */ - mul64_Xsig(&accum, &XSIG_LL(argSignif)); - mul64_Xsig(&accum, &XSIG_LL(argSignif)); - mul64_Xsig(&accum, &XSIG_LL(argSignif)); - mul64_Xsig(&accum, &twothirds); - shr_Xsig(&accum, -2*(exponent+1)); - - /* tan(arg) = arg + accum */ - add_two_Xsig(&accum, &argSignif, &exponent); - - if ( invert ) - { - /* We now have the value of tan(pi_2 - arg) where pi_2 is an - approximation for pi/2 - */ - /* The next step is to fix the answer to compensate for the - error due to the approximation used for pi/2 - */ - - /* This is (approx) delta, the error in our approx for pi/2 - (see above). It has an exponent of -65 - */ - XSIG_LL(fix_up) = 0x898cc51701b839a2LL; - fix_up.lsw = 0; - - if ( exponent == 0 ) - adj = 0xffffffff; /* We want approx 1.0 here, but - this is close enough. */ - else if ( exponent > -30 ) - { - adj = accum.msw >> -(exponent+1); /* tan */ - adj = mul_32_32(adj, adj); /* tan^2 */ + XSIG_LL(argSq) = XSIG_LL(accum); + argSq.lsw = accum.lsw; + mul_Xsig_Xsig(&argSq, &argSq); + XSIG_LL(argSqSq) = XSIG_LL(argSq); + argSqSq.lsw = argSq.lsw; + mul_Xsig_Xsig(&argSqSq, &argSqSq); + + /* Compute the negative terms for the numerator polynomial */ + accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0; + polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, + HiPOWERon - 1); + mul_Xsig_Xsig(&accumulatoro, &argSq); + negate_Xsig(&accumulatoro); + /* Add the positive terms */ + polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, + HiPOWERop - 1); + + /* Compute the positive terms for the denominator polynomial */ + accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0; + polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, + HiPOWERep - 1); + mul_Xsig_Xsig(&accumulatore, &argSq); + negate_Xsig(&accumulatore); + /* Add the negative terms */ + polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, + HiPOWERen - 1); + /* Multiply by arg^2 */ + mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); + mul64_Xsig(&accumulatore, &XSIG_LL(argSignif)); + /* de-normalize and divide by 2 */ + shr_Xsig(&accumulatore, -2 * (1 + exponent) + 1); + negate_Xsig(&accumulatore); /* This does 1 - accumulator */ + + /* Now find the ratio. */ + if (accumulatore.msw == 0) { + /* accumulatoro must contain 1.0 here, (actually, 0) but it + really doesn't matter what value we use because it will + have negligible effect in later calculations + */ + XSIG_LL(accum) = 0x8000000000000000LL; + accum.lsw = 0; + } else { + div_Xsig(&accumulatoro, &accumulatore, &accum); } - else - adj = 0; - adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */ - - fix_up.msw += adj; - if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */ - { - /* Yes, we need to add an msb */ - shr_Xsig(&fix_up, 1); - fix_up.msw |= 0x80000000; - shr_Xsig(&fix_up, 64 + exponent); + + /* Multiply by 1/3 * arg^3 */ + mul64_Xsig(&accum, &XSIG_LL(argSignif)); + mul64_Xsig(&accum, &XSIG_LL(argSignif)); + mul64_Xsig(&accum, &XSIG_LL(argSignif)); + mul64_Xsig(&accum, &twothirds); + shr_Xsig(&accum, -2 * (exponent + 1)); + + /* tan(arg) = arg + accum */ + add_two_Xsig(&accum, &argSignif, &exponent); + + if (invert) { + /* We now have the value of tan(pi_2 - arg) where pi_2 is an + approximation for pi/2 + */ + /* The next step is to fix the answer to compensate for the + error due to the approximation used for pi/2 + */ + + /* This is (approx) delta, the error in our approx for pi/2 + (see above). It has an exponent of -65 + */ + XSIG_LL(fix_up) = 0x898cc51701b839a2LL; + fix_up.lsw = 0; + + if (exponent == 0) + adj = 0xffffffff; /* We want approx 1.0 here, but + this is close enough. */ + else if (exponent > -30) { + adj = accum.msw >> -(exponent + 1); /* tan */ + adj = mul_32_32(adj, adj); /* tan^2 */ + } else + adj = 0; + adj = mul_32_32(0x898cc517, adj); /* delta * tan^2 */ + + fix_up.msw += adj; + if (!(fix_up.msw & 0x80000000)) { /* did fix_up overflow ? */ + /* Yes, we need to add an msb */ + shr_Xsig(&fix_up, 1); + fix_up.msw |= 0x80000000; + shr_Xsig(&fix_up, 64 + exponent); + } else + shr_Xsig(&fix_up, 65 + exponent); + + add_two_Xsig(&accum, &fix_up, &exponent); + + /* accum now contains tan(pi/2 - arg). + Use tan(arg) = 1.0 / tan(pi/2 - arg) + */ + accumulatoro.lsw = accumulatoro.midw = 0; + accumulatoro.msw = 0x80000000; + div_Xsig(&accumulatoro, &accum, &accum); + exponent = -exponent - 1; } - else - shr_Xsig(&fix_up, 65 + exponent); - - add_two_Xsig(&accum, &fix_up, &exponent); - - /* accum now contains tan(pi/2 - arg). - Use tan(arg) = 1.0 / tan(pi/2 - arg) - */ - accumulatoro.lsw = accumulatoro.midw = 0; - accumulatoro.msw = 0x80000000; - div_Xsig(&accumulatoro, &accum, &accum); - exponent = - exponent - 1; - } - - /* Transfer the result */ - round_Xsig(&accum); - FPU_settag0(TAG_Valid); - significand(st0_ptr) = XSIG_LL(accum); - setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */ + + /* Transfer the result */ + round_Xsig(&accum); + FPU_settag0(TAG_Valid); + significand(st0_ptr) = XSIG_LL(accum); + setexponent16(st0_ptr, exponent + EXTENDED_Ebias); /* Result is positive. */ } |