📄 fpu_trig.c
字号:
/*---------------------------------------------------------------------------+
| fpu_trig.c |
| |
| Implementation of the FPU "transcendental" functions. |
| |
| Copyright (C) 1992,1993,1994 |
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
| Australia. E-mail billm@vaxc.cc.monash.edu.au |
| |
| |
+---------------------------------------------------------------------------*/
#include "fpu_system.h"
#include "exception.h"
#include "fpu_emu.h"
#include "status_w.h"
#include "control_w.h"
#include "reg_constant.h"
static void rem_kernel(unsigned long long st0, unsigned long long *y,
unsigned long long st1,
unsigned long long q, int n);
#define BETTER_THAN_486
#define FCOS 4
#define FPTAN 1
/* Used only by fptan, fsin, fcos, and fsincos. */
/* This routine produces very accurate results, similar to
using a value of pi with more than 128 bits precision. */
/* Limited measurements show no results worse than 64 bit precision
except for the results for arguments close to 2^63, where the
precision of the result sometimes degrades to about 63.9 bits */
static int trig_arg(FPU_REG *X, int even)
{
FPU_REG tmp;
unsigned long long q;
int old_cw = control_word, saved_status = partial_status;
if ( X->exp >= EXP_BIAS + 63 )
{
partial_status |= SW_C2; /* Reduction incomplete. */
return -1;
}
control_word &= ~CW_RC;
control_word |= RC_CHOP;
reg_div(X, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f);
round_to_int(&tmp); /* Fortunately, this can't overflow
to 2^64 */
q = significand(&tmp);
if ( q )
{
rem_kernel(significand(X),
&significand(&tmp),
significand(&CONST_PI2),
q, X->exp - CONST_PI2.exp);
tmp.exp = CONST_PI2.exp;
normalize(&tmp);
reg_move(&tmp, X);
}
if ( even == FPTAN )
{
if ( ((X->exp >= EXP_BIAS) ||
((X->exp == EXP_BIAS-1)
&& (X->sigh >= 0xc90fdaa2))) ^ (q & 1) )
even = FCOS;
else
even = 0;
}
if ( (even && !(q & 1)) || (!even && (q & 1)) )
{
reg_sub(&CONST_PI2, X, X, FULL_PRECISION);
#ifdef BETTER_THAN_486
/* So far, the results are exact but based upon a 64 bit
precision approximation to pi/2. The technique used
now is equivalent to using an approximation to pi/2 which
is accurate to about 128 bits. */
if ( (X->exp <= CONST_PI2extra.exp + 64) || (q > 1) )
{
/* This code gives the effect of having p/2 to better than
128 bits precision. */
significand(&tmp) = q + 1;
tmp.exp = EXP_BIAS + 63;
tmp.tag = TW_Valid;
normalize(&tmp);
reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
reg_add(X, &tmp, X, FULL_PRECISION);
if ( X->sign == SIGN_NEG )
{
/* CONST_PI2extra is negative, so the result of the addition
can be negative. This means that the argument is actually
in a different quadrant. The correction is always < pi/2,
so it can't overflow into yet another quadrant. */
X->sign = SIGN_POS;
q++;
}
}
#endif BETTER_THAN_486
}
#ifdef BETTER_THAN_486
else
{
/* So far, the results are exact but based upon a 64 bit
precision approximation to pi/2. The technique used
now is equivalent to using an approximation to pi/2 which
is accurate to about 128 bits. */
if ( ((q > 0) && (X->exp <= CONST_PI2extra.exp + 64)) || (q > 1) )
{
/* This code gives the effect of having p/2 to better than
128 bits precision. */
significand(&tmp) = q;
tmp.exp = EXP_BIAS + 63;
tmp.tag = TW_Valid;
normalize(&tmp);
reg_mul(&CONST_PI2extra, &tmp, &tmp, FULL_PRECISION);
reg_sub(X, &tmp, X, FULL_PRECISION);
if ( (X->exp == CONST_PI2.exp) &&
((X->sigh > CONST_PI2.sigh)
|| ((X->sigh == CONST_PI2.sigh)
&& (X->sigl > CONST_PI2.sigl))) )
{
/* CONST_PI2extra is negative, so the result of the
subtraction can be larger than pi/2. This means
that the argument is actually in a different quadrant.
The correction is always < pi/2, so it can't overflow
into yet another quadrant. */
reg_sub(&CONST_PI, X, X, FULL_PRECISION);
q++;
}
}
}
#endif BETTER_THAN_486
control_word = old_cw;
partial_status = saved_status & ~SW_C2; /* Reduction complete. */
return (q & 3) | even;
}
/* Convert a long to register */
void convert_l2reg(long const *arg, FPU_REG *dest)
{
long num = *arg;
if (num == 0)
{ reg_move(&CONST_Z, dest); return; }
if (num > 0)
dest->sign = SIGN_POS;
else
{ num = -num; dest->sign = SIGN_NEG; }
dest->sigh = num;
dest->sigl = 0;
dest->exp = EXP_BIAS + 31;
dest->tag = TW_Valid;
normalize(dest);
}
static void single_arg_error(void)
{
switch ( FPU_st0_tag )
{
case TW_NaN:
if ( !(FPU_st0_ptr->sigh & 0x40000000) ) /* Signaling ? */
{
EXCEPTION(EX_Invalid);
if ( control_word & CW_Invalid )
FPU_st0_ptr->sigh |= 0x40000000; /* Convert to a QNaN */
}
break; /* return with a NaN in st(0) */
case TW_Empty:
stack_underflow(); /* Puts a QNaN in st(0) */
break;
#ifdef PARANOID
default:
EXCEPTION(EX_INTERNAL|0x0112);
#endif PARANOID
}
}
static void single_arg_2_error(void)
{
FPU_REG *st_new_ptr;
switch ( FPU_st0_tag )
{
case TW_NaN:
if ( !(FPU_st0_ptr->sigh & 0x40000000) ) /* Signaling ? */
{
EXCEPTION(EX_Invalid);
if ( control_word & CW_Invalid )
{
/* The masked response */
/* Convert to a QNaN */
FPU_st0_ptr->sigh |= 0x40000000;
st_new_ptr = &st(-1);
push();
reg_move(&st(1), FPU_st0_ptr);
}
}
else
{
/* A QNaN */
st_new_ptr = &st(-1);
push();
reg_move(&st(1), FPU_st0_ptr);
}
break; /* return with a NaN in st(0) */
#ifdef PARANOID
default:
EXCEPTION(EX_INTERNAL|0x0112);
#endif PARANOID
}
}
/*---------------------------------------------------------------------------*/
static void f2xm1(void)
{
clear_C1();
switch ( FPU_st0_tag )
{
case TW_Valid:
{
FPU_REG rv, tmp;
if ( FPU_st0_ptr->exp >= 0 )
{
/* For an 80486 FPU, the result is undefined. */
}
else if ( FPU_st0_ptr->exp >= -64 )
{
if ( FPU_st0_ptr->sign == SIGN_POS )
{
/* poly_2xm1(x) requires 0 < x < 1. */
poly_2xm1(FPU_st0_ptr, &rv);
reg_mul(&rv, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
}
else
{
/* poly_2xm1(x) doesn't handle negative numbers yet. */
/* So we compute z=poly_2xm1(-x), and the answer is
then -z/(1+z) */
FPU_st0_ptr->sign = SIGN_POS;
poly_2xm1(FPU_st0_ptr, &rv);
reg_mul(&rv, FPU_st0_ptr, &rv, FULL_PRECISION);
reg_add(&rv, &CONST_1, &tmp, FULL_PRECISION);
reg_div(&rv, &tmp, FPU_st0_ptr, FULL_PRECISION);
FPU_st0_ptr->sign = SIGN_NEG;
}
}
else
{
#ifdef DENORM_OPERAND
if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
return;
#endif DENORM_OPERAND
/* For very small arguments, this is accurate enough. */
reg_mul(&CONST_LN2, FPU_st0_ptr, FPU_st0_ptr, FULL_PRECISION);
}
set_precision_flag_up();
return;
}
case TW_Zero:
return;
case TW_Infinity:
if ( FPU_st0_ptr->sign == SIGN_NEG )
{
/* -infinity gives -1 (p16-10) */
reg_move(&CONST_1, FPU_st0_ptr);
FPU_st0_ptr->sign = SIGN_NEG;
}
return;
default:
single_arg_error();
}
}
static void fptan(void)
{
FPU_REG *st_new_ptr;
int q;
char arg_sign = FPU_st0_ptr->sign;
/* Stack underflow has higher priority */
if ( FPU_st0_tag == TW_Empty )
{
stack_underflow(); /* Puts a QNaN in st(0) */
if ( control_word & CW_Invalid )
{
st_new_ptr = &st(-1);
push();
stack_underflow(); /* Puts a QNaN in the new st(0) */
}
return;
}
if ( STACK_OVERFLOW )
{ stack_overflow(); return; }
switch ( FPU_st0_tag )
{
case TW_Valid:
if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
{
FPU_st0_ptr->sign = SIGN_POS;
if ( (q = trig_arg(FPU_st0_ptr, FPTAN)) != -1 )
{
reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr,
FULL_PRECISION);
poly_tan(FPU_st0_ptr, FPU_st0_ptr, q & FCOS);
FPU_st0_ptr->sign = (q & 1) ^ arg_sign;
}
else
{
/* Operand is out of range */
FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
return;
}
}
else
{
/* For a small arg, the result == the argument */
/* Underflow may happen */
if ( FPU_st0_ptr->exp <= EXP_UNDER )
{
#ifdef DENORM_OPERAND
if ( denormal_operand() )
return;
#endif DENORM_OPERAND
/* A denormal result has been produced.
Precision must have been lost, this is always
an underflow. */
if ( arith_underflow(FPU_st0_ptr) )
return;
}
else
set_precision_flag_up(); /* Must be up. */
}
push();
reg_move(&CONST_1, FPU_st0_ptr);
return;
break;
case TW_Infinity:
/* The 80486 treats infinity as an invalid operand */
arith_invalid(FPU_st0_ptr);
if ( control_word & CW_Invalid )
{
st_new_ptr = &st(-1);
push();
arith_invalid(FPU_st0_ptr);
}
return;
case TW_Zero:
push();
reg_move(&CONST_1, FPU_st0_ptr);
setcc(0);
break;
default:
single_arg_2_error();
break;
}
}
static void fxtract(void)
{
FPU_REG *st_new_ptr;
register FPU_REG *st1_ptr = FPU_st0_ptr; /* anticipate */
if ( STACK_OVERFLOW )
{ stack_overflow(); return; }
clear_C1();
if ( !(FPU_st0_tag ^ TW_Valid) )
{
long e;
#ifdef DENORM_OPERAND
if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
return;
#endif DENORM_OPERAND
push();
reg_move(st1_ptr, FPU_st0_ptr);
FPU_st0_ptr->exp = EXP_BIAS;
e = st1_ptr->exp - EXP_BIAS;
convert_l2reg(&e, st1_ptr);
return;
}
else if ( FPU_st0_tag == TW_Zero )
{
char sign = FPU_st0_ptr->sign;
if ( divide_by_zero(SIGN_NEG, FPU_st0_ptr) )
return;
push();
reg_move(&CONST_Z, FPU_st0_ptr);
FPU_st0_ptr->sign = sign;
return;
}
else if ( FPU_st0_tag == TW_Infinity )
{
char sign = FPU_st0_ptr->sign;
FPU_st0_ptr->sign = SIGN_POS;
push();
reg_move(&CONST_INF, FPU_st0_ptr);
FPU_st0_ptr->sign = sign;
return;
}
else if ( FPU_st0_tag == TW_NaN )
{
if ( real_2op_NaN(FPU_st0_ptr, FPU_st0_ptr, FPU_st0_ptr) )
return;
push();
reg_move(st1_ptr, FPU_st0_ptr);
return;
}
else if ( FPU_st0_tag == TW_Empty )
{
/* Is this the correct behaviour? */
if ( control_word & EX_Invalid )
{
stack_underflow();
push();
stack_underflow();
}
else
EXCEPTION(EX_StackUnder);
}
#ifdef PARANOID
else
EXCEPTION(EX_INTERNAL | 0x119);
#endif PARANOID
}
static void fdecstp(void)
{
clear_C1();
top--; /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
}
static void fincstp(void)
{
clear_C1();
top++; /* FPU_st0_ptr will be fixed in math_emulate() before the next instr */
}
static void fsqrt_(void)
{
clear_C1();
if ( !(FPU_st0_tag ^ TW_Valid) )
{
int expon;
if (FPU_st0_ptr->sign == SIGN_NEG)
{
arith_invalid(FPU_st0_ptr); /* sqrt(negative) is invalid */
return;
}
#ifdef DENORM_OPERAND
if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
return;
#endif DENORM_OPERAND
expon = FPU_st0_ptr->exp - EXP_BIAS;
FPU_st0_ptr->exp = EXP_BIAS + (expon & 1); /* make st(0) in [1.0 .. 4.0) */
wm_sqrt(FPU_st0_ptr, control_word); /* Do the computation */
FPU_st0_ptr->exp += expon >> 1;
FPU_st0_ptr->sign = SIGN_POS;
}
else if ( FPU_st0_tag == TW_Zero )
return;
else if ( FPU_st0_tag == TW_Infinity )
{
if ( FPU_st0_ptr->sign == SIGN_NEG )
arith_invalid(FPU_st0_ptr); /* sqrt(-Infinity) is invalid */
return;
}
else
{ single_arg_error(); return; }
}
static void frndint_(void)
{
int flags;
if ( !(FPU_st0_tag ^ TW_Valid) )
{
if (FPU_st0_ptr->exp > EXP_BIAS+63)
return;
#ifdef DENORM_OPERAND
if ( (FPU_st0_ptr->exp <= EXP_UNDER) && (denormal_operand()) )
return;
#endif DENORM_OPERAND
/* Fortunately, this can't overflow to 2^64 */
if ( (flags = round_to_int(FPU_st0_ptr)) )
set_precision_flag(flags);
FPU_st0_ptr->exp = EXP_BIAS + 63;
normalize(FPU_st0_ptr);
return;
}
else if ( (FPU_st0_tag == TW_Zero) || (FPU_st0_tag == TW_Infinity) )
return;
else
single_arg_error();
}
static void fsin(void)
{
char arg_sign = FPU_st0_ptr->sign;
if ( FPU_st0_tag == TW_Valid )
{
FPU_REG rv;
int q;
if ( FPU_st0_ptr->exp > EXP_BIAS - 40 )
{
FPU_st0_ptr->sign = SIGN_POS;
if ( (q = trig_arg(FPU_st0_ptr, 0)) != -1 )
{
reg_div(FPU_st0_ptr, &CONST_PI2, FPU_st0_ptr, FULL_PRECISION);
poly_sine(FPU_st0_ptr, &rv);
if (q & 2)
rv.sign ^= SIGN_POS ^ SIGN_NEG;
rv.sign ^= arg_sign;
reg_move(&rv, FPU_st0_ptr);
/* We do not really know if up or down */
set_precision_flag_up();
return;
}
else
{
/* Operand is out of range */
FPU_st0_ptr->sign = arg_sign; /* restore st(0) */
return;
}
}
else
{
/* For a small arg, the result == the argument */
/* Underflow may happen */
if ( FPU_st0_ptr->exp <= EXP_UNDER )
{
#ifdef DENORM_OPERAND
if ( denormal_operand() )
return;
#endif DENORM_OPERAND
/* A denormal result has been produced.
Precision must have been lost, this is always
an underflow. */
arith_underflow(FPU_st0_ptr);
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -