/*---------------------------------------------------------------------------+
| fpu_trig.c |
| |
| Implementation of the FPU "transcendental" functions. |
| |
| Copyright (C) 1992,1993,1994,1997,1999 |
| W. Metzenthen, 22 Parker St, Ormond, Vic 3163, |
| Australia. E-mail billm@melbpc.org.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
/* 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 *st0_ptr, int even)
{
FPU_REG tmp;
u_char tmptag;
unsigned long long q;
int old_cw = control_word, saved_status = partial_status;
int tag, st0_tag = TAG_Valid;
if (exponent(st0_ptr) >= 63) {
partial_status |= SW_C2; /* Reduction incomplete. */
return -1;
}
control_word &= ~CW_RC;
control_word |= RC_CHOP;
setpositive(st0_ptr);
tag = FPU_u_div(st0_ptr, &CONST_PI2, &tmp, PR_64_BITS | RC_CHOP | 0x3f,
SIGN_POS);
FPU_round_to_int(&tmp, tag); /* Fortunately, this can't overflow
to 2^64 */
q = significand(&tmp);
if (q) {
rem_kernel(significand(st0_ptr),
&significand(&tmp),
significand(&CONST_PI2),
q, exponent(st0_ptr) - exponent(&CONST_PI2));
setexponent16(&tmp, exponent(&CONST_PI2));
st0_tag = FPU_normalize(&tmp);
FPU_copy_to_reg0(&tmp, st0_tag);
}
if ((even && !(q & 1)) || (!even && (q & 1))) {
st0_tag =
FPU_sub(REV | LOADED | TAG_Valid, (int)&CONST_PI2,
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 ((exponent(st0_ptr) <= exponent(&CONST_PI2extra) + 64)
|| (q > 1)) {
/* This code gives the effect of having pi/2 to better than
128 bits precision. */
significand(&tmp) = q + 1;
setexponent16(&tmp, 63);
FPU_normalize(&tmp);
tmptag =
FPU_u_mul(&CONST_PI2extra, &tmp, &tmp,
FULL_PRECISION, SIGN_POS,
exponent(&CONST_PI2extra) +
e
| 