/*---------------------------------------------------------------------------+
|
/*---------------------------------------------------------------------------+
|
| poly_tan.c |
|
| poly_tan.c |
|
| |
|
| |
|
| Compute the tan of a FPU_REG, using a polynomial approximation. |
|
| Compute the tan of a FPU_REG, using a polynomial approximation. |
|
| |
|
| |
|
| Copyright (C) 1992,1993,1994 |
|
| Copyright (C) 1992,1993,1994 |
|
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
|
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
|
| Australia. E-mail billm@vaxc.cc.monash.edu.au |
|
| Australia. E-mail billm@vaxc.cc.monash.edu.au |
|
| |
|
| |
|
| |
|
| |
|
+---------------------------------------------------------------------------*/
|
+---------------------------------------------------------------------------*/
|
|
|
#include "exception.h"
|
#include "exception.h"
|
#include "reg_constant.h"
|
#include "reg_constant.h"
|
#include "fpu_emu.h"
|
#include "fpu_emu.h"
|
#include "control_w.h"
|
#include "control_w.h"
|
#include "poly.h"
|
#include "poly.h"
|
|
|
|
|
#define HiPOWERop 3 /* odd poly, positive terms */
|
#define HiPOWERop 3 /* odd poly, positive terms */
|
static const unsigned long long oddplterm[HiPOWERop] =
|
static const unsigned long long oddplterm[HiPOWERop] =
|
{
|
{
|
0x0000000000000000LL,
|
0x0000000000000000LL,
|
0x0051a1cf08fca228LL,
|
0x0051a1cf08fca228LL,
|
0x0000000071284ff7LL
|
0x0000000071284ff7LL
|
};
|
};
|
|
|
#define HiPOWERon 2 /* odd poly, negative terms */
|
#define HiPOWERon 2 /* odd poly, negative terms */
|
static const unsigned long long oddnegterm[HiPOWERon] =
|
static const unsigned long long oddnegterm[HiPOWERon] =
|
{
|
{
|
0x1291a9a184244e80LL,
|
0x1291a9a184244e80LL,
|
0x0000583245819c21LL
|
0x0000583245819c21LL
|
};
|
};
|
|
|
#define HiPOWERep 2 /* even poly, positive terms */
|
#define HiPOWERep 2 /* even poly, positive terms */
|
static const unsigned long long evenplterm[HiPOWERep] =
|
static const unsigned long long evenplterm[HiPOWERep] =
|
{
|
{
|
0x0e848884b539e888LL,
|
0x0e848884b539e888LL,
|
0x00003c7f18b887daLL
|
0x00003c7f18b887daLL
|
};
|
};
|
|
|
#define HiPOWERen 2 /* even poly, negative terms */
|
#define HiPOWERen 2 /* even poly, negative terms */
|
static const unsigned long long evennegterm[HiPOWERen] =
|
static const unsigned long long evennegterm[HiPOWERen] =
|
{
|
{
|
0xf1f0200fd51569ccLL,
|
0xf1f0200fd51569ccLL,
|
0x003afb46105c4432LL
|
0x003afb46105c4432LL
|
};
|
};
|
|
|
static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
|
static const unsigned long long twothirds = 0xaaaaaaaaaaaaaaabLL;
|
|
|
|
|
/*--- poly_tan() ------------------------------------------------------------+
|
/*--- poly_tan() ------------------------------------------------------------+
|
| |
|
| |
|
+---------------------------------------------------------------------------*/
|
+---------------------------------------------------------------------------*/
|
void poly_tan(FPU_REG const *arg, FPU_REG *result)
|
void poly_tan(FPU_REG const *arg, FPU_REG *result)
|
{
|
{
|
long int exponent;
|
long int exponent;
|
int invert;
|
int invert;
|
Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
|
Xsig argSq, argSqSq, accumulatoro, accumulatore, accum,
|
argSignif, fix_up;
|
argSignif, fix_up;
|
unsigned long adj;
|
unsigned long adj;
|
|
|
exponent = arg->exp - EXP_BIAS;
|
exponent = arg->exp - EXP_BIAS;
|
|
|
#ifdef PARANOID
|
#ifdef PARANOID
|
if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
|
if ( arg->sign != 0 ) /* Can't hack a number < 0.0 */
|
{ arith_invalid(result); return; } /* Need a positive number */
|
{ arith_invalid(result); return; } /* Need a positive number */
|
#endif PARANOID
|
#endif PARANOID
|
|
|
/* Split the problem into two domains, smaller and larger than pi/4 */
|
/* Split the problem into two domains, smaller and larger than pi/4 */
|
if ( (exponent == 0) || ((exponent == -1) && (arg->sigh > 0xc90fdaa2)) )
|
if ( (exponent == 0) || ((exponent == -1) && (arg->sigh > 0xc90fdaa2)) )
|
{
|
{
|
/* The argument is greater than (approx) pi/4 */
|
/* The argument is greater than (approx) pi/4 */
|
invert = 1;
|
invert = 1;
|
accum.lsw = 0;
|
accum.lsw = 0;
|
XSIG_LL(accum) = significand(arg);
|
XSIG_LL(accum) = significand(arg);
|
|
|
if ( exponent == 0 )
|
if ( exponent == 0 )
|
{
|
{
|
/* The argument is >= 1.0 */
|
/* The argument is >= 1.0 */
|
/* Put the binary point at the left. */
|
/* Put the binary point at the left. */
|
XSIG_LL(accum) <<= 1;
|
XSIG_LL(accum) <<= 1;
|
}
|
}
|
/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
|
/* pi/2 in hex is: 1.921fb54442d18469 898CC51701B839A2 52049C1 */
|
XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
|
XSIG_LL(accum) = 0x921fb54442d18469LL - XSIG_LL(accum);
|
|
|
argSignif.lsw = accum.lsw;
|
argSignif.lsw = accum.lsw;
|
XSIG_LL(argSignif) = XSIG_LL(accum);
|
XSIG_LL(argSignif) = XSIG_LL(accum);
|
exponent = -1 + norm_Xsig(&argSignif);
|
exponent = -1 + norm_Xsig(&argSignif);
|
}
|
}
|
else
|
else
|
{
|
{
|
invert = 0;
|
invert = 0;
|
argSignif.lsw = 0;
|
argSignif.lsw = 0;
|
XSIG_LL(accum) = XSIG_LL(argSignif) = significand(arg);
|
XSIG_LL(accum) = XSIG_LL(argSignif) = significand(arg);
|
|
|
if ( exponent < -1 )
|
if ( exponent < -1 )
|
{
|
{
|
/* shift the argument right by the required places */
|
/* shift the argument right by the required places */
|
if ( shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
|
if ( shrx(&XSIG_LL(accum), -1-exponent) >= 0x80000000U )
|
XSIG_LL(accum) ++; /* round up */
|
XSIG_LL(accum) ++; /* round up */
|
}
|
}
|
}
|
}
|
|
|
XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
|
XSIG_LL(argSq) = XSIG_LL(accum); argSq.lsw = accum.lsw;
|
mul_Xsig_Xsig(&argSq, &argSq);
|
mul_Xsig_Xsig(&argSq, &argSq);
|
XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
|
XSIG_LL(argSqSq) = XSIG_LL(argSq); argSqSq.lsw = argSq.lsw;
|
mul_Xsig_Xsig(&argSqSq, &argSqSq);
|
mul_Xsig_Xsig(&argSqSq, &argSqSq);
|
|
|
/* Compute the negative terms for the numerator polynomial */
|
/* Compute the negative terms for the numerator polynomial */
|
accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
|
accumulatoro.msw = accumulatoro.midw = accumulatoro.lsw = 0;
|
polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
|
polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddnegterm, HiPOWERon-1);
|
mul_Xsig_Xsig(&accumulatoro, &argSq);
|
mul_Xsig_Xsig(&accumulatoro, &argSq);
|
negate_Xsig(&accumulatoro);
|
negate_Xsig(&accumulatoro);
|
/* Add the positive terms */
|
/* Add the positive terms */
|
polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
|
polynomial_Xsig(&accumulatoro, &XSIG_LL(argSqSq), oddplterm, HiPOWERop-1);
|
|
|
|
|
/* Compute the positive terms for the denominator polynomial */
|
/* Compute the positive terms for the denominator polynomial */
|
accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
|
accumulatore.msw = accumulatore.midw = accumulatore.lsw = 0;
|
polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
|
polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evenplterm, HiPOWERep-1);
|
mul_Xsig_Xsig(&accumulatore, &argSq);
|
mul_Xsig_Xsig(&accumulatore, &argSq);
|
negate_Xsig(&accumulatore);
|
negate_Xsig(&accumulatore);
|
/* Add the negative terms */
|
/* Add the negative terms */
|
polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
|
polynomial_Xsig(&accumulatore, &XSIG_LL(argSqSq), evennegterm, HiPOWERen-1);
|
/* Multiply by arg^2 */
|
/* Multiply by arg^2 */
|
mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
|
mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
|
mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
|
mul64_Xsig(&accumulatore, &XSIG_LL(argSignif));
|
/* de-normalize and divide by 2 */
|
/* de-normalize and divide by 2 */
|
shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
|
shr_Xsig(&accumulatore, -2*(1+exponent) + 1);
|
negate_Xsig(&accumulatore); /* This does 1 - accumulator */
|
negate_Xsig(&accumulatore); /* This does 1 - accumulator */
|
|
|
/* Now find the ratio. */
|
/* Now find the ratio. */
|
if ( accumulatore.msw == 0 )
|
if ( accumulatore.msw == 0 )
|
{
|
{
|
/* accumulatoro must contain 1.0 here, (actually, 0) but it
|
/* accumulatoro must contain 1.0 here, (actually, 0) but it
|
really doesn't matter what value we use because it will
|
really doesn't matter what value we use because it will
|
have negligible effect in later calculations
|
have negligible effect in later calculations
|
*/
|
*/
|
XSIG_LL(accum) = 0x8000000000000000LL;
|
XSIG_LL(accum) = 0x8000000000000000LL;
|
accum.lsw = 0;
|
accum.lsw = 0;
|
}
|
}
|
else
|
else
|
{
|
{
|
div_Xsig(&accumulatoro, &accumulatore, &accum);
|
div_Xsig(&accumulatoro, &accumulatore, &accum);
|
}
|
}
|
|
|
/* Multiply by 1/3 * arg^3 */
|
/* 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, &XSIG_LL(argSignif));
|
mul64_Xsig(&accum, &XSIG_LL(argSignif));
|
mul64_Xsig(&accum, &XSIG_LL(argSignif));
|
mul64_Xsig(&accum, &twothirds);
|
mul64_Xsig(&accum, &twothirds);
|
shr_Xsig(&accum, -2*(exponent+1));
|
shr_Xsig(&accum, -2*(exponent+1));
|
|
|
/* tan(arg) = arg + accum */
|
/* tan(arg) = arg + accum */
|
add_two_Xsig(&accum, &argSignif, &exponent);
|
add_two_Xsig(&accum, &argSignif, &exponent);
|
|
|
if ( invert )
|
if ( invert )
|
{
|
{
|
/* We now have the value of tan(pi_2 - arg) where pi_2 is an
|
/* We now have the value of tan(pi_2 - arg) where pi_2 is an
|
approximation for pi/2
|
approximation for pi/2
|
*/
|
*/
|
/* The next step is to fix the answer to compensate for the
|
/* The next step is to fix the answer to compensate for the
|
error due to the approximation used for pi/2
|
error due to the approximation used for pi/2
|
*/
|
*/
|
|
|
/* This is (approx) delta, the error in our approx for pi/2
|
/* This is (approx) delta, the error in our approx for pi/2
|
(see above). It has an exponent of -65
|
(see above). It has an exponent of -65
|
*/
|
*/
|
XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
|
XSIG_LL(fix_up) = 0x898cc51701b839a2LL;
|
fix_up.lsw = 0;
|
fix_up.lsw = 0;
|
|
|
if ( exponent == 0 )
|
if ( exponent == 0 )
|
adj = 0xffffffff; /* We want approx 1.0 here, but
|
adj = 0xffffffff; /* We want approx 1.0 here, but
|
this is close enough. */
|
this is close enough. */
|
else if ( exponent > -30 )
|
else if ( exponent > -30 )
|
{
|
{
|
adj = accum.msw >> -(exponent+1); /* tan */
|
adj = accum.msw >> -(exponent+1); /* tan */
|
mul_32_32(adj, adj, &adj); /* tan^2 */
|
mul_32_32(adj, adj, &adj); /* tan^2 */
|
}
|
}
|
else
|
else
|
adj = 0;
|
adj = 0;
|
mul_32_32(0x898cc517, adj, &adj); /* delta * tan^2 */
|
mul_32_32(0x898cc517, adj, &adj); /* delta * tan^2 */
|
|
|
fix_up.msw += adj;
|
fix_up.msw += adj;
|
if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */
|
if ( !(fix_up.msw & 0x80000000) ) /* did fix_up overflow ? */
|
{
|
{
|
/* Yes, we need to add an msb */
|
/* Yes, we need to add an msb */
|
shr_Xsig(&fix_up, 1);
|
shr_Xsig(&fix_up, 1);
|
fix_up.msw |= 0x80000000;
|
fix_up.msw |= 0x80000000;
|
shr_Xsig(&fix_up, 64 + exponent);
|
shr_Xsig(&fix_up, 64 + exponent);
|
}
|
}
|
else
|
else
|
shr_Xsig(&fix_up, 65 + exponent);
|
shr_Xsig(&fix_up, 65 + exponent);
|
|
|
add_two_Xsig(&accum, &fix_up, &exponent);
|
add_two_Xsig(&accum, &fix_up, &exponent);
|
|
|
/* accum now contains tan(pi/2 - arg).
|
/* accum now contains tan(pi/2 - arg).
|
Use tan(arg) = 1.0 / tan(pi/2 - arg)
|
Use tan(arg) = 1.0 / tan(pi/2 - arg)
|
*/
|
*/
|
accumulatoro.lsw = accumulatoro.midw = 0;
|
accumulatoro.lsw = accumulatoro.midw = 0;
|
accumulatoro.msw = 0x80000000;
|
accumulatoro.msw = 0x80000000;
|
div_Xsig(&accumulatoro, &accum, &accum);
|
div_Xsig(&accumulatoro, &accum, &accum);
|
exponent = - exponent - 1;
|
exponent = - exponent - 1;
|
}
|
}
|
|
|
/* Transfer the result */
|
/* Transfer the result */
|
round_Xsig(&accum);
|
round_Xsig(&accum);
|
*(short *)&(result->sign) = 0;
|
*(short *)&(result->sign) = 0;
|
significand(result) = XSIG_LL(accum);
|
significand(result) = XSIG_LL(accum);
|
result->exp = EXP_BIAS + exponent;
|
result->exp = EXP_BIAS + exponent;
|
|
|
}
|
}
|
|
|
